Index: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 11070)
+++ branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 11080)
@@ -759,4 +759,9 @@
       ENDDO
 
+      DO jfld = 1, nb_bdy_fld_sum
+         bf(jfld)%igrd = igrid(jfld)
+         bf(jfld)%ibdy = ibdy(jfld)
+      ENDDO
+
       ! Initialise local boundary data arrays
       ! nn_xxx_dta=0 : allocate space - will be filled from initial conditions later
Index: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 11070)
+++ branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 11080)
@@ -31,4 +31,5 @@
    USE dynadv         ! dynamics: vector invariant versus flux form
    USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
+   USE dynspg
    USE domvvl         ! variable volume
    USE bdy_oce   , ONLY: ln_bdy
@@ -179,4 +180,8 @@
             vn(:,:,jk) = va(:,:,jk)
          END DO
+         ! limit velocities
+         IF (ln_ulimit) THEN
+            call dyn_limit_velocity (kt)
+         ENDIF
          IF(.NOT.ln_linssh ) THEN
             DO jk = 1, jpkm1
@@ -203,4 +208,8 @@
                END DO
             END DO
+            ! limit velocities
+            IF (ln_ulimit) THEN
+               call dyn_limit_velocity (kt)
+            ENDIF
             !                             ! ================!
          ELSE                             ! Variable volume !
@@ -250,4 +259,8 @@
                   END DO
                END DO
+               ! limit velocities
+               IF (ln_ulimit) THEN
+                  call dyn_limit_velocity (kt)
+               ENDIF
                !
             ELSE                          ! Asselin filter applied on thickness weighted velocity
@@ -277,4 +290,8 @@
                   END DO
                END DO
+               ! limit velocities
+               IF (ln_ulimit) THEN
+                  call dyn_limit_velocity (kt)
+               ENDIF
                e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
                e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
@@ -353,4 +370,56 @@
    END SUBROUTINE dyn_nxt
 
+   SUBROUTINE dyn_limit_velocity (kt)
+      ! limits maxming vlaues of un and vn by volume courant number
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) :: zzu,zplim,zmlim,isp,ism,zcn,ze3e1,zzcn,zcnn,idivp,idivm
+
+      ! limit fluxes
+      zcn =cn_ulimit !0.9 ! maximum velocity inverse courant number
+      zcnn = cnn_ulimit !0.54 ! how much to reduce cn by in divergen flow
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               ! U direction
+               zzu = un(ji,jj,jk)
+               ze3e1 = e3u_n(ji  ,jj,jk) * e2u(ji  ,jj)
+               ! ips is 1 if flow is positive othersize zero
+               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp )
+               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )
+               ! idev is 1 if divergent flow otherwise zero
+               idivp = -isp * 0.5 * (sign(1.0_wp, un(ji-1,jj,jk)) - 1.0_wp )
+               idivm =  ism * 0.5 * (sign(1.0_wp, un(ji+1,jj,jk)) + 1.0_wp )
+               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp
+               zzcn = zzcn * zcn
+               zplim =  zzcn * (e3t_n(ji  ,jj,jk) * e1t(ji  ,jj) * e2t(ji  ,jj)) / &
+                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk)
+               zmlim = -zzcn * (e3t_n(ji+1,jj,jk) * e1t(ji+1,jj) * e2t(ji+1,jj)) / &
+                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk)
+               ! limit currents
+               un(ji,jj,jk) = min ( zzu,zplim) * isp + max (zzu,zmlim) *ism
+               ! V  direction
+               zzu = vn(ji,jj,jk)
+               ze3e1 = e3v_n(ji  ,jj,jk) * e1v(ji  ,jj)
+               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp )
+               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )
+               ! idev is 1 if divergent flow otherwise zero
+               idivp = -isp * 0.5 * (sign(1.0_wp, vn(ji,jj-1,jk)) - 1.0_wp )
+               idivm =  ism * 0.5 * (sign(1.0_wp, vn(ji,jj+1,jk)) + 1.0_wp )
+               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp
+               zzcn = zzcn * zcn
+               zplim =  zzcn * (e3t_n(ji,jj  ,jk) * e1t(ji,jj  ) * e2t(ji,jj  )) / &
+                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk)
+               zmlim = -zzcn * (e3t_n(ji,jj+1,jk) * e1t(ji,jj+1) * e2t(ji,jj+1)) / &
+                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk)
+               vn(ji,jj,jk) = min ( zzu,zplim) * isp + max (zzu,zmlim) *ism
+            ENDDO
+         ENDDO
+      ENDDO
+
+    END SUBROUTINE dyn_limit_velocity
+
    !!=========================================================================
 END MODULE dynnxt
