Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90	(revision 5539)
@@ -176,5 +176,5 @@
 
      !open output file
-     IF( lwp ) THEN
+     IF( lwm ) THEN
         CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
         CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
@@ -283,5 +283,5 @@
            DO jsec=1,nb_sec
 
-              IF( lwp )CALL dia_dct_wri(kt,jsec,secs(jsec))
+              IF( lwm )CALL dia_dct_wri(kt,jsec,secs(jsec))
             
               !nullify transports values after writing
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 5539)
@@ -166,16 +166,12 @@
          CASE ( 1 )                                  !  ORCA_R1 configurations
             !                                        ! =======================
-            ! This dirty section will be suppressed by simplification process: all this will come back in input files
-            ! Currently these hard-wired indices relate to the original (pre-v3.6) configuration 
-            ! which had a grid-size of 362x292.
-            ! This grid has been extended southwards for use with the under ice-shelf options (isf) introduced in v3.6. 
-            ! The original domain can still be used optionally if the isf code is not activated. 
-            ! An adjustment (isrow) is made to the hard-wired indices if the extended domain (362x332) is being used.
-            ! 
-            IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  ! Using pre-v3.6 files or adjusted start row from isf-extended grid
-            ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 40 ! Using full isf­extended domain. 
-            ENDIF                                       ! Adjust j­indices to account for more southerly starting latitude
+            ! This dirty section will be suppressed by simplification process:
+            ! all this will come back in input files
+            ! Currently these hard-wired indices relate to configuration with
+            ! extend grid (jpjglo=332)
+            isrow = 332 - jpjglo
+            !
             ii0 = 283           ;   ii1 = 283
-            ij0 = 201 + isrow   ;   ij1 = 201 + isrow
+            ij0 = 241 - isrow   ;   ij1 = 241 - isrow
             !                                        ! =======================
          CASE DEFAULT                                !    ORCA R05 or R025
@@ -223,11 +219,11 @@
          CASE ( 1 )                                  !  ORCA_R1 configurations
             !                                        ! =======================
-            !  This dirty section will be suppressed by simplification process:
-            !  all this will come back in input files
-            IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  
-            ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 39
-            ENDIF                                      
+            ! This dirty section will be suppressed by simplification process:
+            ! all this will come back in input files
+            ! Currently these hard-wired indices relate to configuration with
+            ! extend grid (jpjglo=332)
+            isrow = 332 - jpjglo
             ii0 = 282           ;   ii1 = 282
-            ij0 = 201 + isrow   ;   ij1 = 201 + isrow
+            ij0 = 240 - isrow   ;   ij1 = 240 - isrow
             !                                        ! =======================
          CASE DEFAULT                                !    ORCA R05 or R025
@@ -275,11 +271,11 @@
          CASE ( 1 )                                  !  ORCA_R1 configurations
             !                                        ! =======================
-            !  This dirty section will be suppressed by simplification process:
-            !  all this will come back in input files
-            IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  
-            ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 39
-            ENDIF                                      
+            ! This dirty section will be suppressed by simplification process:
+            ! all this will come back in input files
+            ! Currently these hard-wired indices relate to configuration with
+            ! extend grid (jpjglo=332)
+            isrow = 332 - jpjglo
             ii0 = 331           ;   ii1 = 331
-            ij0 = 176 + isrow   ;   ij1 = 176 + isrow
+            ij0 = 215 - isrow   ;   ij1 = 215 - isrow
             !                                        ! =======================
          CASE DEFAULT                                !    ORCA R05 or R025
@@ -327,11 +323,11 @@
          CASE ( 1 )                                  !  ORCA_R1 configurations
             !                                        ! =======================
-            !  This dirty section will be suppressed by simplification process:
-            !  all this will come back in input files
-            IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  
-            ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 39
-            ENDIF                                      
+            ! This dirty section will be suppressed by simplification process:
+            ! all this will come back in input files
+            ! Currently these hard-wired indices relate to configuration with
+            ! extend grid (jpjglo=332)
+            isrow = 332 - jpjglo
             ii0 = 297           ;   ii1 = 297
-            ij0 = 230 + isrow   ;   ij1 = 230 + isrow
+            ij0 = 269 - isrow   ;   ij1 = 269 - isrow
             !                                        ! =======================
          CASE DEFAULT                                !    ORCA R05 or R025
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 5539)
@@ -46,4 +46,6 @@
    USE iom
    USE ioipsl
