Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 4934)
@@ -728,4 +728,5 @@
    ln_traadv_qck    =  .false.  !  QUICKEST scheme
    ln_traadv_msc_ups=  .false.  !  use upstream scheme within muscl
+   ln_traadv_tvd_zts=  .false.  !  TVD scheme with sub-timestepping of vertical tracer advection
 /
 !-----------------------------------------------------------------------
@@ -802,4 +803,5 @@
    ln_dynadv_cen2= .false. !  flux form - 2nd order centered scheme
    ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme
+   ln_dynzad_zts = .false. !  Use (T) sub timestepping for vertical momentum advection
 /
 !-----------------------------------------------------------------------
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/CONFIG/makenemo
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/CONFIG/makenemo	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/CONFIG/makenemo	(revision 4934)
@@ -38,7 +38,8 @@
 #
 # - NEW_CONF    : configuration to be created
-# - REF_CONF    : reference configuration to build the new one
+# - REF_CONF    : reference configuration to build the new one from
 # - CMP_NAM     : compiler name
 # - NBR_PRC     : number of processes used to compile
+# - USP_CONF    : unsupported (external) configuration to build the new one from
 # - NEM_SUBDIR  : NEMO subdirectory used (specified)
 #
@@ -51,4 +52,9 @@
 # - TOOLS_DIR   :   "    "    "
 # - NEMO_DIR    :   "    "    "
+# - REMOTE_CTL  : URL link to a remote resource list for an external configuration 
+#                 which is not part of the reference suite
+# - LOCAL_REF   : Nearest reference configuration to an external configuration 
+#                 which is not part of the reference suite
+#                 (used to populate work directories if remote access is not available)
 #
 # EXAMPLES
@@ -83,4 +89,5 @@
 x_n="";
 x_r="";
+x_u="";
 x_m="";
 x_t="";
@@ -106,4 +113,6 @@
 export AGRIFUSE=10
 declare -a TAB
+declare -a REMOTE_CTL
+declare -a LOCAL_REF
 list_key=0
 chk_key=1
@@ -114,10 +123,10 @@
 #-
 #- Choice of the options ---
-while getopts :hd:n:r:m:j:e:s:v:t:k: V
+while getopts :hd:n:r:u:m:j:e:s:v:t:k: V
 do
     case $V in
 	(h) x_h=${OPTARG};
         echo "Usage	 : "${b_n} \
-	    " [-h] [-n name] [-m arch] [-d "dir1 dir2"] [-r conf] [-s Path] [-e Path] [-j No] [-v No] [-k 0/1]";
+	    " [-h] [-n name] [-m arch] [-d "dir1 dir2"] [-r conf] [-u conf] [-s Path] [-e Path] [-j No] [-v No] [-k 0/1]";
 	echo " -h	        : help";
 	echo " -h institute : specific help for consortium members";
@@ -126,4 +135,5 @@
 	echo " -d dir       : choose NEMO sub-directories";
 	echo " -r conf      : choose reference configuration";
+	echo " -u conf      : choose an unsupported (external) configuration";
         echo " -s Path      : choose alternative location for NEMO main directory";
         echo " -e Path      : choose alternative location for MY_SRC directory";
@@ -139,4 +149,6 @@
 	echo "Available configurations :"; cat ${CONFIG_DIR}/cfg.txt;
 	echo "";
+        echo "Available unsupported (external) configurations :"; cat ${CONFIG_DIR}/uspcfg.txt;
+	echo "";
 	echo "Example to remove bad configuration ";
 	echo "./makenemo -n MY_CONFIG clean_config";
@@ -161,4 +173,5 @@
 	(n)  x_n=${OPTARG};;
 	(r)  x_r=${OPTARG};;
+	(u)  x_u=${OPTARG};;
 	(m)  x_m=${OPTARG};;
 	(j)  x_j=${OPTARG};;
@@ -220,4 +233,5 @@
 NEM_SUBDIR=${x_d}
 REF_CONF=${x_r}
+USP_CONF=${x_u}
 NEMO_TDIR=${x_t:-$NEMO_TDIR}
 export NEMO_TDIR=${NEMO_TDIR:-$CONFIG_DIR}
@@ -228,4 +242,7 @@
     echo "Available configurations :"
     cat ${CONFIG_DIR}/cfg.txt
+    echo
+    echo "Available unsupported (external) configurations :" 
+    cat ${CONFIG_DIR}/uspcfg.txt
     exit
 fi
