Index: branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Biblio/Biblio.bib
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Biblio/Biblio.bib	(revision 2985)
+++ branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Biblio/Biblio.bib	(revision 2986)
@@ -1275,4 +1275,47 @@
   url = {http://dx.doi.org/10.1016/j.ocemod.2009.12.003},
   issn = {1463-5003},
+}
+
+@ARTICLE{HollowayOM86,
+  author = {Greg Holloway},
+  title = {A Shelf Wave/Topographic Pump Drives Mean Coastal Circulation (part I)},
+  journal = OM,
+  year = {1986},
+  volume = {68},  
+}
+
+@ARTICLE{HollowayJPO92,
+  author = {Greg Holloway},
+  title = {Representing Topographic Stress for Large-Scale Ocean Models},
+  journal = JPO,
+  year = {1992},
+  volume = {22},  
+  pages = {1033--1046},
+}
+
+@ARTICLE{HollowayJPO94,
+  author = {Michael Eby and Greg Holloway},
+  title = {Sensitivity of a Large-Scale Ocean Model to a Parameterization of Topographic Stress},
+  journal = JPO,
+  year = {1994},
+  volume = {24},  
+  pages = {2577--2587},
+}
+
+@ARTICLE{HollowayJGR09,
+  author = {Greg Holloway and Zeliang Wang},
+  title = {Representing eddy stress in an Arctic Ocean model},
+  journal = JGR,
+  year = {2009},
+  doi = {10.1029/2008JC005169},  
+}
+
+@ARTICLE{HollowayOM08,
+  author = {Mathew Maltrud and Greg Holloway},
+  title = {Implementing biharmonic neptune in a global eddying ocean model},
+  journal = OM,
+  year = {2008},
+  volume = {21},  
+  pages = {22--34},
 }
 
Index: branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Chapters/Chap_DYN.tex
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 2985)
+++ branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 2986)
@@ -1162,2 +1162,85 @@
 
 % ================================================================
+% Neptune effect 
+% ================================================================
+\section  [Neptune effect (\textit{dynnept})]
+                {Neptune effect (\mdl{dynnept})}
+\label{DYN_nept}
+
+The "Neptune effect" (thus named in \citep{HollowayOM86}) is a
+parameterisation of the potentially large effect of topographic form stress
+(caused by eddies) in driving the ocean circulation. Originally developed for
+low-resolution models, in which it was applied via a Laplacian (second-order)
+diffusion-like term in the momentum equation, it can also be applied in eddy
+permitting or resolving models, in which a more scale-selective bilaplacian
+(fourth-order) implementation is preferred. This mechanism has a
+significant effect on boundary currents (including undercurrents), and the
+upwelling of deep water near continental shelves.
+
+The theoretical basis for the method can be found in 
+\citep{HollowayJPO92}, including the explanation of why form stress is not
+necessarily a drag force, but may actually drive the flow. 
+\citep{HollowayJPO94} demonstrate the effects of the parameterisation in
+the GFDL-MOM model, at a horizontal resolution of about 1.8 degrees. 
+\citep{HollowayOM08} demonstrate the biharmonic version of the
+parameterisation in a global run of the POP model, with an average horizontal
+grid spacing of about 32km.
+
+The NEMO implementation is a simplified form of that supplied by
+Greg Holloway, the testing of which was described in \citep{HollowayJGR09}.
+The major simplification is that a time invariant Neptune velocity
+field is assumed.  This is computed only once, during start-up, and
+made available to the rest of the code via a module.  Vertical
+diffusive terms are also ignored, and the model topography itself
+is used, rather than a separate topographic dataset as in
+\citep{HollowayOM08}.  This implementation is only in the iso-level
+formulation, as is the case anyway for the bilaplacian operator.
+
+The velocity field is derived from a transport stream function given by:
+
+\begin{equation} \label{Eq_dynnept_sf}
+\psi = -fL^2H
+\end{equation}
+
+where $L$ is a latitude-dependant length scale given by:
+
+\begin{equation} \label{Eq_dynnept_ls}
+L = l_1 + (l_2 -l_1)\left ( {1 + \cos 2\phi \over 2 } \right )
+\end{equation}
+
+where $\phi$ is latitude and $l_1$ and $l_2$ are polar and equatorial length scales respectively.
+Neptune velocity components, $u^*$, $v^*$ are derived from the stremfunction as:
+
+\begin{equation} \label{Eq_dynnept_vel}
+u^* = -{1\over H} {\partial \psi \over \partial y}\ \ \  ,\ \ \ v^* = {1\over H} {\partial \psi \over \partial x}
+\end{equation}
+
+\smallskip
+%----------------------------------------------namdom----------------------------------------------------
+\namdisplay{namdyn_nept}
+%--------------------------------------------------------------------------------------------------------
+\smallskip
+
+The Neptune effect is enabled when \np{ln\_neptsimp}=true (default=false).
+\np{ln\_smooth\_neptvel} controls whether a scale-selective smoothing is applied
+to the Neptune effect flow field (default=false) (this smoothing method is as
+used by Holloway).  \np{rn\_tslse} and \np{rn\_tslsp} are the equatorial and
+polar values respectively of the length-scale parameter $L$ used in determining
+the Neptune stream function \eqref{Eq_dynnept_sf} and \eqref{Eq_dynnept_ls}.
+Values at intermediate latitudes are given by a cosine fit, mimicking the
+variation of the deformation radius with latitude.  The default values of 12km
+and 3km are those given in \citep{HollowayJPO94}, appropriate for a coarse
+resolution model. The finer resolution study of \citep{HollowayOM08} increased
+the values of L by a factor of $\sqrt 2$ to 17km and 4.2km, thus doubling the
+stream function for a given topography.
+
+The simple formulation for ($u^*$, $v^*$) can give unacceptably large velocities
+in shallow water, and \citep{HollowayOM08} add an offset to the depth in the
+denominator to control this problem. In this implementation we offer instead (at
+the suggestion of G. Madec) the option of ramping down the Neptune flow field to
+zero over a finite depth range. The switch \np{ln\_neptramp} activates this
+option (default=false), in which case velocities at depths greater than
+\np{rn\_htrmax} are unaltered, but ramp down linearly with depth to zero at a
+depth of \np{rn\_htrmin} (and shallower).
+
+% ================================================================
Index: branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Namelist/namdyn_nept
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Namelist/namdyn_nept	(revision 2986)
+++ branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Namelist/namdyn_nept	(revision 2986)
@@ -0,0 +1,15 @@
+!-----------------------------------------------------------------------
+&nam_dynnept  !   Neptune effect (simplified: lateral and vertical diffusions removed)
+!-----------------------------------------------------------------------
+   ! Suggested lengthscale values are those of Eby & Holloway (1994) for a coarse model
+   ln_neptsimp       = .false.  ! yes/no use simplified neptune
+
+   ln_smooth_neptvel = .false.  ! yes/no smooth zunep, zvnep
+   rn_tslse          =  1.2e4   ! value of lengthscale L at the equator
+   rn_tslsp          =  3.0e3   ! value of lengthscale L at the pole
+   ! Specify whether to ramp down the Neptune velocity in shallow
+   ! water, and if so the depth range controlling such ramping down
+   ln_neptramp       = .false.  ! ramp down Neptune velocity in shallow water
+   rn_htrmin         =  100.0   ! min. depth of transition range
+   rn_htrmax         =  200.0   ! max. depth of transition range
+/
Index: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/ARCH/arch-ALTIX_NAUTILUS4.fcm
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/ARCH/arch-ALTIX_NAUTILUS4.fcm	(revision 2985)
+++ branches/2011/dev_NOC_2011_MERGE/NEMOGCM/ARCH/arch-ALTIX_NAUTILUS4.fcm	(revision 2986)
@@ -22,10 +22,10 @@
 # Note use of -Bstatic because the library root directories are not accessible to the back-end compute nodes
 %NCDF_LIB            -L%HDF5_HOME/lib -L%NCDF_HOME/lib -Bstatic -lnetcdf -lhdf5_fortran -lhdf5_hl -lhdf5 -Bdynamic -lz