+   USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities     
+
 #if defined key_lim2
    USE limwri_2 
@@ -125,4 +127,5 @@
       !!
       INTEGER                      ::   ji, jj, jk              ! dummy loop indices
+      INTEGER                      ::   jkbot                   !
       REAL(wp)                     ::   zztmp, zztmpx, zztmpy   ! 
       !!
@@ -148,4 +151,7 @@
          CALL iom_put( "e3w" , fse3w_n(:,:,:) )
       ENDIF
+
+      CALL iom_put( "ssh" , sshn )                 ! sea surface height
+      if( iom_use('ssh2') )   CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
       
       CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
@@ -154,5 +160,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               z2d(ji,jj) = tsn(ji,jj,MAX(mbathy(ji,jj),1),jp_tem)
+               jkbot = mbkt(ji,jj)
+               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem)
             END DO
          END DO
@@ -165,8 +172,25 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               z2d(ji,jj) = tsn(ji,jj,MAX(mbathy(ji,jj),1),jp_sal)
+               jkbot = mbkt(ji,jj)
+               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_sal)
             END DO
          END DO
          CALL iom_put( "sbs", z2d )                ! bottom salinity
+      ENDIF
+
+      IF ( iom_use("taubot") ) THEN                ! bottom stress
+         z2d(:,:) = 0._wp
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  &
+                      &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  )      
+               zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  &
+                      &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  ) 
+               z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 
+               !
+            ENDDO
+         ENDDO
+         CALL lbc_lnk( z2d, 'T', 1. )
+         CALL iom_put( "taubot", z2d )           
       ENDIF
          
@@ -176,9 +200,15 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               z2d(ji,jj) = un(ji,jj,MAX(mbathy(ji,jj),1))
+               jkbot = mbku(ji,jj)
+               z2d(ji,jj) = un(ji,jj,jkbot)
             END DO
          END DO
          CALL iom_put( "sbu", z2d )                ! bottom i-current
       ENDIF
+#if defined key_dynspg_ts
+      CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current
+#else
+      CALL iom_put(  "ubar", un_b(:,:)        )    ! barotropic i-current
+#endif
       
       CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current
@@ -187,8 +217,25 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               z2d(ji,jj) = vn(ji,jj,MAX(mbathy(ji,jj),1))
+               jkbot = mbkv(ji,jj)
+               z2d(ji,jj) = vn(ji,jj,jkbot)
             END DO
          END DO
          CALL iom_put( "sbv", z2d )                ! bottom j-current
+      ENDIF
+#if defined key_dynspg_ts
+      CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic j-current
+#else
+      CALL iom_put(  "vbar", vn_b(:,:)        )    ! barotropic j-current
+#endif
+
+      CALL iom_put( "woce", wn )                   ! vertical velocity
+      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
+         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
+         z2d(:,:) = rau0 * e12t(:,:)
+         DO jk = 1, jpk
+            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
+         END DO
+         CALL iom_put( "w_masstr" , z3d )  
+         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
       ENDIF
 
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 5539)
@@ -72,4 +72,5 @@
       !!----------------------------------------------------------------------
       INTEGER ::   jc            ! dummy loop indices
+      INTEGER :: isrow           ! local index
       !!----------------------------------------------------------------------
       
@@ -91,7 +92,13 @@
          CASE ( 1 )                                  ! ORCA_R1 configuration
             !                                        ! =======================
+            ! This dirty section will be suppressed by simplification process:
+            ! all this will come back in input files
+            ! Currently these hard-wired indices relate to configuration with
+            ! extend grid (jpjglo=332)
+            isrow = 332 - jpjglo
+            !
             ncsnr(1)   = 1    ; ncstt(1)   = 0           ! Caspian Sea
-            ncsi1(1)   = 332  ; ncsj1(1)   = 203
-            ncsi2(1)   = 344  ; ncsj2(1)   = 235
+            ncsi1(1)   = 332  ; ncsj1(1)   = 243 - isrow
+            ncsi2(1)   = 344  ; ncsj2(1)   = 275 - isrow
             ncsir(1,1) = 1    ; ncsjr(1,1) = 1
             !                                        
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 5539)
@@ -162,51 +162,47 @@
             !                                             ! =====================
             ! This dirty section will be suppressed by simplification process: all this will come back in input files