@@ -238,6 +255,6 @@
 
 if [ ${#NEW_CONF} -eq 0 ] ; then
-    if [ ${#NEM_SUBDIR} -eq 0 -a ${#REF_CONF} -eq 0 ]; then
-	echo "You are  installing a new configuration"
+    if [ ${#NEM_SUBDIR} -eq 0 ] && [ ${#REF_CONF} -eq 0 ] && [ ${#USP_CONF} -eq 0 ] ; then
+	echo "You are installing a new default (ORCA2_LIM) configuration"
 	ind=0
 	. ${COMPIL_DIR}/Fread_dir.sh OPA_SRC	 YES
@@ -248,6 +265,6 @@
 	. ${COMPIL_DIR}/Fread_dir.sh OFF_SRC	 NO
 	REF_CONF=ORCA2_LIM
-    elif [ ${#NEM_SUBDIR} -gt 0 ] && [ ${#REF_CONF} -eq 0 ]; then
-	echo "You are  installing a new configuration"
+    elif [ ${#NEM_SUBDIR} -gt 0 ] && [ ${#REF_CONF} -eq 0 ] && [ ${#USP_CONF} -eq 0 ] ; then
+	echo "You are  installing a new configuration based on ORCA2_LIM"
 	TAB=( ${NEM_SUBDIR} )
 	REF_CONF=ORCA2_LIM
@@ -255,7 +272,19 @@
 	echo "You are  installing a new configuration based on ${REF_CONF}"
 	. ${COMPIL_DIR}/Fcopy_dir.sh ${REF_CONF}
+    elif [ ${#NEM_SUBDIR} -eq 0 ] && [ ${#USP_CONF} -gt 0 ]; then
+	echo "You are  installing a new configuration based on the unsupported (external) ${USP_CONF}"
+	. ${COMPIL_DIR}/Fcopy_extdir.sh ${USP_CONF}  
+        #echo "TTT " ${TAB}
+        #echo "RRR " ${REMOTE_CTL}
+        #echo "LLL " ${LOCAL_REF}
     fi
     NEW_CONF=${x_n}
-    . ${COMPIL_DIR}/Fmake_config.sh ${NEW_CONF} ${REF_CONF}
+
+    if [ ${#USP_CONF} -gt 0 ]; then
+      . ${COMPIL_DIR}/Fmake_extconfig.sh ${NEW_CONF} ${LOCAL_REF}
+      . ${COMPIL_DIR}/Ffetch_extdir.sh ${NEW_CONF} ${REMOTE_CTL}  
+    else
+      . ${COMPIL_DIR}/Fmake_config.sh ${NEW_CONF} ${REF_CONF}
+    fi
 else
     sed -e "/${NEW_CONF} /d"  ${CONFIG_DIR}/cfg.txt >  ${COMPIL_DIR}/cfg.tmp
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/CONFIG/uspcfg.txt
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/CONFIG/uspcfg.txt	(revision 4934)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/CONFIG/uspcfg.txt	(revision 4934)
@@ -0,0 +1,1 @@
+ORCA1_CICE # ORCA2_LIM # OPA_SRC TOP_SRC  # http://gws-access.ceda.ac.uk/public/nemo/uspconfigs/ORCA1_CICE/v3.6.0/ORCA1_CICE_ctl.txt
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 4934)
@@ -30,4 +30,5 @@
    LOGICAL, PUBLIC ::   ln_dynadv_cen2  !: flux form - 2nd order centered scheme flag
    LOGICAL, PUBLIC ::   ln_dynadv_ubs   !: flux form - 3rd order UBS scheme flag
+   LOGICAL, PUBLIC ::   ln_dynzad_zts   !: vertical advection with sub-timestepping (requires vector form)
    
    INTEGER ::   nadv   ! choice of the formulation and scheme for the advection
@@ -64,7 +65,10 @@
                       CALL dyn_keg     ( kt )    ! vector form : horizontal gradient of kinetic energy
                       CALL dyn_zad     ( kt )    ! vector form : vertical advection
-      CASE ( 1 ) 
+      CASE ( 1 )     
+                      CALL dyn_keg     ( kt )    ! vector form : horizontal gradient of kinetic energy
+                      CALL dyn_zad_zts ( kt )    ! vector form : vertical advection with sub-timestepping
+      CASE ( 2 ) 
                       CALL dyn_adv_cen2( kt )    ! 2nd order centered scheme
-      CASE ( 2 )   
+      CASE ( 3 )   
                       CALL dyn_adv_ubs ( kt )    ! 3rd order UBS      scheme
       !
@@ -91,5 +95,5 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs
+      NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts
       !!----------------------------------------------------------------------
 
@@ -111,4 +115,5 @@
          WRITE(numout,*) '          2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2
          WRITE(numout,*) '          3rd order UBS advection scheme     ln_dynadv_ubs  = ', ln_dynadv_ubs
+         WRITE(numout,*) '      Sub timestepping of vertical advection ln_dynzad_zts  = ', ln_dynzad_zts
       ENDIF
 
@@ -120,9 +125,12 @@
 
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' )
+      IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec )   &
+          CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' )
 
       !                               ! Set nadv
       IF( ln_dynadv_vec  )   nadv =  0 
-      IF( ln_dynadv_cen2 )   nadv =  1
-      IF( ln_dynadv_ubs  )   nadv =  2
+      IF( ln_dynzad_zts  )   nadv =  1
+      IF( ln_dynadv_cen2 )   nadv =  2
+      IF( ln_dynadv_ubs  )   nadv =  3
       IF( lk_esopa       )   nadv = -1
 
@@ -130,6 +138,7 @@
          WRITE(numout,*)
          IF( nadv ==  0 )   WRITE(numout,*) '         vector form : keg + zad + vor is used'
-         IF( nadv ==  1 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used'
-         IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : UBS       scheme is used'
+         IF( nadv ==  1 )   WRITE(numout,*) '         vector form : keg + zad_zts + vor is used'
+         IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used'
+         IF( nadv ==  3 )   WRITE(numout,*) '         flux form   : UBS       scheme is used'
          IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection formulation'
       ENDIF
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 4934)
@@ -27,5 +27,6 @@
    PRIVATE
    
-   PUBLIC   dyn_zad   ! routine called by step.F90
+   PUBLIC   dyn_zad       ! routine called by dynadv.F90
+   PUBLIC   dyn_zad_zts   ! routine called by dynadv.F90
 
    !! * Substitutions
@@ -83,5 +84,5 @@
          DO jj = 2, jpj                   ! vertical fluxes 
             DO ji = fs_2, jpi             ! vector opt.
-               zww(ji,jj) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
+               zww(ji,jj) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
             END DO
          END DO
@@ -95,8 +96,8 @@
       DO jj = 2, jpjm1              ! Surface and bottom values set to zero
          DO ji = fs_2, fs_jpim1           ! vector opt.
-            zwuw(ji,jj, 1:miku(ji,jj) ) = 0.e0
-            zwvw(ji,jj, 1:mikv(ji,jj) ) = 0.e0
-            zwuw(ji,jj,jpk) = 0.e0
-            zwvw(ji,jj,jpk) = 0.e0
+            zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp
+            zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp
+            zwuw(ji,jj,jpk) = 0._wp
+            zwvw(ji,jj,jpk) = 0._wp
          END DO  
       END DO
@@ -132,4 +133,146 @@
    END SUBROUTINE dyn_zad
 
+   SUBROUTINE dyn_zad_zts ( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dynzad_zts  ***
+      !! 
+      !! ** Purpose :   Compute the now vertical momentum advection trend and 
+      !!      add it to the general trend of momentum equation. This version
+      !!      uses sub-timesteps for improved numerical stability with small
+      !!      vertical grid sizes. This is especially relevant when using 
+      !!      embedded ice with thin surface boxes.
+      !!
+      !! ** Method  :   The now vertical advection of momentum is given by:
+      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
+      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
+      !!      Add this trend to the general trend (ua,va):
+      !!         (ua,va) = (ua,va) + w dz(u,v)
+      !!
+      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
+      !!              - Save the trends in (ztrdu,ztrdv) ('key_trddyn')
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
+      !
+      INTEGER  ::   ji, jj, jk, jl  ! dummy loop indices
+      INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection
+      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
+      REAL(wp) ::   zua, zva        ! temporary scalars
+      REAL(wp) ::   zr_rdt          ! temporary scalar
+      REAL(wp) ::   z2dtzts         ! length of Euler forward sub-timestep for vertical advection
+      REAL(wp) ::   zts             ! length of sub-timestep for vertical advection
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zwuw , zwvw, zww
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztrdu, ztrdv
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zus , zvs
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('dyn_zad_zts')
+      !
+      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 
+      CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs ) 
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp)WRITE(numout,*)
+         IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps'
+      ENDIF
+
+      IF( l_trddyn )   THEN         ! Save ua and va trends
+         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
+         ztrdu(:,:,:) = ua(:,:,:) 
+         ztrdv(:,:,:) = va(:,:,:) 
+      ENDIF
+      
+      IF( neuler == 0 .AND. kt == nit000 ) THEN
+          z2dtzts =         rdt / REAL( jnzts, wp )   ! = rdt (restart with Euler time stepping)
+      ELSE
+          z2dtzts = 2._wp * rdt / REAL( jnzts, wp )   ! = 2 rdt (leapfrog)
+      ENDIF
+      
+      DO jk = 2, jpkm1                    ! Calculate and store vertical fluxes
+         DO jj = 2, jpj                   
+            DO ji = fs_2, jpi             ! vector opt.
+               zww(ji,jj,jk) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      DO jj = 2, jpjm1                    ! Surface and bottom advective fluxes set to zero
+         DO ji = fs_2, fs_jpim1           ! vector opt.
+            zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp
+            zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp
+            zwuw(ji,jj,jpk) = 0._wp
+            zwvw(ji,jj,jpk) = 0._wp
+         END DO  
+      END DO
+
+! Start with before values and use sub timestepping to reach after values
+
+      zus(:,:,:,1) = ub(:,:,:)
+      zvs(:,:,:,1) = vb(:,:,:)
+
+      DO jl = 1, jnzts                   ! Start of sub timestepping loop
+
+         IF( jl == 1 ) THEN              ! Euler forward to kick things off
+           jtb = 1   ;   jtn = 1   ;   jta = 2
+           zts = z2dtzts
+         ELSEIF( jl == 2 ) THEN          ! First leapfrog step
+           jtb = 1   ;   jtn = 2   ;   jta = 3
+           zts = 2._wp * z2dtzts
+         ELSE                            ! Shuffle pointers for subsequent leapfrog steps
+           jtb = MOD(jtb,3) + 1
+           jtn = MOD(jtn,3) + 1
+           jta = MOD(jta,3) + 1
+         ENDIF
+
+         DO jk = 2, jpkm1           ! Vertical momentum advection at level w and u- and v- vertical
+            DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
+               DO ji = fs_2, fs_jpim1        ! vector opt.
+                  zwuw(ji,jj,jk) = ( zww(ji+1,jj  ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) )
+                  zwvw(ji,jj,jk) = ( zww(ji  ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) )
+               END DO  
+            END DO   
+         END DO
+         DO jk = 1, jpkm1           ! Vertical momentum advection at u- and v-points
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1       ! vector opt.
+                  !                         ! vertical momentum advective trends
+                  zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
+                  zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
+                  zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts
+                  zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts
+               END DO  
+            END DO  
+         END DO
+
+      END DO      ! End of sub timestepping loop
+
+      zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts )
+      DO jk = 1, jpkm1              ! Recover trends over the outer timestep
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1       ! vector opt.
+               !                         ! vertical momentum advective trends
+               !                         ! add the trends to the general momentum trends
+               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt
+               va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt
+            END DO  
+         END DO  
+      END DO
+
+      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
+         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         CALL trd_mod(ztrdu, ztrdv, jpdyn_trd_zad, 'DYN', kt)
+         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
+      ENDIF
+      !                             ! Control print
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
+         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 
+      CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad_zts')
+      !
+   END SUBROUTINE dyn_zad_zts
+
    !!======================================================================
 END MODULE dynzad
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icb_oce.F90
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icb_oce.F90	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icb_oce.F90	(revision 4934)
@@ -96,5 +96,5 @@
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ua_e, va_e
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_e
-#if defined key_lim2 || defined key_lim3
+#if defined key_lim2 || defined key_lim3 || defined key_cice
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ui_e, vi_e
 #endif
@@ -144,4 +144,7 @@
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)       ::   nicbflddest                      !: nfold destination proc
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)       ::   nicbfldproc                      !: nfold destination proc
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)       ::   nicbfldnsend                     !: nfold number of bergs to send to nfold neighbour
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)       ::   nicbfldexpect                    !: nfold expected number of bergs
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)       ::   nicbfldreq                       !: nfold message handle (immediate send)
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: griddata                           !: work array for icbrst
@@ -162,4 +165,6 @@
       !
       icb_alloc = 0
+      ALLOCATE( berg_grid, STAT=ill )
+      icb_alloc = icb_alloc + ill
       ALLOCATE( berg_grid%calving    (jpi,jpj) , berg_grid%calving_hflx (jpi,jpj)          ,   &
          &      berg_grid%stored_heat(jpi,jpj) , berg_grid%floating_melt(jpi,jpj)          ,   &
@@ -171,5 +176,5 @@
       ALLOCATE( uo_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) ,   &
          &      vo_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) ,   &
-#if defined key_lim2 || defined key_lim3
+#if defined key_lim2 || defined key_lim3 || defined key_cice
          &      ui_e(0:jpi+1,0:jpj+1) ,                            &
          &      vi_e(0:jpi+1,0:jpj+1) ,                            &
@@ -181,5 +186,6 @@
       icb_alloc = icb_alloc + ill
 
-      ALLOCATE( nicbfldpts(jpi) , nicbflddest(jpi) , nicbfldproc(jpni) , STAT=ill)
+      ALLOCATE( nicbfldpts(jpi) , nicbflddest(jpi) , nicbfldproc(jpni) , &
+         &      nicbfldnsend(jpni), nicbfldexpect(jpni) , nicbfldreq(jpni), STAT=ill)
       icb_alloc = icb_alloc + ill
 
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90	(revision 4934)
@@ -33,5 +33,5 @@
 CONTAINS
 
-   SUBROUTINE icb_dyn()
+   SUBROUTINE icb_dyn( kt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_dyn  ***
@@ -50,4 +50,5 @@
       TYPE(iceberg), POINTER          ::   berg
       TYPE(point)  , POINTER          ::   pt
+      INTEGER                         ::   kt
       !!----------------------------------------------------------------------
 
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90	(revision 4934)
@@ -172,13 +172,16 @@
          DO ji = nicbdi, nicbei
             ii = nicbflddest(ji)
-            DO jn = 1, jpni
-               ! work along array until we find an empty slot
-               IF( nicbfldproc(jn) == -1 ) THEN
-                  nicbfldproc(jn) = ii
-                  EXIT                             !!gm EXIT should be avoided: use DO WHILE expression instead
-               ENDIF
-               ! before we find an empty slot, we may find processor number is already here so we exit
-               IF( nicbfldproc(jn) == ii ) EXIT
-            END DO
+            IF( ii .GT. 0 ) THEN     ! Needed because land suppression can mean
+                                     ! that unused points are not set in edge haloes
+               DO jn = 1, jpni
+                  ! work along array until we find an empty slot
+                  IF( nicbfldproc(jn) == -1 ) THEN
+                     nicbfldproc(jn) = ii
+                     EXIT                             !!gm EXIT should be avoided: use DO WHILE expression instead
+                  ENDIF
+                  ! before we find an empty slot, we may find processor number is already here so we exit
+                  IF( nicbfldproc(jn) == ii ) EXIT
+               END DO
+            ENDIF
          END DO
       ENDIF
@@ -210,4 +213,6 @@
             WRITE(numicb,*) 'north fold destination procs  '
             WRITE(numicb,*) nicbflddest
+            WRITE(numicb,*) 'north fold destination proclist  '
+            WRITE(numicb,*) nicbfldproc
          ENDIF
          CALL flush(numicb)
@@ -397,7 +402,7 @@
       ENDIF
 
-      IF( lk_lim3 .AND. ln_icebergs ) THEN
-         CALL ctl_stop( 'icb_nam: the use of ICB with LIM3 not allowed. ice thickness missing in ICB' )
-      ENDIF
+!     IF( lk_lim3 .AND. ln_icebergs ) THEN
+!        CALL ctl_stop( 'icb_nam: the use of ICB with LIM3 not allowed. ice thickness missing in ICB' )
+!     ENDIF
 
       IF(lwp) THEN                  ! control print
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icblbc.F90
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icblbc.F90	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icblbc.F90	(revision 4934)
@@ -280,5 +280,5 @@
          zwebergs(1) = ibergs_to_send_e
          CALL mppsend( 12, zwebergs(1), 1, ipe_E, iml_req1)
-         CALL mpprecv( 11, zewbergs(2), 1 )
+         CALL mpprecv( 11, zewbergs(2), 1, ipe_E )
          IF( l_isend ) CALL mpi_wait( iml_req1, iml_stat, iml_err )
          ibergs_rcvd_from_e = INT( zewbergs(2) )
@@ -288,6 +288,6 @@
          CALL mppsend( 11, zewbergs(1), 1, ipe_W, iml_req2)
          CALL mppsend( 12, zwebergs(1), 1, ipe_E, iml_req3)
-         CALL mpprecv( 11, zewbergs(2), 1 )
-         CALL mpprecv( 12, zwebergs(2), 1 )
+         CALL mpprecv( 11, zewbergs(2), 1, ipe_E )
+         CALL mpprecv( 12, zwebergs(2), 1, ipe_W )
          IF( l_isend ) CALL mpi_wait( iml_req2, iml_stat, iml_err )
          IF( l_isend ) CALL mpi_wait( iml_req3, iml_stat, iml_err )
@@ -297,5 +297,5 @@
          zewbergs(1) = ibergs_to_send_w
          CALL mppsend( 11, zewbergs(1), 1, ipe_W, iml_req4)
-         CALL mpprecv( 12, zwebergs(2), 1 )
+         CALL mpprecv( 12, zwebergs(2), 1, ipe_W )
          IF( l_isend ) CALL mpi_wait( iml_req4, iml_stat, iml_err )
          ibergs_rcvd_from_w = INT( zwebergs(2) )
@@ -411,5 +411,5 @@
          zsnbergs(1) = ibergs_to_send_n
          CALL mppsend( 16, zsnbergs(1), 1, ipe_N, iml_req1)
-         CALL mpprecv( 15, znsbergs(2), 1 )
+         CALL mpprecv( 15, znsbergs(2), 1, ipe_N )
          IF( l_isend ) CALL mpi_wait( iml_req1, iml_stat, iml_err )
          ibergs_rcvd_from_n = INT( znsbergs(2) )
@@ -419,6 +419,6 @@
          CALL mppsend( 15, znsbergs(1), 1, ipe_S, iml_req2)
          CALL mppsend( 16, zsnbergs(1), 1, ipe_N, iml_req3)
-         CALL mpprecv( 15, znsbergs(2), 1 )
-         CALL mpprecv( 16, zsnbergs(2), 1 )
+         CALL mpprecv( 15, znsbergs(2), 1, ipe_N )
+         CALL mpprecv( 16, zsnbergs(2), 1, ipe_S )
          IF( l_isend ) CALL mpi_wait( iml_req2, iml_stat, iml_err )
          IF( l_isend ) CALL mpi_wait( iml_req3, iml_stat, iml_err )
@@ -428,5 +428,5 @@
          znsbergs(1) = ibergs_to_send_s
          CALL mppsend( 15, znsbergs(1), 1, ipe_S, iml_req4)
-         CALL mpprecv( 16, zsnbergs(2), 1 )
+         CALL mpprecv( 16, zsnbergs(2), 1, ipe_S )
          IF( l_isend ) CALL mpi_wait( iml_req4, iml_stat, iml_err )
          ibergs_rcvd_from_s = INT( zsnbergs(2) )
@@ -581,5 +581,6 @@
       INTEGER                             :: ifldproc, iproc, ipts
       INTEGER                             :: iine, ijne
-      REAL(wp), DIMENSION(2)              :: zsbergs, znbergs
+      INTEGER                             :: jjn
+      REAL(wp), DIMENSION(0:3)            :: zsbergs, znbergs
       INTEGER                             :: iml_req1, iml_req2, iml_err
       INTEGER, DIMENSION(MPI_STATUS_SIZE) :: iml_stat
@@ -591,4 +592,93 @@
       ! its of fixed size, the first -1 marks end of list of processors
       !
+      nicbfldnsend(:) = 0
+      nicbfldexpect(:) = 0
+      nicbfldreq(:) = 0
+      !
+      ! Since each processor may be communicating with more than one northern
+      ! neighbour, cycle through the sends so that the receive order can be
+      ! controlled.
+      !
+      ! First compute how many icebergs each active neighbour should expect
+      DO jn = 1, jpni
+         IF( nicbfldproc(jn) /= -1 ) THEN
+            ifldproc = nicbfldproc(jn)
+            nicbfldnsend(jn) = 0
+
+            ! Find number of bergs that need to be exchanged
+            ! Pick out exchanges with processor ifldproc
+            ! if ifldproc is this processor then don't send
+            !
+            IF( ASSOCIATED(first_berg) ) THEN
+               this => first_berg
+               DO WHILE (ASSOCIATED(this))
+                  pt => this%current_point
+                  iine = INT( pt%xi + 0.5 )
+                  ijne = INT( pt%yj + 0.5 )
+                  iproc = nicbflddest(mi1(iine))
+                  IF( ijne .GT. mjg(nicbej) ) THEN
+                     IF( iproc == ifldproc ) THEN
+                        !
+                        IF( iproc /= narea ) THEN
+                           tmpberg => this
+                           nicbfldnsend(jn) = nicbfldnsend(jn) + 1
+                        ENDIF
+                        !
+                     ENDIF
+                  ENDIF
+                  this => this%next
+               END DO
+            ENDIF
+            !
+         ENDIF
+         !
+      END DO
+      !
+      ! Now tell each active neighbour how many icebergs to expect
+      DO jn = 1, jpni
+         IF( nicbfldproc(jn) /= -1 ) THEN
+            ifldproc = nicbfldproc(jn)
+            IF( ifldproc == narea ) CYCLE
+   
+            zsbergs(0) = narea
+            zsbergs(1) = nicbfldnsend(jn)
+            !IF ( nicbfldnsend(jn) .GT. 0) write(numicb,*) 'ICB sending ',nicbfldnsend(jn),' to ', ifldproc
+            CALL mppsend( 21, zsbergs(0:1), 2, ifldproc-1, nicbfldreq(jn))
+         ENDIF
+         !
+      END DO
+      !
+      ! and receive the heads-up from active neighbours preparing to send
+      DO jn = 1, jpni
+         IF( nicbfldproc(jn) /= -1 ) THEN
+            ifldproc = nicbfldproc(jn)
+            IF( ifldproc == narea ) CYCLE
+
+            CALL mpprecv( 21, znbergs(1:2), 2 )
+            DO jjn = 1,jpni
+             IF( nicbfldproc(jjn) .eq. INT(znbergs(1)) ) EXIT
+            END DO
+            IF( jjn .GT. jpni ) write(numicb,*) 'ICB ERROR'
+            nicbfldexpect(jjn) = INT( znbergs(2) )
+            !IF ( nicbfldexpect(jjn) .GT. 0) write(numicb,*) 'ICB expecting ',nicbfldexpect(jjn),' from ', nicbfldproc(jjn)
+            !CALL FLUSH(numicb)
+         ENDIF
+         !
+      END DO
+      !
+      ! post the mpi waits if using immediate send protocol
+      DO jn = 1, jpni
+         IF( nicbfldproc(jn) /= -1 ) THEN
+            ifldproc = nicbfldproc(jn)
+            IF( ifldproc == narea ) CYCLE
+
+            IF( l_isend ) CALL mpi_wait( nicbfldreq(jn), iml_stat, iml_err )
+         ENDIF
+         !
+      END DO
+   
+         !
+         ! Cycle through the icebergs again, this time packing and sending any
+         ! going through the north fold. They will be expected.
       DO jn = 1, jpni
          IF( nicbfldproc(jn) /= -1 ) THEN
@@ -646,19 +736,25 @@
             IF( ifldproc == narea ) CYCLE
    
-            zsbergs(1) = ibergs_to_send
-            CALL mppsend( 21, zsbergs(1), 1, ifldproc-1, iml_req1)
-            CALL mpprecv( 21, znbergs(2), 1 )
-            IF( l_isend ) CALL mpi_wait( iml_req1, iml_stat, iml_err )
-            ibergs_to_rcv = INT( znbergs(2) )
-   
             ! send bergs
    
             IF( ibergs_to_send > 0 )  &
-                CALL mppsend( 12, obuffer_f%data, ibergs_to_send*jp_buffer_width, ifldproc-1, iml_req2 )
+                CALL mppsend( 12, obuffer_f%data, ibergs_to_send*jp_buffer_width, ifldproc-1, nicbfldreq(jn) )
+            !
+         ENDIF
+         !
+      END DO
+      !
+      ! Now receive the expected number of bergs from the active neighbours
+      DO jn = 1, jpni
+         IF( nicbfldproc(jn) /= -1 ) THEN
+            ifldproc = nicbfldproc(jn)
+            IF( ifldproc == narea ) CYCLE
+            ibergs_to_rcv = nicbfldexpect(jn)
+
             IF( ibergs_to_rcv  > 0 ) THEN
                CALL icb_increase_ibuffer(ibuffer_f, ibergs_to_rcv)
-               CALL mpprecv( 12, ibuffer_f%data, ibergs_to_rcv*jp_buffer_width )
-            ENDIF
-            IF( ibergs_to_send > 0 .AND. l_isend ) CALL mpi_wait( iml_req2, iml_stat, iml_err )
+               CALL mpprecv( 12, ibuffer_f%data, ibergs_to_rcv*jp_buffer_width, ifldproc-1 )
+            ENDIF
+            !
             DO jk = 1, ibergs_to_rcv
                IF( nn_verbose_level >= 4 ) THEN
@@ -668,5 +764,15 @@
                CALL icb_unpack_from_buffer(first_berg, ibuffer_f, jk )
             END DO
-            !
+         ENDIF
+         !
+      END DO
+      !
+      ! Finally post the mpi waits if using immediate send protocol
+      DO jn = 1, jpni
+         IF( nicbfldproc(jn) /= -1 ) THEN
+            ifldproc = nicbfldproc(jn)
+            IF( ifldproc == narea ) CYCLE
+
+            IF( l_isend ) CALL mpi_wait( nicbfldreq(jn), iml_stat, iml_err )
          ENDIF
          !
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 4934)
@@ -64,5 +64,5 @@
                                                                                             ! start and count arrays
       LOGICAL                      ::   ll_found_restart
-      CHARACTER(len=80)            ::   cl_filename
+      CHARACTER(len=256)           ::   cl_filename
       CHARACTER(len=NF90_MAX_NAME) ::   cl_dname
       TYPE(iceberg)                ::   localberg ! NOT a pointer but an actual local variable
@@ -228,5 +228,5 @@
       INTEGER ::   jn   ! dummy loop index
       INTEGER ::   ix_dim, iy_dim, ik_dim, in_dim
-      CHARACTER(len=80)      :: cl_filename
+      CHARACTER(len=256)     :: cl_filename
       TYPE(iceberg), POINTER :: this
       TYPE(point)  , POINTER :: pt
@@ -256,4 +256,19 @@
       IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim k failed')
 
+      ! global attributes
+      IF( lk_mpp ) THEN
+         ! Set domain parameters (assume jpdom_local_full)
+         nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_number_total'   , jpnij              )
+         nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_number'         , narea-1            )
+         nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_dimensions_ids' , (/1     , 2     /) )
+         nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_size_global'    , (/jpiglo, jpjglo/) )
+         nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_size_local'     , (/jpi   , jpj   /) )
+         nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_position_first' , (/nimpp , njmpp /) )
+         nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_position_last'  , (/nimpp + jpi - 1 , njmpp + jpj - 1  /) )
+         nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_halo_size_start', (/nldi - 1        , nldj - 1         /) )
+         nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_halo_size_end'  , (/jpi - nlei      , jpj - nlej       /) )
+         nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_type'           , 'BOX'              )
+      ENDIF
+      
       IF (associated(first_berg)) then
          nret = NF90_DEF_DIM(ncid, 'n', NF90_UNLIMITED, in_dim)
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbstp.F90
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbstp.F90	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbstp.F90	(revision 4934)
@@ -105,5 +105,5 @@
 !                               !==  For each berg, evolve  ==!
       !
-      IF( ASSOCIATED(first_berg) )   CALL icb_dyn()           ! ice berg dynamics
+      IF( ASSOCIATED(first_berg) )   CALL icb_dyn( kt )       ! ice berg dynamics
 
       IF( lk_mpp ) THEN          ;   CALL icb_lbc_mpp()       ! Send bergs to other PEs
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90	(revision 4934)
@@ -76,9 +76,9 @@
       va_e(:,:) = 0._wp ;   va_e(1:jpi, 1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk
 
-      CALL lbc_lnk_e( uo_e, 'U', -1._wp, 1, 1 )
-      CALL lbc_lnk_e( vo_e, 'V', -1._wp, 1, 1 )
-      CALL lbc_lnk_e( ff_e, 'F', +1._wp, 1, 1 )
-      CALL lbc_lnk_e( ua_e, 'U', -1._wp, 1, 1 )
-      CALL lbc_lnk_e( va_e, 'V', -1._wp, 1, 1 )
+      CALL lbc_lnk_icb( uo_e, 'U', -1._wp, 1, 1 )
+      CALL lbc_lnk_icb( vo_e, 'V', -1._wp, 1, 1 )
+      CALL lbc_lnk_icb( ff_e, 'F', +1._wp, 1, 1 )
+      CALL lbc_lnk_icb( ua_e, 'U', -1._wp, 1, 1 )
+      CALL lbc_lnk_icb( va_e, 'V', -1._wp, 1, 1 )
 
 #if defined key_lim2 || defined key_lim3
@@ -86,6 +86,6 @@
       vi_e(:,:) = 0._wp ;   vi_e(1:jpi, 1:jpj) = v_ice(:,:)
 
-      CALL lbc_lnk_e( ui_e, 'U', -1._wp, 1, 1 )
-      CALL lbc_lnk_e( vi_e, 'V', -1._wp, 1, 1 )
+      CALL lbc_lnk_icb( ui_e, 'U', -1._wp, 1, 1 )
+      CALL lbc_lnk_icb( vi_e, 'V', -1._wp, 1, 1 )
 #endif
 
@@ -102,5 +102,5 @@
       ssh_e(0,jpj+1)     = ssh_e(1,jpj)
       ssh_e(jpi+1,jpj+1) = ssh_e(jpi,jpj)
-      CALL lbc_lnk_e( ssh_e, 'T', +1._wp, 1, 1 )
+      CALL lbc_lnk_icb( ssh_e, 'T', +1._wp, 1, 1 )
       !
    END SUBROUTINE icb_utl_copy
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90	(revision 4934)
@@ -34,7 +34,12 @@
    END INTERFACE
 
+   INTERFACE lbc_lnk_icb
+      MODULE PROCEDURE mpp_lnk_2d_icb
+   END INTERFACE
+
    PUBLIC lbc_lnk       ! ocean lateral boundary conditions
    PUBLIC lbc_lnk_e
    PUBLIC lbc_bdy_lnk   ! ocean lateral BDY boundary conditions
+   PUBLIC lbc_lnk_icb
 
    !!----------------------------------------------------------------------
@@ -73,7 +78,12 @@
    END INTERFACE
 
+   INTERFACE lbc_lnk_icb
+      MODULE PROCEDURE lbc_lnk_2d_e
+   END INTERFACE
+
    PUBLIC   lbc_lnk       ! ocean/ice  lateral boundary conditions
    PUBLIC   lbc_lnk_e 
    PUBLIC   lbc_bdy_lnk   ! ocean lateral BDY boundary conditions
+   PUBLIC   lbc_lnk_icb
    
    !!----------------------------------------------------------------------
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 4934)
@@ -42,4 +42,5 @@
    !!   mpp_lnk_3d_gather :  Message passing manadgement for two 3D arrays
    !!   mpp_lnk_e     : interface (defined in lbclnk) for message passing of 2d array with extra halo (mpp_lnk_2d_e)
+   !!   mpp_lnk_icb   : interface for message passing of 2d arrays with extra halo for icebergs (mpp_lnk_2d_icb)
    !!   mpprecv         :
    !!   mppsend       :   SUBROUTINE mpp_ini_znl
@@ -56,4 +57,5 @@
    !!   mpp_lbc_north : north fold processors gathering
    !!   mpp_lbc_north_e : variant of mpp_lbc_north for extra outer halo
+   !!   mpp_lbc_north_icb : variant of mpp_lbc_north for extra outer halo with icebergs
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
@@ -74,4 +76,5 @@
    PUBLIC   mppsend, mpprecv                          ! needed by TAM and ICB routines
    PUBLIC   mpp_lnk_bdy_2d, mpp_lnk_bdy_3d
+   PUBLIC   mpp_lbc_north_icb, mpp_lnk_2d_icb
 
    !! * Interfaces
@@ -2893,4 +2896,267 @@
    END SUBROUTINE DDPDD_MPI
 
+   SUBROUTINE mpp_lbc_north_icb( pt2d, cd_type, psgn, pr2dj)
+      !!---------------------------------------------------------------------
+      !!                   ***  routine mpp_lbc_north_icb  ***
+      !!
+      !! ** Purpose :   Ensure proper north fold horizontal bondary condition
+      !!              in mpp configuration in case of jpn1 > 1 and for 2d
+      !!              array with outer extra halo
+      !!
+      !! ** Method  :   North fold condition and mpp with more than one proc
+      !!              in i-direction require a specific treatment. We gather
+      !!              the 4+2*jpr2dj northern lines of the global domain on 1
+      !!              processor and apply lbc north-fold on this sub array.
+      !!              Then we scatter the north fold array back to the processors.
+      !!              This version accounts for an extra halo with icebergs.
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt2d     ! 2D array with extra halo
+      CHARACTER(len=1)        , INTENT(in   ) ::   cd_type  ! nature of pt3d grid-points
+      !                                                     !   = T ,  U , V , F or W -points
+      REAL(wp)                , INTENT(in   ) ::   psgn     ! = -1. the sign change across the
+      !!                                                    ! north fold, =  1. otherwise
+      INTEGER, OPTIONAL       , INTENT(in   ) ::   pr2dj
+      INTEGER ::   ji, jj, jr
+      INTEGER ::   ierr, itaille, ildi, ilei, iilb
+      INTEGER ::   ijpj, ij, iproc, ipr2dj
+      !
+      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE  ::  ztab_e, znorthloc_e
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE  ::  znorthgloio_e
+
+      !!----------------------------------------------------------------------
+      !
+      ijpj=4
+      IF( PRESENT(pr2dj) ) THEN           ! use of additional halos
+         ipr2dj = pr2dj
+      ELSE
+         ipr2dj = 0
+      ENDIF
+      ALLOCATE( ztab_e(jpiglo,4+2*ipr2dj), znorthloc_e(jpi,4+2*ipr2dj), znorthgloio_e(jpi,4+2*ipr2dj,jpni) )
+
+      !
+      ztab_e(:,:) = 0.e0
+
+      ij=0
+      ! put in znorthloc_e the last 4 jlines of pt2d
+      DO jj = nlcj - ijpj + 1 - ipr2dj, nlcj +ipr2dj
+         ij = ij + 1
+         DO ji = 1, jpi
+            znorthloc_e(ji,ij)=pt2d(ji,jj)
+         END DO
+      END DO
+      !
+      itaille = jpi * ( ijpj + 2 * ipr2dj )
+      CALL MPI_ALLGATHER( znorthloc_e(1,1)  , itaille, MPI_DOUBLE_PRECISION,    &
+         &                znorthgloio_e(1,1,1), itaille, MPI_DOUBLE_PRECISION, ncomm_north, ierr )
+      !
+      DO jr = 1, ndim_rank_north            ! recover the global north array
+         iproc = nrank_north(jr) + 1
+         ildi = nldit (iproc)
+         ilei = nleit (iproc)
+         iilb = nimppt(iproc)
+         DO jj = 1, ijpj+2*ipr2dj
+            DO ji = ildi, ilei
+               ztab_e(ji+iilb-1,jj) = znorthgloio_e(ji,jj,jr)
+            END DO
+         END DO
+      END DO
+
+
+      ! 2. North-Fold boundary conditions
+      ! ----------------------------------
+      CALL lbc_nfd( ztab_e(:,:), cd_type, psgn, pr2dj = ipr2dj )
+
+      ij = ipr2dj
+      !! Scatter back to pt2d
+      DO jj = nlcj - ijpj + 1 , nlcj +ipr2dj
+      ij  = ij +1
+         DO ji= 1, nlci
+            pt2d(ji,jj) = ztab_e(ji+nimpp-1,ij)
+         END DO
+      END DO
+      !
+      DEALLOCATE( ztab_e, znorthloc_e, znorthgloio_e )
+      !
+   END SUBROUTINE mpp_lbc_north_icb
+
+   SUBROUTINE mpp_lnk_2d_icb( pt2d, cd_type, psgn, jpri, jprj )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mpp_lnk_2d_icb  ***
+      !!
+      !! ** Purpose :   Message passing manadgement for 2d array (with extra halo and icebergs)
+      !!
+      !! ** Method  :   Use mppsend and mpprecv function for passing mask
+      !!      between processors following neighboring subdomains.
+      !!            domain parameters
+      !!                    nlci   : first dimension of the local subdomain
+      !!                    nlcj   : second dimension of the local subdomain
+      !!                    jpri   : number of rows for extra outer halo
+      !!                    jprj   : number of columns for extra outer halo
+      !!                    nbondi : mark for "east-west local boundary"
+      !!                    nbondj : mark for "north-south local boundary"
+      !!                    noea   : number for local neighboring processors
+      !!                    nowe   : number for local neighboring processors
+      !!                    noso   : number for local neighboring processors
+      !!                    nono   : number for local neighboring processors
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                                             , INTENT(in   ) ::   jpri
+      INTEGER                                             , INTENT(in   ) ::   jprj
+      REAL(wp), DIMENSION(1-jpri:jpi+jpri,1-jprj:jpj+jprj), INTENT(inout) ::   pt2d     ! 2D array with extra halo
+      CHARACTER(len=1)                                    , INTENT(in   ) ::   cd_type  ! nature of ptab array grid-points
+      !                                                                                 ! = T , U , V , F , W and I points
+      REAL(wp)                                            , INTENT(in   ) ::   psgn     ! =-1 the sign change across the
+      !!                                                                                ! north boundary, =  1. otherwise
+      INTEGER  ::   jl   ! dummy loop indices
+      INTEGER  ::   imigr, iihom, ijhom        ! temporary integers
+      INTEGER  ::   ipreci, iprecj             ! temporary integers
+      INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend
+      INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   ml_stat   ! for key_mpi_isend
+      !!
+      REAL(wp), DIMENSION(1-jpri:jpi+jpri,jprecj+jprj,2) :: r2dns
+      REAL(wp), DIMENSION(1-jpri:jpi+jpri,jprecj+jprj,2) :: r2dsn
+      REAL(wp), DIMENSION(1-jprj:jpj+jprj,jpreci+jpri,2) :: r2dwe
+      REAL(wp), DIMENSION(1-jprj:jpj+jprj,jpreci+jpri,2) :: r2dew
+      !!----------------------------------------------------------------------
+
+      ipreci = jpreci + jpri      ! take into account outer extra 2D overlap area
+      iprecj = jprecj + jprj
+
+
+      ! 1. standard boundary treatment
+      ! ------------------------------
+      ! Order matters Here !!!!
+      !
+      !                                      ! East-West boundaries
+      !                                           !* Cyclic east-west
+      IF( nbondi == 2 .AND. (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN
+         pt2d(1-jpri:     1    ,:) = pt2d(jpim1-jpri:  jpim1 ,:)       ! east
+         pt2d(   jpi  :jpi+jpri,:) = pt2d(     2      :2+jpri,:)       ! west
+         !
+      ELSE                                        !* closed
+         IF( .NOT. cd_type == 'F' )   pt2d(  1-jpri   :jpreci    ,:) = 0.e0    ! south except at F-point
+                                      pt2d(nlci-jpreci+1:jpi+jpri,:) = 0.e0    ! north
+      ENDIF
+      !
+
+      ! north fold treatment
+      ! -----------------------
+      IF( npolj /= 0 ) THEN
+         !
+         SELECT CASE ( jpni )
+         CASE ( 1 )     ;   CALL lbc_nfd        ( pt2d(1:jpi,1:jpj+jprj), cd_type, psgn, pr2dj=jprj )
+         CASE DEFAULT   ;   CALL mpp_lbc_north_icb( pt2d(1:jpi,1:jpj+jprj)  , cd_type, psgn , pr2dj=jprj  )
+         END SELECT
+         !
+      ENDIF
+
+      ! 2. East and west directions exchange
+      ! ------------------------------------
+      ! we play with the neigbours AND the row number because of the periodicity
+      !
+      SELECT CASE ( nbondi )      ! Read Dirichlet lateral conditions
+      CASE ( -1, 0, 1 )                ! all exept 2 (i.e. close case)
+         iihom = nlci-nreci-jpri
+         DO jl = 1, ipreci
+            r2dew(:,jl,1) = pt2d(jpreci+jl,:)
+            r2dwe(:,jl,1) = pt2d(iihom +jl,:)
+         END DO
+      END SELECT
+      !
+      !                           ! Migrations
+      imigr = ipreci * ( jpj + 2*jprj)
+      !
+      SELECT CASE ( nbondi )
+      CASE ( -1 )
+         CALL mppsend( 2, r2dwe(1-jprj,1,1), imigr, noea, ml_req1 )
+         CALL mpprecv( 1, r2dew(1-jprj,1,2), imigr, noea )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+      CASE ( 0 )
+         CALL mppsend( 1, r2dew(1-jprj,1,1), imigr, nowe, ml_req1 )
+         CALL mppsend( 2, r2dwe(1-jprj,1,1), imigr, noea, ml_req2 )
+         CALL mpprecv( 1, r2dew(1-jprj,1,2), imigr, noea )
+         CALL mpprecv( 2, r2dwe(1-jprj,1,2), imigr, nowe )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+         IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
+      CASE ( 1 )
+         CALL mppsend( 1, r2dew(1-jprj,1,1), imigr, nowe, ml_req1 )
+         CALL mpprecv( 2, r2dwe(1-jprj,1,2), imigr, nowe )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+      END SELECT
+      !
+      !                           ! Write Dirichlet lateral conditions
+      iihom = nlci - jpreci
+      !
+      SELECT CASE ( nbondi )
+      CASE ( -1 )
+         DO jl = 1, ipreci
+            pt2d(iihom+jl,:) = r2dew(:,jl,2)
+         END DO
+      CASE ( 0 )
+         DO jl = 1, ipreci
+            pt2d(jl-jpri,:) = r2dwe(:,jl,2)
+            pt2d( iihom+jl,:) = r2dew(:,jl,2)
+         END DO
+      CASE ( 1 )
+         DO jl = 1, ipreci
+            pt2d(jl-jpri,:) = r2dwe(:,jl,2)
+         END DO
+      END SELECT
+
+
+      ! 3. North and south directions
+      ! -----------------------------
+      ! always closed : we play only with the neigbours
+      !
+      IF( nbondj /= 2 ) THEN      ! Read Dirichlet lateral conditions
+         ijhom = nlcj-nrecj-jprj
+         DO jl = 1, iprecj
+            r2dsn(:,jl,1) = pt2d(:,ijhom +jl)
+            r2dns(:,jl,1) = pt2d(:,jprecj+jl)
+         END DO
+      ENDIF
+      !
+      !                           ! Migrations
+      imigr = iprecj * ( jpi + 2*jpri )
+      !
+      SELECT CASE ( nbondj )
+      CASE ( -1 )
+         CALL mppsend( 4, r2dsn(1-jpri,1,1), imigr, nono, ml_req1 )
+         CALL mpprecv( 3, r2dns(1-jpri,1,2), imigr, nono )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+      CASE ( 0 )
+         CALL mppsend( 3, r2dns(1-jpri,1,1), imigr, noso, ml_req1 )
+         CALL mppsend( 4, r2dsn(1-jpri,1,1), imigr, nono, ml_req2 )
+         CALL mpprecv( 3, r2dns(1-jpri,1,2), imigr, nono )
+         CALL mpprecv( 4, r2dsn(1-jpri,1,2), imigr, noso )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+         IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
+      CASE ( 1 )
+         CALL mppsend( 3, r2dns(1-jpri,1,1), imigr, noso, ml_req1 )
+         CALL mpprecv( 4, r2dsn(1-jpri,1,2), imigr, noso )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+      END SELECT
+      !
+      !                           ! Write Dirichlet lateral conditions
+      ijhom = nlcj - jprecj
+      !
+      SELECT CASE ( nbondj )
+      CASE ( -1 )
+         DO jl = 1, iprecj
+            pt2d(:,ijhom+jl) = r2dns(:,jl,2)
+         END DO
+      CASE ( 0 )
+         DO jl = 1, iprecj
+            pt2d(:,jl-jprj) = r2dsn(:,jl,2)
+            pt2d(:,ijhom+jl ) = r2dns(:,jl,2)
+         END DO
+      CASE ( 1 )
+         DO jl = 1, iprecj
+            pt2d(:,jl-jprj) = r2dsn(:,jl,2)
+         END DO
+      END SELECT
+
+   END SUBROUTINE mpp_lnk_2d_icb
 #else
    !!----------------------------------------------------------------------
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 4934)
@@ -208,29 +208,6 @@
       !-----------------------------------------------------------------------
 
-      !Initalise all values in namelist arrays
-      enactfiles(:) = ''
-      coriofiles(:) = ''
-      profbfiles(:) = ''
-      slafilesact(:) = ''
-      slafilespas(:) = ''
-      slafbfiles(:) = ''
-      sstfiles(:)   = ''
-      sstfbfiles(:) = ''
-      seaicefiles(:) = ''
       velcurfiles(:) = ''
       veladcpfiles(:) = ''
-      velavcurfiles(:) = ''
-      velhrcurfiles(:) = ''
-      velavadcpfiles(:) = ''
-      velhradcpfiles(:) = ''
-      velfbfiles(:) = ''
-      velcurfiles(:) = ''
-      veladcpfiles(:) = ''
-      endailyavtypes(:) = -1
-      endailyavtypes(1) = 820
-      ln_profb_ena(:) = .FALSE.
-      ln_profb_enatim(:) = .TRUE.
-      ln_velfb_av(:) = .FALSE.
-      ln_ignmis = .FALSE.
       CALL ini_date( dobsini )
       CALL fin_date( dobsend )
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 4934)
@@ -44,4 +44,5 @@
    LOGICAL ::   ln_traadv_cen2     ! 2nd order centered scheme flag
    LOGICAL ::   ln_traadv_tvd      ! TVD scheme flag
+   LOGICAL ::   ln_traadv_tvd_zts  ! TVD scheme flag with vertical sub time-stepping
    LOGICAL ::   ln_traadv_muscl    ! MUSCL scheme flag
    LOGICAL ::   ln_traadv_muscl2   ! MUSCL2 scheme flag
@@ -121,4 +122,5 @@
       CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS 
       CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  QUICKEST 
+      CASE ( 7 )   ;   CALL tra_adv_tvd_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD ZTS
       !
       CASE (-1 )                                      !==  esopa: test all possibility with control print  ==!
@@ -167,5 +169,5 @@
          &                 ln_traadv_muscl, ln_traadv_muscl2,  &
          &                 ln_traadv_ubs  , ln_traadv_qck,     &
-         &                 ln_traadv_msc_ups
+         &                 ln_traadv_msc_ups, ln_traadv_tvd_zts
       !!----------------------------------------------------------------------
 
@@ -191,4 +193,5 @@
          WRITE(numout,*) '      QUICKEST advection scheme      ln_traadv_qck     = ', ln_traadv_qck
          WRITE(numout,*) '      upstream scheme within muscl   ln_traadv_msc_ups = ', ln_traadv_msc_ups
+         WRITE(numout,*) '      TVD advection scheme with zts  ln_traadv_tvd_zts = ', ln_traadv_tvd_zts
       ENDIF
 
@@ -200,4 +203,5 @@
       IF( ln_traadv_ubs    )   ioptio = ioptio + 1
       IF( ln_traadv_qck    )   ioptio = ioptio + 1
+      IF( ln_traadv_tvd_zts)   ioptio = ioptio + 1
       IF( lk_esopa         )   ioptio =          1
 
@@ -214,4 +218,5 @@
       IF( ln_traadv_ubs    )   nadv =  5
       IF( ln_traadv_qck    )   nadv =  6
+      IF( ln_traadv_tvd_zts)   nadv =  7
       IF( lk_esopa         )   nadv = -1
 
@@ -224,4 +229,5 @@
          IF( nadv ==  5 )   WRITE(numout,*) '         UBS       scheme is used'
          IF( nadv ==  6 )   WRITE(numout,*) '         QUICKEST  scheme is used'
+         IF( nadv ==  7 )   WRITE(numout,*) '         TVD ZTS   scheme is used'
          IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection scheme'
       ENDIF
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 4934)
@@ -37,5 +37,6 @@
    PRIVATE
 
-   PUBLIC   tra_adv_tvd    ! routine called by step.F90
+   PUBLIC   tra_adv_tvd        ! routine called by traadv.F90
+   PUBLIC   tra_adv_tvd_zts    ! routine called by traadv.F90
 
    LOGICAL ::   l_trd   ! flag to compute trends
@@ -262,4 +263,279 @@
    END SUBROUTINE tra_adv_tvd
 
+   SUBROUTINE tra_adv_tvd_zts ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
+      &                                       ptb, ptn, pta, kjpt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_adv_tvd_zts  ***
+      !! 
+      !! **  Purpose :   Compute the now trend due to total advection of 
+      !!       tracers and add it to the general trend of tracer equations
+      !!
+      !! **  Method  :   TVD ZTS scheme, i.e. 2nd order centered scheme with
+      !!       corrected flux (monotonic correction). This version use sub-
+      !!       timestepping for the vertical advection which increases stability
+      !!       when vertical metrics are small.
+      !!       note: - this advection scheme needs a leap-frog time scheme
+      !!
+      !! ** Action : - update (pta) with the now advective tracer trends
+      !!             - save the trends 
+      !!----------------------------------------------------------------------
+      USE oce     , ONLY:   zwx => ua        , zwy => va          ! (ua,va) used as workspace
+      !
+      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
+      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
+      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
+      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend 
+      !
+      REAL(wp), DIMENSION( jpk )                           ::   zts             ! length of sub-timestep for vertical advection
+      REAL(wp), DIMENSION( jpk )                           ::   zr_p2dt         ! reciprocal of tracer timestep
+      INTEGER  ::   ji, jj, jk, jl, jn       ! dummy loop indices  
+      INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection
+      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
+      INTEGER  ::   jtaken          ! toggle for collecting appropriate fluxes from sub timesteps
+      REAL(wp) ::   z_rzts          ! Fractional length of Euler forward sub-timestep for vertical advection
+      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar
+      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      -
+      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zwx_sav , zwy_sav
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwz, zhdiv, zwz_sav, zwzts
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrs
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tra_adv_tvd_zts')
+      !
+      CALL wrk_alloc( jpi, jpj, zwx_sav, zwy_sav )
+      CALL wrk_alloc( jpi, jpj, jpk, zwi, zwz , zhdiv, zwz_sav, zwzts )
+      CALL wrk_alloc( jpi, jpj, jpk, 3, ztrs )
+      !
+      IF( kt == kit000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_adv_tvd_zts : TVD ZTS advection scheme on ', cdtype
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+      !
+      l_trd = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
+      !
+      IF( l_trd )  THEN
+         CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
+         ztrdx(:,:,:) = 0._wp  ;    ztrdy(:,:,:) = 0._wp  ;   ztrdz(:,:,:) = 0._wp
+      ENDIF
+      !
+      zwi(:,:,:) = 0._wp
+      z_rzts = 1._wp / REAL( jnzts, wp )
+      zr_p2dt(:) = 1._wp / p2dt(:)
+      !
+      !                                                          ! ===========
+      DO jn = 1, kjpt                                            ! tracer loop
+         !                                                       ! ===========
+         ! 1. Bottom value : flux set to zero
+         ! ----------------------------------
+         zwx(:,:,jpk) = 0._wp   ;    zwz(:,:,jpk) = 0._wp
+         zwy(:,:,jpk) = 0._wp   ;    zwi(:,:,jpk) = 0._wp
+
+         ! 2. upstream advection with initial mass fluxes & intermediate update
+         ! --------------------------------------------------------------------
+         ! upstream tracer flux in the i and j direction
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  ! upstream scheme
+                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
+                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
+                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
+                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
+                  zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) )
+                  zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) )
+               END DO
+            END DO
+         END DO
+
+         ! upstream tracer flux in the k direction
+         ! Surface value
+         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0._wp                        ! volume variable
+         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface 
+         ENDIF
+         ! Interior value
+         DO jk = 2, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
+                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
+                  zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )
+               END DO
+            END DO
+         END DO
+
+         ! total advective trend
+         DO jk = 1, jpkm1
+            z2dtt = p2dt(jk)
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  ! total intermediate advective trends
+                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
+                  ! update and guess with monotonic sheme
+                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra
+                  zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !                             ! Lateral boundary conditions on zwi  (unchanged sign)
+         CALL lbc_lnk( zwi, 'T', 1. )  
+
+         !                                 ! trend diagnostics (contribution of upstream fluxes)
+         IF( l_trd )  THEN 
+            ! store intermediate advective trends
+            ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:)
+         END IF
+         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN  
+           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
+           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
+         ENDIF
+
+         ! 3. antidiffusive flux : high order minus low order
+         ! --------------------------------------------------
+         ! antidiffusive flux on i and j
+
+
+         DO jk = 1, jpkm1
+
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwx_sav(ji,jj) = zwx(ji,jj,jk)
+                  zwy_sav(ji,jj) = zwy(ji,jj,jk)
+
+                  zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) )
+                  zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) )
+               END DO
+            END DO
+
+            DO jj = 2, jpjm1         ! partial horizontal divergence
+               DO ji = fs_2, fs_jpim1
+                  zhdiv(ji,jj,jk) = (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
+                     &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
+               END DO
+            END DO
+
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwx(ji,jj,jk) = zwx(ji,jj,jk)  - zwx_sav(ji,jj)
+                  zwy(ji,jj,jk) = zwy(ji,jj,jk)  - zwy_sav(ji,jj)
+               END DO
+            END DO
+         END DO
+      
+         ! antidiffusive flux on k
+         zwz(:,:,1) = 0._wp        ! Surface value
+         zwz_sav(:,:,:) = zwz(:,:,:)
+         !
+         ztrs(:,:,:,1) = ptb(:,:,:,jn)
+         zwzts(:,:,:) = 0._wp
+
+         DO jl = 1, jnzts                   ! Start of sub timestepping loop
+
+            IF( jl == 1 ) THEN              ! Euler forward to kick things off
+              jtb = 1   ;   jtn = 1   ;   jta = 2
+              zts(:) = p2dt(:) * z_rzts
+              jtaken = MOD( jnzts + 1 , 2)  ! Toggle to collect every second flux
+                                            ! starting at jl =1 if jnzts is odd; 
+                                            ! starting at jl =2 otherwise
+            ELSEIF( jl == 2 ) THEN          ! First leapfrog step
+              jtb = 1   ;   jtn = 2   ;   jta = 3
+              zts(:) = 2._wp * p2dt(:) * z_rzts
+            ELSE                            ! Shuffle pointers for subsequent leapfrog steps
+              jtb = MOD(jtb,3) + 1
+              jtn = MOD(jtn,3) + 1
+              jta = MOD(jta,3) + 1
+            ENDIF
+            DO jk = 2, jpkm1          ! Interior value
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) )
+                     IF( jtaken == 0 ) zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk)*zts(jk)           ! Accumulate time-weighted vertcal flux
+                  END DO
+               END DO
+            END DO
+
+            jtaken = MOD( jtaken + 1 , 2 )
+
+            DO jk = 2, jpkm1          ! Interior value
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                     ! total advective trends
+                     ztra = - zbtr * (  zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
+                     ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) + zts(jk) * ztra
+                  END DO
+               END DO
+            END DO
+
+         END DO
+
+         DO jk = 2, jpkm1          ! Anti-diffusive vertical flux using average flux from the sub-timestepping
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  zwz(ji,jj,jk) = zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
+         CALL lbc_lnk( zwz, 'W',  1. )
+
+         ! 4. monotonicity algorithm
+         ! -------------------------
+         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
+
+
+         ! 5. final trend with corrected fluxes
+         ! ------------------------------------
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.  
+                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  ! total advective trends
+                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
+                  ! add them to the general tracer trends
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
+               END DO
+            END DO
+         END DO
+
+         !                                 ! trend diagnostics (contribution of upstream fluxes)
+         IF( l_trd )  THEN 
+            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
+            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
+            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
+            
+            CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, ztrdx, pun, ptn(:,:,:,jn) )   
+            CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, ztrdy, pvn, ptn(:,:,:,jn) )  
+            CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
+         END IF
+         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN  
+           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:)
+           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:)
+         ENDIF
+         !
+      END DO
+      !
+                   CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz, zhdiv, zwz_sav, zwzts )
+                   CALL wrk_dealloc( jpi, jpj, jpk, 3, ztrs )
+                   CALL wrk_dealloc( jpi, jpj, zwx_sav, zwy_sav )
+      IF( l_trd )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_tvd_zts')
+      !
+   END SUBROUTINE tra_adv_tvd_zts
 
    SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fadd_keys.sh
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fadd_keys.sh	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fadd_keys.sh	(revision 4934)
@@ -33,6 +33,6 @@
 #
 #