-%FC                  mpif90
+%FC                  ifort
 %FCFLAGS             -r8 -O3 -xT -ip -vec-report0
 %FFLAGS              -r8 -O3 -xT -ip -vec-report0
-%LD                  mpif90
+%LD                  ifort
 %FPPFLAGS            -P -C -traditional
-%LDFLAGS
+%LDFLAGS             -lmpi
 %AR                  ar 
 %ARFLAGS             -r
Index: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/GYRE/EXP00/namelist
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 2985)
+++ branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 2986)
@@ -865,2 +865,17 @@
     salfixmin = -9999      !  Minimum salinity after applying the increments
 /
+!-----------------------------------------------------------------------
+&namdyn_nept  !   Neptune effect (simplified: lateral and vertical diffusions removed)
+!-----------------------------------------------------------------------
+   ! Suggested lengthscale values are those of Eby & Holloway (1994) for a coarse model
+   ln_neptsimp       = .false.  ! yes/no use simplified neptune
+
+   ln_smooth_neptvel = .false.  ! yes/no smooth zunep, zvnep
+   rn_tslse          =  1.2e4   ! value of lengthscale L at the equator
+   rn_tslsp          =  3.0e3   ! value of lengthscale L at the pole
+   ! Specify whether to ramp down the Neptune velocity in shallow
+   ! water, and if so the depth range controlling such ramping down
+   ln_neptramp       = .false.  ! ramp down Neptune velocity in shallow water
+   rn_htrmin         =  100.0   ! min. depth of transition range
+   rn_htrmax         =  200.0   ! max. depth of transition range
+/
Index: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 2985)
+++ branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 2986)
@@ -865,2 +865,17 @@
     salfixmin = -9999      !  Minimum salinity after applying the increments
 /
+!-----------------------------------------------------------------------
+&namdyn_nept  !   Neptune effect (simplified: lateral and vertical diffusions removed)
+!-----------------------------------------------------------------------
+   ! Suggested lengthscale values are those of Eby & Holloway (1994) for a coarse model
+   ln_neptsimp       = .true.   ! yes/no use simplified neptune
+
+   ln_smooth_neptvel = .false.  ! yes/no smooth zunep, zvnep
+   rn_tslse          =  1.2e4   ! value of lengthscale L at the equator
+   rn_tslsp          =  3.0e3   ! value of lengthscale L at the pole
+   ! Specify whether to ramp down the Neptune velocity in shallow
+   ! water, and if so the depth range controlling such ramping down
+   ln_neptramp       = .true.   ! ramp down Neptune velocity in shallow water
+   rn_htrmin         =  100.0   ! min. depth of transition range
+   rn_htrmax         =  200.0   ! max. depth of transition range
+/
Index: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 2985)
+++ branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 2986)
@@ -880,2 +880,17 @@
     salfixmin = -9999      !  Minimum salinity after applying the increments
 /