-            ! Currently these hard-wired indices relate to the original (pre-v3.6) configuration 
+            ! Currently these hard-wired indices relate to configuration with
+            ! extend grid (jpjglo=332)
             ! which had a grid-size of 362x292.
-            ! This grid has been extended southwards for use with the under ice-shelf options (isf) introduced in v3.6. 
-            ! The original domain can still be used optionally if the isf code is not activated. 
-            ! An adjustment (isrow) is made to the hard-wired indices if the extended domain (362x332) is being used.
             ! 
-            IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  ! Using pre-v3.6 files or adjusted start row from isf-extended grid
-            ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 40 ! Using full isf­extended domain. 
-            ENDIF     
+            isrow = 332 - jpjglo
             !
             ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u = 20 km)
-            ij0 = 201 + isrow   ;   ij1 = 201 + isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
+            ij0 = 201 + isrow   ;   ij1 = 241 - isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: Gibraltar : e2u reduced to 20 km'
 
             ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u = 10 km)
-            ij0 = 208 + isrow   ;   ij1 = 208 + isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
+            ij0 = 208 + isrow   ;   ij1 = 248 - isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: Bhosporus : e2u reduced to 10 km'
 
             ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v = 13 km)
-            ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  13.e3
+            ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  13.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: Lombok : e1v reduced to 10 km'
 
             ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on]
-            ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  8.e3
+            ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  8.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: Sumba : e1v reduced to 8 km'
 
             ii0 =  53           ;   ii1 =  53        ! Ombai Strait (e1v = 13 km)
-            ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3
+            ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: Ombai : e1v reduced to 13 km'
 
             ii0 =  56           ;   ii1 =  56        ! Timor Passage (e1v = 20 km)
-            ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3
+            ij0 = 124 + isrow   ;   ij1 = 145 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: Timor Passage : e1v reduced to 20 km'
 
             ii0 =  55           ;   ii1 =  55        ! West Halmahera Strait (e1v = 30 km)
-            ij0 = 141 + isrow   ;   ij1 = 142 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3
+            ij0 = 141 + isrow   ;   ij1 = 182 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: W Halmahera : e1v reduced to 30 km'
 
             ii0 =  58           ;   ii1 =  58        ! East Halmahera Strait (e1v = 50 km)
-            ij0 = 141 + isrow   ;   ij1 = 142 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3
+            ij0 = 141 + isrow   ;   ij1 = 182 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: E Halmahera : e1v reduced to 50 km'
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 5539)
@@ -402,48 +402,44 @@
       IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
          !                                                 ! Increased lateral friction near of some straits
-         ! This dirty section will be suppressed by simplification process: all this will come back in input files
-         ! Currently these hard-wired indices relate to the original (pre-v3.6) configuration 
-         ! which had a grid-size of 362x292.
-         ! This grid has been extended southwards for use with the under ice-shelf options (isf) introduced in v3.6. 
-         ! The original domain can still be used optionally if the isf code is not activated. 
-         ! An adjustment (isrow) is made to the hard-wired indices if the extended domain (362x332) is being used.
-         ! 
-         IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  ! Using pre-v3.6 files or adjusted start row from isf-extended grid
-         ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 40 ! Using full isf­extended domain. 
-         ENDIF     
-
+         ! This dirty section will be suppressed by simplification process:
+         ! all this will come back in input files
+         ! Currently these hard-wired indices relate to configuration with
+         ! extend grid (jpjglo=332)
+         !
+         isrow = 332 - jpjglo
+         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
          IF(lwp) WRITE(numout,*) '      Gibraltar '
          ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait 
-         ij0 = 201 + isrow   ;   ij1 = 201 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+         ij0 = 201 + isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
 
          IF(lwp) WRITE(numout,*) '      Bhosporus '
          ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait 
-         ij0 = 208 + isrow   ;   ij1 = 208 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+         ij0 = 208 + isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
 
          IF(lwp) WRITE(numout,*) '      Makassar (Top) '
          ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top) 
-         ij0 = 149 + isrow   ;   ij1 = 150 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
+         ij0 = 149 + isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
 
          IF(lwp) WRITE(numout,*) '      Lombok '
          ii0 =  44           ;   ii1 =  44        ! Lombok Strait 