-# Script to add a set off key when compiling a configuration.
-# The list off key to be added has to be enclosed with " ". 
+# Script to add a set of key when compiling a configuration.
+# The list of key to be added has to be enclosed with " ". 
 # A 'sed' is performed to modify the CONFIG_NAME/cpp.fcm file to    
 # add the new key(s). 
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fcheck_config.sh
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fcheck_config.sh	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fcheck_config.sh	(revision 4934)
@@ -39,5 +39,5 @@
 # - Nothing set, use the previous in use 
 #
-# We use TOOLS/CONFIG_DIR/cfg.txt to check if the onfiguration exists.
+# We use TOOLS/CONFIG_DIR/cfg.txt to check if the configuration exists.
 #
 # EXAMPLES
@@ -76,5 +76,5 @@
    echo "Use makenemo -n MYCONFIG"
    echo "or  makenemo -h for help"
-   echo "Using defaut configuration : ${NEW_CONF}"
+   echo "Using default configuration : ${NEW_CONF}"
 fi
 if [ "$1" == cfg.txt ]; then
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fclean_var.sh
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fclean_var.sh	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fclean_var.sh	(revision 4934)
@@ -11,5 +11,5 @@
 #
 # ----------------------------
-# Clean environement variables
+# Clean environment variables
 # ----------------------------
 #