Index: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 11070)
+++ branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 11080)
@@ -39,4 +39,7 @@
    INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from lk_dynspg_... 
 
+   LOGICAL, PUBLIC ::  ln_ulimit
+   REAL(wp), PUBLIC :: cn_ulimit,cnn_ulimit
+
    !                       ! Parameter to control the surface pressure gradient scheme
    INTEGER, PARAMETER ::   np_TS  = 1   ! split-explicit time stepping (Time-Splitting)
@@ -183,5 +186,6 @@
       NAMELIST/namdyn_spg/ ln_dynspg_exp       , ln_dynspg_ts,   &
       &                    ln_bt_fw, ln_bt_av  , ln_bt_auto  ,   &
-      &                    nn_baro , rn_bt_cmax, nn_bt_flt
+      &                    nn_baro , rn_bt_cmax, nn_bt_flt   ,   &
+                           ln_ulimit, cn_ulimit, cnn_ulimit
       !!----------------------------------------------------------------------
       !
@@ -203,4 +207,7 @@
          WRITE(numout,*) '     Explicit free surface                  ln_dynspg_exp = ', ln_dynspg_exp
          WRITE(numout,*) '     Free surface with time splitting       ln_dynspg_ts  = ', ln_dynspg_ts
+         WRITE(numout,*) '     Limit velocities                       ln_ulimit     = ', ln_ulimit
+         WRITE(numout,*) '     Limit velocities      max inverse Courant number     = ', cn_ulimit
+         WRITE(numout,*) '     Limit velocities   multiplier for divergant flow     = ', cnn_ulimit
       ENDIF
       !                          ! Control of surface pressure gradient scheme options
Index: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
===================================================================
--- branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 11070)
+++ branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 11080)
@@ -810,4 +810,5 @@
                CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
                CALL iom_get ( num, jpdom_unknown, 'e3v',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
+            CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in fld_map' )
          END SELECT
 
@@ -909,7 +910,7 @@
                      WRITE(ibstr,"(I10.10)") map%ptr(ib) 
                      CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
-                     IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_n(zij,zjj,:), mask=umask(zij,zjj,:)==1),  sum(umask(zij,zjj,:)), &
-                       &                hu_n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1  , &
-                        &                dta_read(map%ptr(ib),1,:)
+                     !IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_n(zij,zjj,:), mask=umask(zij,zjj,:)==1),  sum(umask(zij,zjj,:)), &
+                     !  &                hu_n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1  , &
+                     !   &                dta_read(map%ptr(ib),1,:)
                   ENDIF
                CASE(3)
Index: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_istate.F90
===================================================================
--- branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_istate.F90	(revision 11070)
+++ branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_istate.F90	(revision 11080)
@@ -62,18 +62,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               pts(ji,jj,jk,jp_tem) =  (  (  16. - 12. * TANH( (pdept(ji,jj,jk) - 400) / 700 ) )   &
-                    &           * (-TANH( (500. - pdept(ji,jj,jk)) / 150. ) + 1.) / 2.             &
-                    &           + ( 15. * ( 1. - TANH( (pdept(ji,jj,jk)-50.) / 1500.) )            &
-                    &           - 1.4 * TANH((pdept(ji,jj,jk)-100.) / 100.)                        &
-                    &           + 7.  * (1500. - pdept(ji,jj,jk) ) / 1500.)                        &
-                    &           * (-TANH( (pdept(ji,jj,jk) - 500.) / 150.) + 1.) / 2.  ) * ptmask(ji,jj,jk)
-
-               pts(ji,jj,jk,jp_sal) =  (  (  36.25 - 1.13 * TANH( (pdept(ji,jj,jk) - 305) / 460 ) )  &
-                    &         * (-TANH((500. - pdept(ji,jj,jk)) / 150.) + 1.) / 2                  &
-                    &         + ( 35.55 + 1.25 * (5000. - pdept(ji,jj,jk)) / 5000.                 &
-                    &         - 1.62 * TANH( (pdept(ji,jj,jk) - 60.  ) / 650. )                    &
-                    &         + 0.2  * TANH( (pdept(ji,jj,jk) - 35.  ) / 100. )                    &
-                    &         + 0.2  * TANH( (pdept(ji,jj,jk) - 1000.) / 5000.) )                  &
-                    &         * (-TANH( (pdept(ji,jj,jk) - 500.) / 150.) + 1.) / 2  ) * ptmask(ji,jj,jk)
+               pts(ji,jj,jk,jp_tem) = 20._wp * ptmask(ji,jj,jk)
+               pts(ji,jj,jk,jp_sal) = 36.25_wp * ptmask(ji,jj,jk)
             END DO
          END DO