-         ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+         ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
 
          IF(lwp) WRITE(numout,*) '      Ombai '
          ii0 =  53           ;   ii1 =  53        ! Ombai Strait 
-         ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+         ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
 
          IF(lwp) WRITE(numout,*) '      Timor Passage '
          ii0 =  56           ;   ii1 =  56        ! Timor Passage 
-         ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+         ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
 
          IF(lwp) WRITE(numout,*) '      West Halmahera '
          ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait 
-         ij0 = 141 + isrow   ;   ij1 = 142 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
+         ij0 = 141 + isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
 
          IF(lwp) WRITE(numout,*) '      East Halmahera '
          ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait 
-         ij0 = 141 + isrow   ;   ij1 = 142 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
+         ij0 = 141 + isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
          !
       ENDIF
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 5539)
@@ -1123,17 +1123,13 @@
       IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
          !                                             ! =====================
-         ! This dirty section will be suppressed by simplification process: all this will come back in input files
-         ! Currently these hard-wired indices relate to the original (pre-v3.6) configuration 
+         ! This dirty section will be suppressed by simplification process:
+         ! all this will come back in input files
+         ! Currently these hard-wired indices relate to configuration with
+         ! extend grid (jpjglo=332)
          ! which had a grid-size of 362x292.
-         ! This grid has been extended southwards for use with the under ice-shelf options (isf) introduced in v3.6. 
-         ! The original domain can still be used optionally if the isf code is not activated. 
-         ! An adjustment (isrow) is made to the hard-wired indices if the extended domain (362x332) is being used.
-         ! 
-         IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  ! Using pre-v3.6 files or adjusted start row from isf-extended grid
-         ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 40 ! Using full isf­extended domain. 
-         ENDIF     
+         isrow = 332 - jpjglo
          !
          ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u was modified)
-         ij0 = 201 + isrow   ;   ij1 = 201 + isrow
+         ij0 = 241 - isrow   ;   ij1 = 241 - isrow
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
@@ -1156,5 +1152,5 @@
          !
          ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
-         ij0 = 208 + isrow   ;   ij1 = 208 + isrow
+         ij0 = 248 - isrow   ;   ij1 = 248 - isrow
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
@@ -1177,5 +1173,5 @@
          !
          ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v was modified)
-         ij0 = 124 + isrow   ;   ij1 = 125 + isrow
+         ij0 = 164 - isrow   ;   ij1 = 165 - isrow
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
@@ -1193,5 +1189,5 @@
          !
          ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
-         ij0 = 124 + isrow   ;   ij1 = 125 + isrow
+         ij0 = 164 - isrow   ;   ij1 = 165 - isrow
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
@@ -1209,5 +1205,5 @@
          !
          ii0 =  53          ;   ii1 =  53        ! Ombai Strait (e1v was modified)
-         ij0 = 124 + isrow  ;   ij1 = 125  + isrow  
+         ij0 = 164 - isrow  ;   ij1 = 165  - isrow  
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
@@ -1224,6 +1220,6 @@
          END DO
          !
-         ii0 =  56             ;   ii1 =  56        ! Timor Passage (e1v was modified)
-         ij0 = 124  + isrow    ;   ij1 = 125  + isrow  
+         ii0 =  56            ;   ii1 =  56        ! Timor Passage (e1v was modified)
+         ij0 = 164 - isrow    ;   ij1 = 165  - isrow  
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
@@ -1240,6 +1236,6 @@
          END DO
          !
-         ii0 =  55             ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
-         ij0 = 141  + isrow    ;   ij1 = 142  + isrow  
+         ii0 =  55            ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
+         ij0 = 181 - isrow    ;   ij1 = 182 - isrow  
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
@@ -1256,6 +1252,6 @@
          END DO
          !
-         ii0 =  58             ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
-         ij0 = 141  + isrow    ;   ij1 = 142  + isrow  
+         ii0 =  58            ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
+         ij0 = 181 - isrow    ;   ij1 = 182 - isrow  
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5539)
@@ -484,19 +484,4 @@
             CALL iom_close( inum )
             mbathy(:,:) = INT( bathy(:,:) )