@@ -26,5 +26,5 @@
 #
 #
-# Clean environement variables
+# Clean environment variables
 #
 # EXAMPLES
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fcopy_dir.sh
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fcopy_dir.sh	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fcopy_dir.sh	(revision 4934)
@@ -26,5 +26,5 @@
 #
 #
-# When a refenrence configuration is set, 
+# When a reference configuration is set, 
 # Copy NEMO sub-directories needed (OPA_SRC, TOP_SRC ...)
 #
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fcopy_extdir.sh
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fcopy_extdir.sh	(revision 4934)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fcopy_extdir.sh	(revision 4934)
@@ -0,0 +1,39 @@
+#!/bin/bash
+#set -x
+set -o posix
+#set -u
+#set -e
+#+
+# ============
+# Fcopy_extdir.sh
+# ============
+# --------------------------
+# Copy a reference directory
+# --------------------------
+# SYNOPSIS
+# ========
+# ::
+#  $ Fcopy_extdir.sh
+# DESCRIPTION
+# ===========
+# When an unsupported configuration is requested, 
+# Prepare sources for the NEMO sub-directories needed (OPA_SRC, TOP_SRC ...)
+# EXAMPLES
+# ========
+# ::
+#  $ ./Fcopy_extdir.sh ORCA2_LIM
+# TODO
+# ====
+# option debug
+# EVOLUTIONS
+# ==========
+# $Id: Fcopy_extdir.sh 3294 2012-01-28 16:44:18Z rblod $
+#   * creation
+#-
+grep "$1 " ${CONFIG_DIR}/uspcfg.txt > ${CONFIG_DIR}/cfg.tmp
+#
+LOCAL_REF=$(cat cfg.tmp | awk 'BEGIN {FS = "#" }{print $2}')
+TAB=$(cat cfg.tmp | awk 'BEGIN {FS = "#" }{print $3}')
+REMOTE_CTL=$(cat cfg.tmp | awk 'BEGIN {FS = "#" }{print $4}')
+#
+\rm ${CONFIG_DIR}/cfg.tmp
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fdel_keys.sh
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fdel_keys.sh	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fdel_keys.sh	(revision 4934)
@@ -34,5 +34,5 @@
 #
 # Add cpp keys when compiling a configuration, key list has to be enclosed with " ".