+!-----------------------------------------------------------------------
+&namdyn_nept  !   Neptune effect (simplified: lateral and vertical diffusions removed)
+!-----------------------------------------------------------------------
+   ! Suggested lengthscale values are those of Eby & Holloway (1994) for a coarse model
+   ln_neptsimp       = .false.  ! yes/no use simplified neptune
+
+   ln_smooth_neptvel = .false.  ! yes/no smooth zunep, zvnep
+   rn_tslse          =  1.2e4   ! value of lengthscale L at the equator
+   rn_tslsp          =  3.0e3   ! value of lengthscale L at the pole
+   ! Specify whether to ramp down the Neptune velocity in shallow
+   ! water, and if so the depth range controlling such ramping down
+   ln_neptramp       = .false.  ! ramp down Neptune velocity in shallow water
+   rn_htrmin         =  100.0   ! min. depth of transition range
+   rn_htrmax         =  200.0   ! max. depth of transition range
+/
Index: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/POMME/EXP00/namelist
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/POMME/EXP00/namelist	(revision 2985)
+++ branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/POMME/EXP00/namelist	(revision 2986)
@@ -870,2 +870,17 @@
     salfixmin = -9999      !  Minimum salinity after applying the increments
 /
+!-----------------------------------------------------------------------
+&namdyn_nept  !   Neptune effect (simplified: lateral and vertical diffusions removed)
+!-----------------------------------------------------------------------
+   ! Suggested lengthscale values are those of Eby & Holloway (1994) for a coarse model
+   ln_neptsimp       = .false.  ! yes/no use simplified neptune
+
+   ln_smooth_neptvel = .false.  ! yes/no smooth zunep, zvnep
+   rn_tslse          =  1.2e4   ! value of lengthscale L at the equator
+   rn_tslsp          =  3.0e3   ! value of lengthscale L at the pole
+   ! Specify whether to ramp down the Neptune velocity in shallow
+   ! water, and if so the depth range controlling such ramping down
+   ln_neptramp       = .false.  ! ramp down Neptune velocity in shallow water
+   rn_htrmin         =  100.0   ! min. depth of transition range
+   rn_htrmax         =  200.0   ! max. depth of transition range
+/
Index: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90	(revision 2986)
+++ branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90	(revision 2986)
@@ -0,0 +1,621 @@
+MODULE dynnept
+   !!======================================================================
+   !!                       ***  MODULE  dynnept  ***
+   !! Ocean dynamics: Neptune effect as proposed by Greg Holloway,
+   !!                 recoded version of simplest case (u*, v* only)
+   !!======================================================================
+   !! History :  1.0  !  2007-06  (Michael Dunphy)  Modular form: - new namelist parameters
+   !!                                                             - horizontal diffusion for Neptune
+   !!                                                             - vertical diffusion for gm in momentum eqns
+   !!                                                             - option to use Neptune in Coriolis eqn
+   !!                    2011-08  (Jeff Blundell, NOCS) Simplified form for temporally invariant u*, v*
+   !!                                               Horizontal and vertical diffusivity formulations removed
+   !!                                               Dynamic allocation of storage added
+   !!                                               Option of ramping Neptune vel. down
+   !!                                               to zero added in shallow depths added
+   !!----------------------------------------------------------------------
+   !! dynnept_alloc        :
+   !! dyn_nept_init        :
+   !! dyn_nept_div_cur_init:
+   !! dyn_nept_cor         :
+   !! dyn_nept_vel         :
+   !! dyn_nept_smooth_vel  :
+   !!----------------------------------------------------------------------
+   USE oce              ! ocean dynamics and tracers
+   USE dom_oce          ! ocean space and time domain
+   USE in_out_manager   ! I/O manager
+   USE lib_mpp          ! distributed memory computing
+   USE prtctl           ! Print control
+   USE phycst
+   USE lbclnk
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC dyn_nept_init      ! routine called by nemogcm.F90
+   PUBLIC dyn_nept_cor       ! routine called by step.F90
+   !! dynnept_alloc()         is called only by dyn_nept_init, within this module
+   !! dyn_nept_div_cur_init   is called only by dyn_nept_init, within this module
+   !! dyn_nept_vel            is called only by dyn_nept_cor,  within this module
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: zunep, zvnep  ! Neptune u and v
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: zhdivnep      ! hor. div for Neptune vel.
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: zmrotnep      ! curl for Neptune vel.
+
+
+   !! * Namelist namdyn_nept variables
+   LOGICAL, PUBLIC  ::  ln_neptsimp        = .FALSE.  ! yes/no simplified neptune
+
+   LOGICAL          ::  ln_smooth_neptvel  = .FALSE.  ! yes/no smooth zunep, zvnep
+   REAL(wp)         ::  rn_tslse           =  1.2e4   ! value of lengthscale L at the equator
+   REAL(wp)         ::  rn_tslsp           =  3.0e3   ! value of lengthscale L at the pole
+!! Specify whether to ramp down the Neptune velocity in shallow
+!! water, and the depth range controlling such ramping down
+   LOGICAL          ::  ln_neptramp        = .FALSE.  ! ramp down Neptune velocity in shallow water
+   REAL(wp)         ::  rn_htrmin          =  100.0   ! min. depth of transition range
+   REAL(wp)         ::  rn_htrmax          =  200.0   ! max. depth of transition range
+
+   !! * Module variables
+
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!   OPA 9.0 , implemented by Bedford Institute of Oceanography
+   !!----------------------------------------------------------------------
+
+ CONTAINS
+
+   INTEGER FUNCTION dynnept_alloc()
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE dynnept_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( zunep(jpi,jpj) , zvnep(jpi,jpj) ,     &
+         &      zhdivnep(jpi,jpj,jpk) , zmrotnep(jpi,jpj,jpk) , STAT=dynnept_alloc )
+         !
+      IF( dynnept_alloc /= 0 )   CALL ctl_warn('dynnept_alloc: array allocate failed.')
+   END FUNCTION dynnept_alloc
+
+
+   SUBROUTINE dyn_nept_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_nept_init  ***
+      !!
+      !! ** Purpose :   Read namelist parameters, initialise arrays
+      !!                and compute the arrays zunep and zvnep
+      !!
+      !! ** Method  :   zunep =
+      !!                zvnep =
+      !!
+      !! ** History :  1.0  !   07-05  (Zeliang Wang)   Original code for zunep, zvnep
+      !!               1.1  !   07-06  (Michael Dunphy) namelist and  initialisation
+      !!               2.0  ! 2011-07  (Jeff Blundell, NOCS)
+      !!                    ! Simplified form for temporally invariant u*, v*
+      !!                    ! Horizontal and vertical diffusivity formulations removed
+      !!                    ! Includes optional tapering-off in shallow depths
+      !!----------------------------------------------------------------------
+      USE iom
+      !!
+!!    Local dynamic allocation
+      INTEGER  ::   ji, jj, jk    ! dummy loop indices
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   ht       ! temporary 2D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   htn      ! temporary 2D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   tscale   ! temporary 2D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   tsp      ! temporary 2D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   hur_n    ! temporary 2D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   hvr_n    ! temporary 2D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   hu_n     ! temporary 2D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   hv_n     ! temporary 2D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   znmask   ! temporary 3D array for nmask
+      REAL(wp) :: unemin,unemax,vnemin,vnemax   ! extrema of (u*, v*) fields
+      REAL(wp) :: zhdivmin,zhdivmax             ! extrema of horizontal divergence of (u*, v*) fields
+      REAL(wp) :: zmrotmin,zmrotmax             ! extrema of the curl of the (u*, v*) fields
+      REAL(wp) :: ustar,vstar                   ! (u*, v*) before tapering in shallow water
+      REAL(wp) :: hramp                         ! depth over which Neptune vel. is ramped down
+      !!
+      NAMELIST/namdyn_nept/ ln_neptsimp,      &
+                            ln_smooth_neptvel,&
+			    rn_tslse,	      &
+			    rn_tslsp,	      &
+                            ln_neptramp,      &
+                            rn_htrmin,        &
+			    rn_htrmax
+      !!----------------------------------------------------------------------
+
+      ! Define the (simplified) Neptune parameters
+      ! ==========================================
+
+!!    WRITE(numout,*) ' start dynnept namelist'
+!!    CALL FLUSH(numout)
+      REWIND( numnam )                  ! Read Namelist namdyn_nept:  Simplified Neptune
+      READ  ( numnam, namdyn_nept )
+!!    WRITE(numout,*) ' dynnept namelist done'
+!!    CALL FLUSH(numout)
+
+      IF(lwp) THEN                      ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dyn_nept_init : Simplified Neptune module enabled'
+         WRITE(numout,*) '~~~~~~~~~~~~~'
+         WRITE(numout,*) ' -->   Reading namelist namdyn_nept parameters:'
+         WRITE(numout,*) '       ln_neptsimp          = ', ln_neptsimp
+         WRITE(numout,*)
+         WRITE(numout,*) '       ln_smooth_neptvel    = ', ln_smooth_neptvel
+         WRITE(numout,*) '       rn_tslse             = ', rn_tslse
+         WRITE(numout,*) '       rn_tslsp             = ', rn_tslsp
+         WRITE(numout,*)
+         WRITE(numout,*) '       ln_neptramp          = ', ln_neptramp
+         WRITE(numout,*) '       rn_htrmin            = ', rn_htrmin
+         WRITE(numout,*) '       rn_htrmax            = ', rn_htrmax
+         WRITE(numout,*)
+         CALL FLUSH(numout)
+      ENDIF
+
+      IF( ln_smooth_neptvel ) THEN
+         IF(lwp) WRITE(numout,*) ' -->   neptune velocities will be smoothed'
+      ELSE
+         IF(lwp) WRITE(numout,*) ' -->   neptune velocities will not be smoothed'
+      ENDIF
+
+      IF( ln_neptsimp ) THEN
+          IF(lwp) WRITE(numout,*) ' -->   ln_neptsimp enabled, solving for U-UN'
+      ELSE
+          IF(lwp) WRITE(numout,*) ' -->   ln_neptsimp disabled'
+          RETURN
+      ENDIF
+
+      IF( ln_neptramp ) THEN
+          IF(lwp) WRITE(numout,*) ' -->   ln_neptramp enabled, ramp down Neptune'
+          IF(lwp) WRITE(numout,*) ' -->   velocity components in shallow water'
+      ELSE
+          IF(lwp) WRITE(numout,*) ' -->   ln_neptramp disabled'
+      ENDIF
+
+
+!!    Perform dynamic allocation of shared module variables
+      IF( dynnept_alloc() /= 0 )   CALL ctl_warn('dynnept_alloc: array allocate failed.')
+
+!!    Dynamically allocate local work arrays
+      ALLOCATE( ht(jpi,jpj), htn(jpi,jpj), tscale(jpi,jpj), tsp(jpi,jpj),      &
+         &      hur_n(jpi,jpj), hvr_n(jpi,jpj), hu_n(jpi,jpj), hv_n(jpi,jpj),  &
+         &      znmask(jpi,jpj,jpk) )
+
+      IF( .not. ln_rstart ) THEN      ! If restarting, these arrays are read from the restart file
+         zhdivnep(:,:,:) = 0.0_wp
+         zmrotnep(:,:,:) = 0.0_wp
+      END IF
+
+      ! Computation of nmask: same as fmask, but fmask cannot be used
+      ! because it is modified after it is computed in dom_msk
+      ! (this can be optimised to save memory, such as merge into next loop)
+      DO jk = 1, jpk
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector loop
+               znmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
+                   &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
+            END DO
+         END DO
+      END DO
+
+      CALL lbc_lnk( znmask, 'F', 1.0_wp )
+
+
+      ! now compute zunep, zvnep (renamed from earlier versions)
+
+      zunep(:,:) = 0.0_wp
+      zvnep(:,:) = 0.0_wp
+
+      htn(:,:) = 0.0_wp            ! ocean depth at F-point
+      DO jk = 1, jpk
+         htn(:,:) = htn(:,:) + fse3f(:,:,jk) * znmask(:,:,jk)
+      END DO
+
+      IF( ln_smooth_neptvel ) THEN
+         CALL dyn_nept_smooth_vel( htn, ht, .TRUE. )
+      !! overwrites ht with a smoothed version of htn
+      ELSE
+         ht(:,:) = htn(:,:)
+      !! use unsmoothed version of htn
+      ENDIF
+      CALL lbc_lnk( ht, 'F', 1.0_wp )
+
+      !! Compute tsp, a stream function for the Neptune velocity,
+      !! with the usual geophysical sign convention
+      !! Then zunep = -latitudinal derivative "-(1/H)*d(tsp)/dy"
+      !!      zvnep = longitudinal derivative " (1/H)*d(tsp)/dx"
+
+      tsp(:,:)    = 0.0_wp
+      tscale(:,:) = 0.0_wp
+
+      tscale(:,:) = rn_tslsp + (rn_tslse - rn_tslsp) *   &
+                   ( 0.5_wp + 0.5_wp * COS( 2.0_wp * rad * gphif(:,:) )  )
+      tsp   (:,:) = -2.0_wp * omega * SIN( rad * gphif(:,:) ) * tscale(:,:) * tscale(:,:) * ht(:,:)
+
+
+      IF( ln_smooth_neptvel ) THEN
+         CALL dyn_nept_smooth_vel( hu, hu_n, .TRUE. )
+      !! overwrites hu_n with a smoothed version of hu
+      ELSE
+         hu_n(:,:) = hu(:,:)
+      !! use unsmoothed version of hu
+      ENDIF
+      CALL lbc_lnk( hu_n, 'U', 1.0_wp )
+      hu_n(:,:) = hu_n(:,:) * umask(:,:,1)
+
+      WHERE( hu_n(:,:) == 0.0_wp )
+         hur_n(:,:) = 0.0_wp
+      ELSEWHERE
+         hur_n(:,:) = 1.0_wp / hu_n(:,:)
+      END WHERE
+
+
+      IF( ln_smooth_neptvel ) THEN
+         CALL dyn_nept_smooth_vel( hv, hv_n, .TRUE. )
+      !! overwrites hv_n with a smoothed version of hv
+      ELSE
+         hv_n(:,:) = hv(:,:)
+      !! use unsmoothed version of hv
+      ENDIF
+      CALL lbc_lnk( hv_n, 'V', 1.0_wp )
+      hv_n(:,:) = hv_n(:,:) * vmask(:,:,1)
+
+      WHERE( hv_n == 0.0_wp )
+         hvr_n(:,:) = 0.0_wp
+      ELSEWHERE
+         hvr_n(:,:) = 1.0_wp / hv_n(:,:)
+      END WHERE
+
+
+      unemin =  1.0e35
+      unemax = -1.0e35
+      vnemin =  1.0e35
+      vnemax = -1.0e35
+      hramp = rn_htrmax - rn_htrmin
+      DO jj = 2, jpj-1
+         DO ji = 2, jpi-1
+            if ( umask(ji,jj,1) /= 0.0_wp ) then
+               ustar =-1.0_wp/e2u(ji,jj) * hur_n(ji,jj) * ( tsp(ji,jj)-tsp(ji,jj-1) ) * umask(ji,jj,1)
+               if ( ln_neptramp ) then
+!!                Apply ramp down to velocity component
+                  if ( hu_n(ji,jj) <= rn_htrmin ) then
+                    zunep(ji,jj) = 0.0_wp
+                   else if ( hu_n(ji,jj) >= rn_htrmax ) then
+                    zunep(ji,jj) = ustar
+                   else if ( hramp > 0.0_wp ) then
+                    zunep(ji,jj) = ( hu_n(ji,jj) - rn_htrmin) * ustar/hramp
+                  endif
+                else
+                 zunep(ji,jj) = ustar
+               endif 
+             else
+              zunep(ji,jj) = 0.0_wp
+            endif
+            if ( vmask(ji,jj,1) /= 0.0_wp ) then
+               vstar = 1.0_wp/e1v(ji,jj) * hvr_n(ji,jj) * ( tsp(ji,jj)-tsp(ji-1,jj) ) * vmask(ji,jj,1)
+               if ( ln_neptramp ) then
+!!                Apply ramp down to velocity component
+                  if ( hv_n(ji,jj) <= rn_htrmin ) then
+                    zvnep(ji,jj) = 0.0_wp
+                   else if ( hv_n(ji,jj) >= rn_htrmax ) then
+                    zvnep(ji,jj) = vstar
+                   else if ( hramp > 0.0_wp ) then
+                    zvnep(ji,jj) = ( hv_n(ji,jj) - rn_htrmin) * vstar/hramp
+                  endif
+                else
+                  zvnep(ji,jj) = vstar
+               endif
+             else
+              zvnep(ji,jj) = 0.0_wp
+            endif
+            unemin = min( unemin, zunep(ji,jj) )
+            unemax = max( unemax, zunep(ji,jj) )
+            vnemin = min( vnemin, zvnep(ji,jj) )
+            vnemax = max( vnemax, zvnep(ji,jj) )
+         END DO
+      END DO
+      CALL lbc_lnk( zunep, 'U', -1.0_wp )
+      CALL lbc_lnk( zvnep, 'V', -1.0_wp )
+      WRITE(numout,*) '      zunep: min, max       = ', unemin,unemax
+      WRITE(numout,*) '      zvnep: min, max       = ', vnemin,vnemax
+      WRITE(numout,*)
+
+      !!  Compute, once and for all, the horizontal divergence (zhdivnep)
+      !!  and the curl (zmrotnep) of the Neptune velocity field (zunep, zvnep)
+      CALL dyn_nept_div_cur_init
+
+      !! Check the ranges of the computed divergence & vorticity
+      zhdivmin =  1.0e35
+      zhdivmax = -1.0e35
+      zmrotmin =  1.0e35
+      zmrotmax = -1.0e35
+      hramp = rn_htrmax - rn_htrmin
+      DO jk = 1, jpkm1                                 ! Horizontal slab
+         DO jj = 2, jpj-1
+            DO ji = 2, jpi-1
+               zhdivmin = min( zhdivmin, zhdivnep(ji,jj,jk) )
+               zhdivmax = max( zhdivmax, zhdivnep(ji,jj,jk) )
+               zmrotmin = min( zmrotmin, zmrotnep(ji,jj,jk) )
+               zmrotmax = max( zmrotmax, zmrotnep(ji,jj,jk) )
+            END DO
+         END DO
+      END DO
+      WRITE(numout,*) '   zhdivnep: min, max       = ', zhdivmin,zhdivmax
+      WRITE(numout,*) '   zmrotnep: min, max       = ', zmrotmin,zmrotmax
+      WRITE(numout,*)
+
+!!    Deallocate temporary workspace arrays, which are all local to
+!!    this routine, except where passed as arguments to other routines
+      DEALLOCATE( ht, htn, tscale, tsp, hur_n, hvr_n, hu_n, hv_n, znmask )
+
+   END SUBROUTINE dyn_nept_init
+
+
+   SUBROUTINE dyn_nept_div_cur_init
+      !!----------------------------------------------------------------------
+      !!             ***  ROUTINE dyn_nept_div_cur_init  ***
+      !!
+      !! ** Purpose :   compute the horizontal divergence and the relative
+      !!                vorticity of the time-invariant u* and v* Neptune
+      !!                effect velocities (called zunep, zvnep)
+      !!
+      !! ** Method  : - Divergence:
+      !!      - compute the divergence given by :
+      !!         zhdivnep = 1/(e1t*e2t*e3t) ( di[e2u*e3u zunep] + dj[e1v*e3v zvnep] )
+      !!      - compute the curl in tensorial formalism:
+      !!         zmrotnep = 1/(e1f*e2f) ( di[e2v zvnep] - dj[e1u zunep] )
+      !!      Note: Coastal boundary condition: lateral friction set through
+      !!      the value of fmask along the coast (see dommsk.F90) and shlat
+      !!      (namelist parameter)
+      !!
+      !! ** Action  : - compute zhdivnep, the hor. divergence of (u*, v*)
+      !!              - compute zmrotnep, the rel. vorticity  of (u*, v*)
+      !!
+      !! History :  OPA  ! 1987-06  (P. Andrich, D. L Hostis)  Original code
+      !!            4.0  ! 1991-11  (G. Madec)
+      !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions
+      !!            7.0  ! 1996-01  (G. Madec)  s-coordinates
+      !!            8.0  ! 1997-06  (G. Madec)  lateral boundary cond., lbc
+      !!            8.1  ! 1997-08  (J.M. Molines)  Open boundaries
+      !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
+      !!  NEMO      1.0  ! 2002-09  (G. Madec, E. Durand)  Free form, F90
+      !!             -   ! 2005-01  (J. Chanut) Unstructured open boundaries
+      !!             -   ! 2003-08  (G. Madec)  merged of cur and div, free form, F90
+      !!             -   ! 2005-01  (J. Chanut, A. Sellar) unstructured open boundaries
+      !!            3.3  ! 2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
+      !!                 ! 2011-06  (Jeff Blundell, NOCS) Adapt code from divcur.F90
+      !!                 !           to compute Neptune effect fields only
+      !!----------------------------------------------------------------------
+      !
+      INTEGER  ::   ji, jj, jk          ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'dyn_nept_div_cur_init :'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~'
+      IF(lwp) WRITE(numout,*) 'horizontal velocity divergence and'
+      IF(lwp) WRITE(numout,*) 'relative vorticity of Neptune flow'
+#if defined key_noslip_accurate
+   !!----------------------------------------------------------------------
+   !!   'key_noslip_accurate'                     2nd order centered scheme
+   !!                                                4th order at the coast
+   !!----------------------------------------------------------------------
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'WARNING: key_noslip_accurate option'
+      IF(lwp) WRITE(numout,*) 'not implemented in simplified Neptune'
+      CALL ctl_warn( ' noslip_accurate option not implemented' )
+#endif
+
+   !!----------------------------------------------------------------------
+   !!   Default option                           2nd order centered schemes
+   !!----------------------------------------------------------------------
+
+      ! Apply the div and curl operators to the depth-dependent velocity
+      ! field produced by multiplying (zunep, zvnep) by (umask, vmask), exactly
+      ! equivalent to the equivalent calculation in the unsimplified code
+      !                                                ! ===============
+      DO jk = 1, jpkm1                                 ! Horizontal slab
+         !                                             ! ===============
+         !                                             ! --------
+         ! Horizontal divergence                       !   div
+         !                                             ! --------
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zhdivnep(ji,jj,jk) =   &
+               &   (  e2u(ji  ,jj  )*fse3u(ji  ,jj  ,jk) * zunep(ji  ,jj  ) * umask(ji  ,jj  ,jk)    &
+               &    - e2u(ji-1,jj  )*fse3u(ji-1,jj  ,jk) * zunep(ji-1,jj  ) * umask(ji-1,jj  ,jk)    &
+               &    + e1v(ji  ,jj  )*fse3v(ji  ,jj  ,jk) * zvnep(ji  ,jj  ) * vmask(ji  ,jj  ,jk)    &
+               &    - e1v(ji  ,jj-1)*fse3v(ji  ,jj-1,jk) * zvnep(ji  ,jj-1) * vmask(ji  ,jj-1,jk) )  &
+               &  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+            END DO
+         END DO
+
+#if defined key_obc
+         IF( Agrif_Root() ) THEN
+            ! open boundaries (div must be zero behind the open boundary)
+            !  mpp remark: The zeroing of zhdivnep can probably be extended to 1->jpi/jpj for the correct row/column
+            IF( lp_obc_east  )  zhdivnep(nie0p1:nie1p1,nje0  :nje1  ,jk) = 0.0_wp      ! east
+            IF( lp_obc_west  )  zhdivnep(niw0  :niw1  ,njw0  :njw1  ,jk) = 0.0_wp      ! west
+            IF( lp_obc_north )  zhdivnep(nin0  :nin1  ,njn0p1:njn1p1,jk) = 0.0_wp      ! north
+            IF( lp_obc_south )  zhdivnep(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.0_wp      ! south
+         ENDIF
+#endif
+         IF( .NOT. AGRIF_Root() ) THEN
+            IF ((nbondi ==  1).OR.(nbondi == 2))  zhdivnep(nlci-1 , :     ,jk) = 0.0_wp   ! east
+            IF ((nbondi == -1).OR.(nbondi == 2))  zhdivnep(2      , :     ,jk) = 0.0_wp   ! west
+            IF ((nbondj ==  1).OR.(nbondj == 2))  zhdivnep(:      ,nlcj-1 ,jk) = 0.0_wp   ! north
+            IF ((nbondj == -1).OR.(nbondj == 2))  zhdivnep(:      ,2      ,jk) = 0.0_wp   ! south
+         ENDIF
+
+         !                                             ! --------
+         ! relative vorticity                          !   rot
+         !                                             ! --------
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zmrotnep(ji,jj,jk) =   &
+                  &       (  e2v(ji+1,jj  ) * zvnep(ji+1,jj  ) * vmask(ji+1,jj  ,jk)     &
+                  &        - e2v(ji  ,jj  ) * zvnep(ji  ,jj  ) * vmask(ji  ,jj  ,jk)     &
+                  &        - e1u(ji  ,jj+1) * zunep(ji  ,jj+1) * umask(ji  ,jj+1,jk)     &
+                  &        + e1u(ji  ,jj  ) * zunep(ji  ,jj  ) * umask(ji  ,jj  ,jk)  )  &
+                  &       * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
+            END DO
+         END DO
+         !                                             ! ===============
+      END DO                                           !   End of slab
+      !                                                ! ===============
+
+      ! 4. Lateral boundary conditions on zhdivnep and zmrotnep
+      ! ----------------------------------=======-----=======
+      CALL lbc_lnk( zhdivnep, 'T', 1. )   ;   CALL lbc_lnk( zmrotnep , 'F', 1. )     ! lateral boundary cond. (no sign change)
+      !
+   END SUBROUTINE dyn_nept_div_cur_init
+
+
+   SUBROUTINE dyn_nept_cor( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_nept_cor  ***
+      !!
+      !! ** Purpose :  Add or subtract the Neptune velocity from the now velocities
+      !!
+      !! ** Method  :  First call : kt not equal to lastkt -> subtract zunep, zvnep
+      !!               Second call: kt     equal to lastkt -> add zunep, zvnep
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
+      !!
+      INTEGER, SAVE :: lastkt             ! store previous kt
+      DATA lastkt/-1/                     ! initialise previous kt
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_neptsimp ) THEN
+         !
+         IF( lastkt /= kt ) THEN          ! 1st call for this kt: subtract the Neptune velocities zunep, zvnep from un, vn
+            CALL dyn_nept_vel( -1 )      ! -1 = subtract
+            !
+         ELSE                              ! 2nd call for this kt: add the Neptune velocities zunep, zvnep to un, vn
+            CALL dyn_nept_vel(  1 )      !  1 = add
+            !
+         ENDIF
+         !
+	 lastkt = kt	     ! Store kt
+	 !
+      ENDIF
+      !
+   END SUBROUTINE dyn_nept_cor
+
+
+   SUBROUTINE dyn_nept_vel( ksign )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_nept_vel  ***
+      !!
+      !! ** Purpose :  Add or subtract the Neptune velocity from the now
+      !!               velocities based on ksign
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   ksign    ! 1 or -1 to add or subtract neptune velocities
+      !!
+      INTEGER :: jk                       ! dummy loop index
+      !!----------------------------------------------------------------------
+      !
+      ! Adjust the current velocity un, vn by adding or subtracting the
+      ! Neptune velocities zunep, zvnep, as determined by argument ksign
+      DO jk=1, jpk
+         un(:,:,jk) = un(:,:,jk) + ksign * zunep(:,:) * umask(:,:,jk)
+         vn(:,:,jk) = vn(:,:,jk) + ksign * zvnep(:,:) * vmask(:,:,jk)
+      END DO
+      !
+   END SUBROUTINE dyn_nept_vel
+
+
+   SUBROUTINE dyn_nept_smooth_vel( htold, htnew, option )
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_nept_smooth_vel  ***
+      !!
+      !! ** Purpose :  Compute smoothed topography field.
+      !!
+      !! ** Action : - Updates the array htnew (output) with a smoothed
+      !!               version of the (input) array htold. Form of smoothing
+      !!               algorithm is controlled by the (logical) argument option.
+      !!----------------------------------------------------------------------
+
+      INTEGER                                   ::   ji, jj  ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN)  ::   htold   ! temporary 2D workspace
+      LOGICAL, INTENT(IN)                       ::   option  ! temporary 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)              ::   htnew   ! temporary 2D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   work    ! temporary 2D workspace
+      INTEGER,  ALLOCATABLE, DIMENSION(:,:)     ::   nb      ! temporary 2D workspace
+      INTEGER,  ALLOCATABLE, DIMENSION(:,:)     ::   iwork   ! temporary 2D workspace
+
+!!    Dynamically allocate local work arrays
+      ALLOCATE( work(jpi,jpj), nb(jpi,jpj), iwork(jpi,jpj) )
+
+      iwork(:,:) = 0
+
+      !! iwork is a mask of gridpoints: iwork = 1 => ocean, iwork = 0 => land
+      WHERE( htold(:,:) > 0 )
+         iwork(:,:) = 1
+         htnew(:,:) = htold(:,:)
+      ELSEWHERE
+         iwork(:,:) = 0
+         htnew(:,:) = 0.0_wp
+      END WHERE
+      !! htnew contains valid ocean depths from htold, or zero
+
+      !! set work to a smoothed/averaged version of htnew; choice controlled by option
+      !! nb is set to the sum of the weights of the valid values used in work
+      IF( option ) THEN
+
+         !! Apply scale-selective smoothing in determining work from htnew
+         DO jj=2,jpj-1
+            DO ji=2,jpi-1
+               work(ji,jj) = 4.0*htnew( ji , jj ) +                        &
+                           & 2.0*htnew(ji+1, jj ) + 2.0*htnew(ji-1, jj ) + &
+                           & 2.0*htnew( ji ,jj+1) + 2.0*htnew( ji ,jj-1) + &
+                           &     htnew(ji+1,jj+1) +     htnew(ji+1,jj-1) + &
+                           &     htnew(ji-1,jj+1) +     htnew(ji-1,jj-1)
+
+               nb(ji,jj)   = 4 * iwork( ji , jj ) +                        &
+                           & 2 * iwork(ji+1, jj ) + 2 * iwork(ji-1, jj ) + &
+                           & 2 * iwork( ji ,jj+1) + 2 * iwork( ji ,jj-1) + &
+                           &     iwork(ji+1,jj+1) +     iwork(ji+1,jj-1) + &
+                           &     iwork(ji-1,jj+1) +     iwork(ji-1,jj-1)
+            END DO
+         END DO
+
+      ELSE
+
+         !! Apply simple 9-point averaging in determining work from htnew
+         DO jj=2,jpj-1
+            DO ji=2,jpi-1
+               work(ji,jj) =  htnew( ji , jj ) +                    &
+                           &  htnew(ji+1, jj ) + htnew(ji-1, jj ) + &
+                           &  htnew( ji ,jj+1) + htnew( ji ,jj-1) + &
+                           &  htnew(ji+1,jj+1) + htnew(ji+1,jj-1) + &
+                           &  htnew(ji-1,jj+1) + htnew(ji-1,jj-1)
+
+               nb(ji,jj) =    iwork( ji , jj ) +                    &
+                           &  iwork(ji+1, jj ) + iwork(ji-1, jj ) + &
+                           &  iwork( ji ,jj+1) + iwork( ji ,jj-1) + &
+                           &  iwork(ji+1,jj+1) + iwork(ji+1,jj-1) + &
+                           &  iwork(ji-1,jj+1) + iwork(ji-1,jj-1)
+            END DO
+         END DO
+
+      ENDIF
+
+      !! write averaged value of work into htnew,
+      !! if average is valid and point is unmasked
+      WHERE( (htold(:,:) /= 0.0_wp ) .AND. ( nb(:,:) /= 0 ) )
+         htnew(:,:) = work(:,:)/real(nb(:,:))
+      ELSEWHERE
+         htnew(:,:) = 0.0_wp
+      END WHERE
+
+!!    Deallocate temporary workspace arrays, all local to this routine
+      DEALLOCATE( work, nb, iwork )
+
+   END SUBROUTINE dyn_nept_smooth_vel
+
+END MODULE dynnept
Index: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 2985)
+++ branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 2986)
@@ -66,4 +66,5 @@
    USE c1d             ! 1D configuration
    USE step_c1d        ! Time stepping loop for the 1D configuration