Index: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_sbc.F90
===================================================================
--- branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_sbc.F90	(revision 11070)
+++ branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_sbc.F90	(revision 11080)
@@ -1,7 +1,7 @@
 MODULE usrdef_sbc
    !!======================================================================
-   !!                     ***  MODULE  usrdef_sbc  ***
-   !!
-   !!                     ===  GYRE configuration  ===
+   !!                       ***  MODULE usrdef_sbc  ***
+   !! 
+   !!                  ===  WAD_TEST_CASES configuration  ===
    !!
    !! User defined :   surface forcing of a user configuration
@@ -11,5 +11,5 @@
 
    !!----------------------------------------------------------------------
-   !!   usrdef_sbc    : user defined surface bounday conditions in GYRE case
+   !!   usrdef_sbc     : user defined surface bounday conditions in WAD_TEST_CASES case
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
@@ -21,5 +21,5 @@
    USE lib_mpp         ! distribued memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE lib_fortran     !
+   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
 
    IMPLICIT NONE
@@ -41,184 +41,34 @@
    SUBROUTINE usrdef_sbc_oce( kt )
       !!---------------------------------------------------------------------
-      !!                    ***  ROUTINE usrdef_sbc  ***
+      !!                    ***  ROUTINE usr_def_sbc  ***
       !!              
-      !! ** Purpose :   provide at each time-step the GYRE surface boundary
+      !! ** Purpose :   provide at each time-step the surface boundary
       !!              condition, i.e. the momentum, heat and freshwater fluxes.
       !!
-      !! ** Method  :   analytical seasonal cycle for GYRE configuration.
+      !! ** Method  :   all 0 fields, for WAD_TEST_CASES case
       !!                CAUTION : never mask the surface stress field !
       !!
-      !! ** Action  : - set the ocean surface boundary condition, i.e.   
+      !! ** Action  : - set to ZERO all the ocean surface boundary condition, i.e.   
       !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx
       !!
-      !! Reference : Hazeleger, W., and S. Drijfhout, JPO, 30, 677-695, 2000.
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
-      !!
-      INTEGER  ::   ji, jj                 ! dummy loop indices
-      INTEGER  ::   zyear0                 ! initial year 
-      INTEGER  ::   zmonth0                ! initial month
-      INTEGER  ::   zday0                  ! initial day
-      INTEGER  ::   zday_year0             ! initial day since january 1st
-      REAL(wp) ::   ztau     , ztau_sais   ! wind intensity and of the seasonal cycle
-      REAL(wp) ::   ztime                  ! time in hour
-      REAL(wp) ::   ztimemax , ztimemin    ! 21th June, and 21th decem. if date0 = 1st january
-      REAL(wp) ::   ztimemax1, ztimemin1   ! 21th June, and 21th decem. if date0 = 1st january
-      REAL(wp) ::   ztimemax2, ztimemin2   ! 21th June, and 21th decem. if date0 = 1st january
-      REAL(wp) ::   ztaun                  ! intensity
-      REAL(wp) ::   zemp_s, zemp_n, zemp_sais, ztstar
-      REAL(wp) ::   zcos_sais1, zcos_sais2, ztrp, zconv, t_star
-      REAL(wp) ::   zsumemp, zsurf
-      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
-      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
-      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
-      REAL(wp) ::   zyydd                 ! number of days in one year
       !!---------------------------------------------------------------------