-            !
-            ! CL :  add Amazon deeper
-            IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration 
-                 ii0 = 230   ;   ii1 = 245                   ! Amazon area 
-                 ij0 = 140   ;   ij1 = 155                   ! no ocean shallower than 30 meters
-                 DO ji = mi0(ii0), mi1(ii1)
-                    DO jj = mj0(ij0), mj1(ij1)
-                       IF( bathy(ji,jj) .LE. 30. .AND. bathy(ji,jj) .GT. 0.0 ) bathy(ji,jj) = 30._wp
-                    END DO
-                 END DO
-                 IF(lwp) WRITE(numout,*)
-                 IF(lwp) WRITE(numout,*) '             orca_r1: Amazon area not shallower than 30 meters for: '
-                 IF(lwp) WRITE(numout,*) '                    Longitude index ',ii0, ii0
-                 IF(lwp) WRITE(numout,*) '                    Latitude index  ',ij0, ij0
-            ENDIF
             !                                                ! =====================
             IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 5539)
@@ -237,5 +237,5 @@
       !
       CALL wrk_dealloc( jpi  , jpj+2, zwu               )
-      CALL wrk_dealloc( jpi+4, jpj  , zwv, kjstart = -1 )
+      CALL wrk_dealloc( jpi+4, jpj  , zwv, kistart = -1 )
       !
       IF( nn_timing == 1 )  CALL timing_stop('div_cur')
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 5539)
@@ -266,6 +266,6 @@
                ! Add volume filter correction: compatibility with tracer advection scheme
                ! => time filter + conservation correction (only at the first level)
-               fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
-            !
+               fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) &
+                              &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
             ENDIF
             !
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 5539)
@@ -21,5 +21,4 @@
    USE domvvl          ! Variable volume
    USE divcur          ! hor. divergence and curl      (div & cur routines)
-   USE iom             ! I/O library
    USE restart         ! only for lrst_oce
    USE in_out_manager  ! I/O manager
@@ -31,5 +30,4 @@
    USE bdy_par         
    USE bdydyn2d        ! bdy_ssh routine
-   USE iom
 #if defined key_agrif
    USE agrif_opa_update
@@ -137,6 +135,4 @@
       !                                           !           outputs            !
       !                                           !------------------------------!
-      CALL iom_put( "ssh" , sshn )   ! sea surface height
-      if( iom_use('ssh2') )   CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
       !
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
@@ -228,22 +224,4 @@
 #endif
       !
-      !                                           !------------------------------!
-      !                                           !           outputs            !
-      !                                           !------------------------------!
-      CALL iom_put( "woce", wn )   ! vertical velocity
-      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
-         CALL wrk_alloc( jpi, jpj, z2d ) 
-         CALL wrk_alloc( jpi, jpj, jpk, z3d ) 
-         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
-         z2d(:,:) = rau0 * e12t(:,:)
-         DO jk = 1, jpk
-            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
-         END DO
-         CALL iom_put( "w_masstr" , z3d )  
-         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
-         CALL wrk_dealloc( jpi, jpj, z2d  ) 
-         CALL wrk_dealloc( jpi, jpj, jpk, z3d ) 
-      ENDIF
-      !
       IF( nn_timing == 1 )  CALL timing_stop('wzv')
 
@@ -290,5 +268,5 @@
       ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
          sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
-         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * ssmask(:,:)
+         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) - rnf_b(:,:) + rnf(:,:) ) * ssmask(:,:)
          sshn(:,:) = ssha(:,:)                           ! now <-- after
       ENDIF
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 5539)
@@ -91,5 +91,6 @@
    INTEGER ::   nitrst                !: time step at which restart file should be written
    LOGICAL ::   lrst_oce              !: logical to control the oce restart write 
-   INTEGER ::   numror, numrow        !: logical unit for cean restart (read and write)
+   INTEGER ::   numror = 0            !: logical unit for ocean restart (read). Init to 0 is needed for SAS (in daymod.F90)
+   INTEGER ::   numrow                !: logical unit for ocean restart (write)
    INTEGER ::   nrst_lst              !: number of restart to output next
 
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 5539)
@@ -663,4 +663,5 @@
 
       ! --- evaporation minus precipitation --- !
+      zsnw(:,:) = 0._wp
       CALL lim_thd_snwblow( pfrld, zsnw )          ! snow redistribution by wind
       emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * ( 1._wp - zsnw )
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 5539)
@@ -613,4 +613,5 @@
 
       ! --- evaporation minus precipitation --- !