-# We perform a 'sed' on the CONFIG_NAME/CPP.fcm file, contianing the list of keys. 
+# We perform a 'sed' on the CONFIG_NAME/CPP.fcm file, containing the list of keys. 
 #
 # EXAMPLES
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Ffetch_extdir.sh
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Ffetch_extdir.sh	(revision 4934)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Ffetch_extdir.sh	(revision 4934)
@@ -0,0 +1,53 @@
+#!/bin/bash
+#set -x
+set -o posix
+#set -u
+#set -e
+#+
+# ===============
+# Ffetch_extdir.sh
+# ===============
+# ---------------
+# Make the config 
+# ---------------
+# SYNOPSIS
+# ========
+# ::
+#  $ Ffetch_extdir.sh
+# DESCRIPTION
+# ===========
+# - Make the config directory 
+# - Create repositories needed :
+#  
+#  - EXP00 for namelist
+#  - MY_SRC for user sources
+#  - BLD for compilation 
+# EXAMPLES
+# ========
+# ::
+#  $ ./Ffetch_extdir.sh CONFIG_NAME REMOTE_CTL 
+# TODO
+# ====
+# option debug
+# EVOLUTIONS
+# ==========
+# $Id: Ffetch_extdir.sh 3715 2012-11-28 16:06:02Z acc $
+#   * creation
+#-
+basedir=$(pwd)
+cd ${1}
+wget ${2} -O remote_file.list
+#
+if [ -f remote_file.list ] ; then
+ cat remote_file.list | grep -v '^#' |
+  while
+   read remfile locfile
+   do
+    if [ $remfile == 'create_directory' ] ;then
+      mkdir $locfile
+    else
+      wget $remfile -O $locfile
+    fi
+   done
+fi
+cd $basedir
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fmake_WORK.sh
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fmake_WORK.sh	(revision 4932)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fmake_WORK.sh	(revision 4934)
@@ -35,5 +35,5 @@
 # Make the WORK directory:
 #