-      zyydd = REAL(nyear_len(1),wp)
-
-      ! ---------------------------- !
-      !  heat and freshwater fluxes  !
-      ! ---------------------------- !
-      !same temperature, E-P as in HAZELEGER 2000
-
-      zyear0     =   ndate0 / 10000                             ! initial year
-      zmonth0    = ( ndate0 - zyear0 * 10000 ) / 100            ! initial month
-      zday0      =   ndate0 - zyear0 * 10000 - zmonth0 * 100    ! initial day betwen 1 and 30
-      zday_year0 = ( zmonth0 - 1 ) * 30.+zday0                  ! initial day betwen 1 and 360
-
-      ! current day (in hours) since january the 1st of the current year
-      ztime = REAL( kt ) * rdt / (rmmss * rhhmm)   &       !  total incrementation (in hours)
-         &      - (nyear  - 1) * rjjhh * zyydd             !  minus years since beginning of experiment (in hours)
-
-      ztimemax1 = ((5.*30.)+21.)* 24.                      ! 21th june     at 24h in hours
-      ztimemin1 = ztimemax1 + rjjhh * zyydd / 2            ! 21th december        in hours
-      ztimemax2 = ((6.*30.)+21.)* 24.                      ! 21th july     at 24h in hours
-      ztimemin2 = ztimemax2 - rjjhh * zyydd / 2            ! 21th january         in hours
-      !                                                    ! NB: rjjhh * zyydd / 4 = one seasonal cycle in hours
-
-      ! amplitudes
-      zemp_S    = 0.7       ! intensity of COS in the South
-      zemp_N    = 0.8       ! intensity of COS in the North
-      zemp_sais = 0.1
-      zTstar    = 28.3      ! intemsity from 28.3 a -5 deg
-
-      ! 1/2 period between 21th June and 21th December and between 21th July and 21th January
-      zcos_sais1 = COS( (ztime - ztimemax1) / (ztimemin1 - ztimemax1) * rpi ) 
-      zcos_sais2 = COS( (ztime - ztimemax2) / (ztimemax2 - ztimemin2) * rpi )
-
-      ztrp= - 40.e0        ! retroaction term on heat fluxes (W/m2/K)
-      zconv = 3.16e-5      ! convertion factor: 1 m/yr => 3.16e-5 mm/s
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3
-            ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
-            ! 64.5 in summer, 42.5 in winter
-            t_star = zTstar * ( 1. + 1. / 50. * zcos_sais2 )                &
-               &                    * COS( rpi * (gphit(ji,jj) - 5.)               &
-               &                    / ( 53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) )
-            ! 23.5 deg : tropics
-            qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) )
-            qns (ji,jj) = ztrp * ( tsb(ji,jj,1,jp_tem) - t_star ) - qsr(ji,jj)
-            IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj) ) THEN    ! zero at 37.8 deg, max at 24.6 deg
-               emp  (ji,jj) =   zemp_S * zconv   &
-                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  &
-                  &         * ( 1 - zemp_sais / zemp_S * zcos_sais1)
-            ELSE
-               emp (ji,jj) =  - zemp_N * zconv   &
-                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  &
-                  &         * ( 1 - zemp_sais / zemp_N * zcos_sais1 )
-            ENDIF
-         END DO
-      END DO
-
-      zsumemp = GLOB_SUM( emp  (:,:)   ) 
-      zsurf   = GLOB_SUM( tmask(:,:,1) ) 
-      zsumemp = zsumemp / zsurf         ! Default GYRE configuration
-
-      ! freshwater (mass flux) and update of qns with heat content of emp
-      emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1)        ! freshwater flux (=0 in domain average)
-      sfx (:,:) = 0.0_wp                                   ! no salt flux
-      qns (:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp   ! evap and precip are at SST
-
-
-      ! ---------------------------- !
-      !       momentum fluxes        !
-      ! ---------------------------- !
-      ! same wind as in Wico
-      !test date0 : ndate0 = 010203
-      zyear0  =   ndate0 / 10000
-      zmonth0 = ( ndate0 - zyear0 * 10000 ) / 100
-      zday0   =   ndate0 - zyear0 * 10000 - zmonth0 * 100
-      !Calculates nday_year, day since january 1st
-      zday_year0 = (zmonth0-1)*30.+zday0
-
-      !accumulates days of previous months of this year
-      ! day (in hours) since january the 1st
-      ztime = FLOAT( kt ) * rdt / (rmmss * rhhmm)  &  ! incrementation in hour
-         &     - (nyear - 1) * rjjhh * zyydd          !  - nber of hours the precedent years
-      ztimemax = ((5.*30.)+21.)* 24.               ! 21th june     in hours
-      ztimemin = ztimemax + rjjhh * zyydd / 2      ! 21th december in hours
-      !                                            ! NB: rjjhh * zyydd / 4 = 1 seasonal cycle in hours
-
-      ! mean intensity at 0.105 ; srqt(2) because projected with 45deg angle
-      ztau = 0.105 / SQRT( 2. )
-      ! seasonal oscillation intensity
-      ztau_sais = 0.015
-      ztaun = ztau - ztau_sais * COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi )
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-           ! domain from 15deg to 50deg and 1/2 period along 14deg
-           ! so 5/4 of half period with seasonal cycle
-           utau(ji,jj) = - ztaun * SIN( rpi * (gphiu(ji,jj) - 15.) / (29.-15.) )
-           vtau(ji,jj) =   ztaun * SIN( rpi * (gphiv(ji,jj) - 15.) / (29.-15.) )
-         END DO
-      END DO
-
-      ! module of wind stress and wind speed at T-point
-      zcoef = 1. / ( zrhoa * zcdrag ) 
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vect. opt.
-            ztx = utau(ji-1,jj  ) + utau(ji,jj) 
-            zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
-            zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
-            taum(ji,jj) = zmod
-            wndm(ji,jj) = SQRT( zmod * zcoef )
-         END DO
-      END DO
-      CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
-
-      ! ---------------------------------- !
-      !  control print at first time-step  !
-      ! ---------------------------------- !
-      IF( kt == nit000 .AND. lwp ) THEN 
-         WRITE(numout,*)
-         WRITE(numout,*)'usrdef_sbc_oce : analytical surface fluxes for GYRE configuration'               
-         WRITE(numout,*)'~~~~~~~~~~~ ' 
-         WRITE(numout,*)'           nyear      = ', nyear
-         WRITE(numout,*)'           nmonth     = ', nmonth
-         WRITE(numout,*)'           nday       = ', nday
-         WRITE(numout,*)'           nday_year  = ', nday_year
-         WRITE(numout,*)'           ztime      = ', ztime
-         WRITE(numout,*)'           ztimemax   = ', ztimemax
-         WRITE(numout,*)'           ztimemin   = ', ztimemin
-         WRITE(numout,*)'           ztimemax1  = ', ztimemax1
-         WRITE(numout,*)'           ztimemin1  = ', ztimemin1
-         WRITE(numout,*)'           ztimemax2  = ', ztimemax2
-         WRITE(numout,*)'           ztimemin2  = ', ztimemin2
-         WRITE(numout,*)'           zyear0     = ', zyear0
-         WRITE(numout,*)'           zmonth0    = ', zmonth0
-         WRITE(numout,*)'           zday0      = ', zday0
-         WRITE(numout,*)'           zday_year0 = ', zday_year0
-         WRITE(numout,*)'           zyydd      = ', zyydd
-         WRITE(numout,*)'           zemp_S     = ', zemp_S
-         WRITE(numout,*)'           zemp_N     = ', zemp_N
-         WRITE(numout,*)'           zemp_sais  = ', zemp_sais
-         WRITE(numout,*)'           zTstar     = ', zTstar
-         WRITE(numout,*)'           zsumemp    = ', zsumemp
-         WRITE(numout,*)'           zsurf      = ', zsurf
-         WRITE(numout,*)'           ztrp       = ', ztrp
-         WRITE(numout,*)'           zconv      = ', zconv
-         WRITE(numout,*)'           ndastp     = ', ndastp
-         WRITE(numout,*)'           adatrj     = ', adatrj
+      !
+      IF( kt == nit000 ) THEN
+         !
+         IF(lwp) WRITE(numout,*)' usr_sbc : WAD_TEST_CASES case: NO surface forcing'
+         IF(lwp) WRITE(numout,*)' ~~~~~~~~~~~   utau = vtau = taum = wndm = qns = qsr = emp = sfx = 0'
+         !
+         utau(:,:) = 0._wp
+         vtau(:,:) = 0._wp
+         taum(:,:) = 0._wp
+         wndm(:,:) = 0._wp
+         !
+         emp (:,:) = 0._wp
+         sfx (:,:) = 0._wp
+         qns (:,:) = 0._wp
+         qsr (:,:) = 0._wp
+         !         
       ENDIF
       !