+   USE dynnept         ! simplified form of Neptune effect
 #if defined key_top
    USE trcini          ! passive tracer initialisation
@@ -296,4 +297,8 @@
       IF( lk_obc        )   CALL     obc_init   ! Open boundaries 
       IF( lk_bdy        )   CALL     bdy_init   ! Unstructured open boundaries
+
+                            CALL flush(numout)
+                            CALL dyn_nept_init  ! simplified form of Neptune effect
+                            CALL flush(numout)
 
                             CALL  istate_init   ! ocean initial state (Dynamics and tracers)
Index: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 2985)
+++ branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 2986)
@@ -36,4 +36,5 @@
 #endif
    USE asminc           ! assimilation increments    (tra_asm_inc, dyn_asm_inc routines)
+   USE dynnept          ! simplified form of Neptune effect
 
    IMPLICIT NONE
@@ -220,7 +221,9 @@
       IF(  ln_asmiau .AND. &
          & ln_dyninc       )   CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment
+      IF( ln_neptsimp )        CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified)
                                CALL dyn_adv( kstp )         ! advection (vector or flux form)
                                CALL dyn_vor( kstp )         ! vorticity term including Coriolis
                                CALL dyn_ldf( kstp )         ! lateral mixing
+      IF( ln_neptsimp )        CALL dyn_nept_cor( kstp )    ! add Neptune velocities (simplified)
 #if defined key_agrif
       IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_dyn        ! momemtum sponge