-# - Create lin in NEW_CONF/WORK
+# - Create line in NEW_CONF/WORK
 # - Use specified sub-directories previously
 # - OPA has to be done first !!!
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fmake_extconfig.sh
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fmake_extconfig.sh	(revision 4934)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/COMPILE/Fmake_extconfig.sh	(revision 4934)
@@ -0,0 +1,48 @@
+#!/bin/bash
+#set -x
+set -o posix
+#set -u
+#set -e
+#+
+# ===============
+# Fmake_extconfig.sh
+# ===============
+# ---------------
+# Make the directories for an externally supported configuration
+# and base the initial versions on the nearest equivalent from the
+# reference set (as named in uspcfg.txt)
+# ---------------
+# SYNOPSIS
+# ========
+# ::
+#  $ Fmake_extconfig.sh
+# DESCRIPTION
+# ===========
+# - Make the config directory 
+# - Create repositories needed :
+#  - EXP00 for namelist
+#  - MY_SRC for user sources
+#  - BLD for compilation 
+# EXAMPLES
+# ========
+# ::
+#  $ ./Fmake_extconfig.sh CONFIG_NAME REF_CONFIG_NAME 
+# TODO
+# ====
+# option debug
+# EVOLUTIONS
+# ==========
+# $Id: Fmake_extconfig.sh 3715 2012-11-28 16:06:02Z acc $
+#   * creation
+#-
+\mkdir  ${1}
+\mkdir  ${1}/EXP00
+\mkdir  ${1}/MY_SRC
+\cp -R  ${2}/cpp_${2}.fcm ${1}/cpp_${1}.fcm
+\cp -R  ${2}/EXP00/*namelist* ${1}/EXP00/.
+[ -f ${2}/EXP00/AGRIF_FixedGrids.in ] &&  \cp -R  ${2}/EXP00/AGRIF_FixedGrids.in ${1}/EXP00/.
+[ -f ${2}/EXP00/iodef.xml ]           &&  \cp -R  ${2}/EXP00/iodef.xml ${1}/EXP00/.
+[ -f ${2}/EXP00/field_def.xml ]       &&  \cp -R  ${2}/EXP00/field_def.xml ${1}/EXP00/.
+[ -f ${2}/EXP00/domain_def.xml ]      &&  \cp -R  ${2}/EXP00/domain_def.xml ${1}/EXP00/.
+[ -f ${2}/EXP00/xmlio_server.def ]    &&  \cp -R  ${2}/EXP00/xmlio_server.def ${1}/EXP00/.
+[ -d    ${2}/MY_SRC ]                 &&  \cp  ${2}/MY_SRC/* ${1}/MY_SRC/. 2> /dev/null
Index: branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/MISCELLANEOUS/icb_pp.py
===================================================================
--- branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/MISCELLANEOUS/icb_pp.py	(revision 4934)
+++ branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/TOOLS/MISCELLANEOUS/icb_pp.py	(revision 4934)
@@ -0,0 +1,159 @@
+from netCDF4 import Dataset
+from argparse import ArgumentParser
+import numpy as np
+import sys
+
+#
+# Basic iceberg trajectory post-processing python script.
+# This script collates iceberg trajectories from the distributed datasets written
+# out by each processing region and rearranges the ragged arrays into contiguous
+# streams for each unique iceberg. The output arrays are 2D (ntraj, ntimes) arrays.
+# Note that some icebergs may only exist for a subset of the possible times. In these
+# cases the missing instances are filled with invalid (NaN) values.
+#
+
+parser = ArgumentParser(description='produce collated trajectory file from distributed output\
+                                     files, e.g. \n python ./icb_pp.py \
+                                     -t  trajectory_icebergs_004248_ -n 296 -o trajsout.nc' )
+
+parser.add_argument('-t',dest='froot',help='fileroot_of_distrbuted_data; root name of \
+                     distributed trajectory output (usually completed with XXXX.nc, where \
+                     XXXX is the 4 digit processor number)', 
+                     default='trajectory_icebergs_004248_')
+
+parser.add_argument('-n',dest='fnum',help='number of distributed files to process', 
+                     type=int, default=None)
+
+parser.add_argument('-o',dest='fout',help='collated_output_file; file name to receive the \
+                     collated trajectory data', default='trajsout.nc')
+
+args = parser.parse_args()
+
+default_used = 0
+if args.froot is None:
+    pathstart = 'trajectory_icebergs_004248_'
+    default_used = 1
+else:
+    pathstart = args.froot
+
+if args.fnum is None:
+    procnum = 0
+    default_used = 1
+else:
+    procnum = args.fnum
+
+if args.fout is None:
+    pathout = 'trajsout.nc'
+    default_used = 1
+else:
+    pathout = args.fout
+
+if default_used == 1:
+   print('At least one default value will be used; command executing is:')
+   print('icb_pp.py -t ',pathstart,' -n ',procnum,' -o ',pathout)
+
+if procnum < 1:
+   print('Need some files to collate! procnum = ',procnum)
+   sys.exit()
+
+icu = []
+times = []
+#
+# Loop through all distributed datasets to obtain the complete list
+# of iceberg identification numbers and timesteps
+#
+for n in range(procnum):
+ nn = '%4.4d' % n
+ fw = Dataset(pathstart+nn+'.nc')
+ if len(fw.dimensions['n']) > 0:
+   print pathstart+nn+'.nc'
+   ic = fw.variables['iceberg_number'][:,0]
+   ts = fw.variables['timestep'][:]
+   icv = np.unique(ic)
+   ts = np.unique(ts)
+   print('Min Max ts: ',ts.min(), ts.max())
+   print('Number unique icebergs= ',icv.shape[0])
+   icu.append(icv)
+   times.append(ts)
+ fw.close()
+#
+# Now flatten the lists and reduce to the unique spanning set
+#
+icu = np.concatenate(icu)
+icu = np.unique(icu)
+times = np.concatenate(times)
+times = np.unique(times)
+ntraj = icu.shape[0]
+print(ntraj, ' unique icebergs found across all datasets')
+print('Icebergs ids range from: ',icu.min(), 'to: ',icu.max())
+print('times range from:        ',times.min(), 'to: ', times.max())
+#
+# Declare 2-D arrays to receive the data from all files
+#
+nt = times.shape[0]
+lons = np.zeros((ntraj, nt))
+lats = np.zeros((ntraj, nt))
+tims = np.zeros((ntraj, nt))
+xis  = np.zeros((ntraj, nt))
+yjs  = np.zeros((ntraj, nt))
+#
+# initially fill with invalid data
+#
+lons.fill(np.nan)
+lats.fill(np.nan)
+xis.fill(np.nan)
+yjs.fill(np.nan)
+tims.fill(np.nan)
+#
+# loop through distributed datasets again, this time
+# checking indices against icu and times lists and
+# inserting data into the correct locations in the 
+# 2-D collated sets.
+#
+for n in range(procnum):
+ nn = '%4.4d' % n
+ fw = Dataset(pathstart+nn+'.nc')
+# Note many distributed datafiles will contain no iceberg data
+# so skip quickly over these
+ m  = len(fw.dimensions['n'])
+ if m > 0:
+   inx = np.zeros(m, dtype=int)
+   tsx = np.zeros(m, dtype=int)
+   print pathstart+nn+'.nc'
+   ic = fw.variables['iceberg_number'][:,0]
+   ts = fw.variables['timestep'][:]
+   lns = fw.variables['lon'][:]
+   lts = fw.variables['lat'][:]
+   xxs = fw.variables['xi'][:]
+   yys = fw.variables['yj'][:]
+   for k in range(m):
+     inxx   = np.where(icu == ic[k])
+     inx[k] = inxx[0]
+   for k in range(m):
+     inxx   = np.where(times == ts[k])
+     tsx[k] = inxx[0]
+   lons[inx[:],tsx[:]] = lns[:]
+   lats[inx[:],tsx[:]] = lts[:]
+   tims[inx[:],tsx[:]] = ts[:]
+   xis[inx[:],tsx[:]] = xxs[:]
+   yjs[inx[:],tsx[:]] = yys[:]
+ fw.close()
+
+# Finally create the output file and write out the collated sets
+#
+fo = Dataset(pathout, 'w', format='NETCDF4')
+ntrj = fo.createDimension('ntraj', ntraj)
+nti  = fo.createDimension('ntime', None)
+olon = fo.createVariable('lon', 'f4',('ntraj','ntime'))
+olat = fo.createVariable('lat', 'f4',('ntraj','ntime'))
+otim = fo.createVariable('ttim', 'f4',('ntraj','ntime'))
+oxis = fo.createVariable('xis', 'f4',('ntraj','ntime'))
+oyjs = fo.createVariable('yjs', 'f4',('ntraj','ntime'))
+icbn = fo.createVariable('icbn', 'f4',('ntraj'))
+olon[:,:] = lons
+olat[:,:] = lats
+otim[:,:] = tims
+oxis[:,:] = xis
+oyjs[:,:] = yjs
+icbn[:] = icu
+fo.close()