Index: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/stpctl.F90
===================================================================
--- branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/stpctl.F90	(revision 11070)
+++ branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/stpctl.F90	(revision 11080)
@@ -52,5 +52,5 @@
       INTEGER  ::   ji, jj, jk             ! dummy loop indices
       INTEGER  ::   ii, ij, ik             ! local integers
-      REAL(wp) ::   zumax, zsmin, zssh2, zsshmax    ! local scalars
+      REAL(wp) ::   velmax2, zsmin, zssh2, zsshmax    ! local scalars
       INTEGER, DIMENSION(3) ::   ilocu     ! 
       INTEGER, DIMENSION(2) ::   ilocs     ! 
@@ -70,22 +70,22 @@
       !                                              !* Test maximum of velocity (zonal only)
       !                                              !  ------------------------
-      !! zumax = MAXVAL( ABS( un(:,:,:) ) )                ! slower than the following loop on NEC SX5
-      zumax = 0.e0
+      !! velmax2 = MAXVAL( ABS( un(:,:,:) ) )                ! slower than the following loop on NEC SX5
+      velmax2 = 0.e0
       DO jk = 1, jpk
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zumax = MAX(zumax,ABS(un(ji,jj,jk)))
+               velmax2 = MAX( velmax2,un(ji,jj,jk)**2 + vn(ji,jj,jk)**2 )
           END DO 
         END DO 
       END DO        