+      zsnw(:,:) = 0._wp
       CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing 
       emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 5539)
@@ -780,5 +780,5 @@
       xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
       !
-      ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'S_QsrOce' ) + cpl_freq( 'S_QsrMix' )
+      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' )
@@ -1379,5 +1379,5 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ! for LIM3
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zevap, zsnw, zqns_oce, zqsr_oce, zqprec_ice, zqemp_oce ! for LIM3
       !!----------------------------------------------------------------------
       !
@@ -1524,4 +1524,5 @@
 
       ! --- heat flux associated with emp --- !
+      zsnw(:,:) = 0._wp
       CALL lim_thd_snwblow( p_frld, zsnw )  ! snow distribution over ice after wind blowing
       zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &      ! evap
@@ -1556,5 +1557,4 @@
 
       CALL wrk_dealloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
-
 #else
 
@@ -1624,4 +1624,17 @@
          ENDDO
       ENDIF
+
+#if defined key_lim3
+      CALL wrk_alloc( jpi,jpj, zqsr_oce ) 
+      ! --- solar flux over ocean --- !
+      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
+      zqsr_oce = 0._wp
+      WHERE( p_frld /= 0._wp )  zqsr_oce(:,:) = ( zqsr_tot(:,:) - SUM( a_i * zqsr_ice, dim=3 ) ) / p_frld(:,:)
+
+      IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
+      ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF
+
+      CALL wrk_dealloc( jpi,jpj, zqsr_oce ) 
+#endif
 
       IF( ln_mixcpl ) THEN
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 5539)
@@ -509,5 +509,5 @@
       CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
 
-      CALL wrk_dealloc( jpi,jpj, ztmp )
+      CALL wrk_dealloc( jpi,jpj, ztmp, zpice )
       CALL wrk_dealloc( jpi,jpj,ncat, ztmpn )
       !
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5539)
@@ -216,8 +216,8 @@
             zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
                                  CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
+            ! clem: evap_ice is forced to 0 in coupled mode for now 
+            !       but it needs to be changed (along with modif in limthd_dh) once heat flux from evap will be avail. from atm. models
+            evap_ice  (:,:,:) = 0._wp   ;   devap_ice (:,:,:) = 0._wp
             IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
-            ! Latent heat flux is forced to 0 in coupled: it is included in qns (non-solar heat flux)
-            evap_ice  (:,:,:) = 0._wp
-            devap_ice (:,:,:) = 0._wp
          END SELECT
          CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 5539)
@@ -53,4 +53,5 @@
    USE timing           ! Timing
    USE sbcwave          ! Wave module
+   USE bdy_par          ! Require lk_bdy
 
    IMPLICIT NONE
@@ -343,5 +344,7 @@
       !                                            ! ---------------------------------------- !
       !
-      IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc
+      IF ( .NOT. lk_bdy ) then
+         IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc
+      ENDIF
                                                          ! (caution called before sbc_ssm)
       !
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 5539)
@@ -209,5 +209,5 @@
       zfact = 0.5_wp
       !
-      IF( ln_rnf_depth ) THEN      !==   runoff distributed over several levels   ==!
+      IF( ln_rnf_depth .OR. ln_rnf_depth_ini ) THEN      !==   runoff distributed over several levels   ==!
          IF( lk_vvl ) THEN             ! variable volume case
             DO jj = 1, jpj                   ! update the depth over which runoffs are distributed
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/STO/stopar.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/STO/stopar.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/STO/stopar.F90	(revision 5539)
@@ -112,7 +112,4 @@
    INTEGER         :: nn_trc_ord = 1          ! order of autoregressive processes
 
-   ! Public array with density correction
-   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: drho_ran
-
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -263,5 +260,4 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
 
-      IF ( lk_agrif ) CALL ctl_stop('EOS stochastic parametrization is not compatible with AGRIF')
       ! Read namsto namelist : stochastic parameterization
       REWIND( numnam_ref )              ! Namelist namdyn_adv in reference namelist : Momentum advection scheme
@@ -341,4 +337,5 @@
       ENDIF
       IF( ln_sto_eos ) THEN