-      IF( lk_mpp )   CALL mpp_max( zumax )                 ! max over the global domain
+      IF( lk_mpp )   CALL mpp_max( velmax2 )               ! max over the global domain
       !
-      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,' abs(U) max: ', zumax
+      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,' 3d speed max: ', velmax2
       !
-      IF( zumax > 20.e0 ) THEN
+      IF( velmax2 > 20.e0**2 ) THEN
          IF( lk_mpp ) THEN
-            CALL mpp_maxloc(ABS(un),umask,zumax,ii,ij,ik)
+            CALL mpp_maxloc( un(:,:,:)**2+vn(:,:,:)**2,umask,velmax2,ii,ij,ik)
          ELSE
-            ilocu = MAXLOC( ABS( un(:,:,:) ) )
+            ilocu = MAXLOC( un(:,:,:)**2 + vn(:,:,:)**2 )
             ii = ilocu(1) + nimpp - 1
             ij = ilocu(2) + njmpp - 1
@@ -94,7 +94,7 @@
          IF(lwp) THEN
             WRITE(numout,cform_err)
-            WRITE(numout,*) ' stpctl: the zonal velocity is larger than 20 m/s'
+            WRITE(numout,*) ' stpctl: the speed is larger than 20 m/s'
             WRITE(numout,*) ' ====== '
-            WRITE(numout,9400) kt, zumax, ii, ij, ik
+            WRITE(numout,9400) kt, velmax2, ii, ij, ik
             WRITE(numout,*)
             WRITE(numout,*) '          output of last fields in numwso'
@@ -102,5 +102,5 @@
          kindic = -3
       ENDIF
-9400  FORMAT (' kt=',i6,' max abs(U): ',1pg11.4,', i j k: ',3i5)
+9400  FORMAT (' kt=',i6,' max abs(vel)**2: ',1pg11.4,', i j k: ',3i5)
       !
       !                                              !* Test minimum of salinity
@@ -180,7 +180,7 @@
       zssh2 = glob_sum( sshn(:,:) * sshn(:,:) )
       !
-      IF(lwp) WRITE(numsol,9700) kt, zssh2, zumax, zsmin      ! ssh statistics
+      IF(lwp) WRITE(numsol,9700) kt, zssh2, velmax2, zsmin      ! ssh statistics
       !
-9700  FORMAT(' it :', i8, ' ssh2: ', d23.16, ' Umax: ',d23.16,' Smin: ',d23.16)
+9700  FORMAT(' it :', i8, ' ssh2: ', d23.16, ' velmax2: ',d23.16,' SSSmin: ',d23.16)
       !
    END SUBROUTINE stp_ctl