+         IF ( lk_agrif ) CALL ctl_stop('EOS stochastic parametrization is not compatible with AGRIF')
          IF(lwp) WRITE(numout,*) '       - stochastic equation of state'
          ALLOCATE(jsto_eosi(nn_sto_eos))
@@ -668,7 +665,4 @@
       ! ------------------------------------------
       IF( ln_rststo ) CALL sto_rst_read
-
-      ! Allocate drho_ran
-      ALLOCATE(drho_ran(jpi,jpj,jpk))
 
    END SUBROUTINE sto_par_init
@@ -781,6 +775,4 @@
                CALL iom_rstput( kt, nitrst, numstow, clsto3d , sto3d(:,:,:,jsto) )
             END DO
-            ! Save drho_ran in restart file
-            CALL iom_rstput( kt, nitrst, numstow, 'drho' , drho_ran(:,:,:) )
             ! close the restart file
             CALL iom_close( numstow )
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 5539)
@@ -27,4 +27,5 @@
    USE dom_oce         ! ocean space and time domain variables 
    USE sbc_oce         ! surface boundary condition: ocean
+   USE sbcrnf          ! river runoffs
    USE zdf_oce         ! ocean vertical mixing
    USE domvvl          ! variable volume
@@ -278,5 +279,5 @@
 
       !!     
-      LOGICAL  ::   ll_tra_hpg, ll_traqsr   ! local logical
+      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf   ! local logical
       INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
       REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
@@ -293,7 +294,9 @@
          ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
          ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
+         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
       ELSE                          
          ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
          ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
+         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs
       ENDIF
       !
@@ -319,18 +322,23 @@
                   !
                   IF( jk == 1 ) THEN           ! first level 
-                     ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) )
+                     ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) + rnf(ji,jj) - rnf_b(ji,jj) )
                      ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
                   ENDIF
+
                   IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )   &     ! solar penetration (temperature only)
                      &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
 
-                   ze3t_f = 1.e0 / ze3t_f
-                   ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
-                   ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
-                   !
-                   IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
-                      ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
-                      pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
-                   ENDIF
+                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )   & 		      ! river runoffs
+                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
+                     &                              * fse3t_n(ji,jj,jk) / h_rnf(ji,jj)
+
+                  ze3t_f = 1.e0 / ze3t_f
+                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
+                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
+                  !
+                  IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
+                     ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
+                     pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
+                  ENDIF
                END DO
             END DO
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5539)
@@ -101,9 +101,10 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( lk_tide    )   CALL sbc_tide( kstp )
-      IF( lk_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
-
+      IF( lk_bdy     )  THEN
+         IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib 
+                         CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
+      ENDIF
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
                                                       ! clem: moved here for bdy ice purpose
-
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Update stochastic parameters and random T/S fluctuations
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5539)
@@ -27,4 +27,5 @@
    USE sbc_oce          ! surface boundary condition: ocean
    USE sbctide          ! Tide initialisation
+   USE sbcapr           ! surface boundary condition: ssh_ib required by bdydta 
 
    USE traqsr           ! solar radiation penetration      (tra_qsr routine)
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/wrk_nemo.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/wrk_nemo.F90	(revision 5538)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/wrk_nemo.F90	(revision 5539)
@@ -121,4 +121,5 @@
 
    LOGICAL ::   linit = .FALSE.
+   LOGICAL ::   ldebug = .FALSE.
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011)
@@ -486,4 +487,5 @@
       
       IF( SUM( tree(ii)%ishape ) == 0 ) THEN                    ! create a new branch 
+         IF(ldebug) PRINT *, 'create new branch ', ii,ishape, isrt, itype
          tree(ii)%itype = itype                                        ! define the type of this branch 
          tree(ii)%ishape(:) = ishape(:)                                ! define the shape of this branch 
@@ -515,4 +517,5 @@
          tree(ii)%current%next%in_use = .FALSE.                        ! this leaf is not yet used
          tree(ii)%current%next%indic = tree(ii)%current%indic + 1      ! number of this leaf
+         IF(ldebug) PRINT *, 'add a leaf ', ii, tree(ii)%current%indic
          tree(ii)%current%next%prev => tree(ii)%current                ! previous leaf of the new leaf is the current leaf
          tree(ii)%current%next%next => NULL()                          ! next leaf is not yet defined
