ecearth3/sources/nemo-3.6/NEMO/OPA_SRC/DYN/dynnept.F90

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
   USE wrk_nemo        ! Memory Allocation

   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          ! yes/no simplified neptune

   LOGICAL          ::  ln_smooth_neptvel    ! yes/no smooth zunep, zvnep
   REAL(wp)         ::  rn_tslse             ! value of lengthscale L at the equator
   REAL(wp)         ::  rn_tslsp             ! 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          ! ramp down Neptune velocity in shallow water
   REAL(wp)         ::  rn_htrmin            ! min. depth of transition range
   REAL(wp)         ::  rn_htrmax            ! 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
   !!----------------------------------------------------------------------

   !! $Id: dynnept.F90 2355 2015-05-20 07:11:50Z ufla $
 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
      !!
      INTEGER  ::   ji, jj, jk    ! dummy loop indices
      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
      !
      REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, htn, tscale, tsp, hur_n, hvr_n, hu_n, hv_n      
      REAL(wp), POINTER, DIMENSION(:,:,:) :: znmask
      !!
      NAMELIST/namdyn_nept/ ln_neptsimp, ln_smooth_neptvel, rn_tslse, rn_tslsp,      &
                            ln_neptramp, rn_htrmin, rn_htrmax
      INTEGER  ::   ios
      !!----------------------------------------------------------------------
      ! Define the (simplified) Neptune parameters
      ! ==========================================

      REWIND( numnam_ref )              ! Namelist namdyn_nept in reference namelist : Simplified Neptune
      READ  ( numnam_ref, namdyn_nept, IOSTAT = ios, ERR = 901)
901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_nept in reference namelist', lwp )

      REWIND( numnam_cfg )              ! Namelist namdyn_nept in reference namelist : Simplified Neptune
      READ  ( numnam_cfg, namdyn_nept, IOSTAT = ios, ERR = 902 )
902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_nept in configuration namelist', lwp )
      IF(lwm) WRITE ( numond, namdyn_nept )

      IF(lwp) THEN                      ! Control print
         WRITE(numout,*)
         WRITE(numout,*) 'dyn_nept_init : Simplified Neptune module'
         WRITE(numout,*) '~~~~~~~~~~~~~'
         WRITE(numout,*) ' -->   Reading namelist namdyn_nept parameters:'
         WRITE(numout,*) '       ln_neptsimp          = ', ln_neptsimp
         WRITE(numout,*)
         IF( ln_neptsimp ) THEN
            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,*)
         ENDIF
      ENDIF
      !
      IF( .NOT. ln_neptsimp ) RETURN
      !                                 ! Dynamically allocate local work arrays
      CALL wrk_alloc( jpi, jpj     , zht, htn, tscale, tsp, hur_n, hvr_n, hu_n, hv_n  ) 
      CALL wrk_alloc( jpi, jpj, jpk, znmask                                          ) 

      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_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.')

      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, zht, .TRUE. )
      !! overwrites zht with a smoothed version of htn
      ELSE
         zht(:,:) = htn(:,:)
      !! use unsmoothed version of htn
      ENDIF
      CALL lbc_lnk( zht, '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(:,:) * zht(:,:)


      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
      CALL wrk_dealloc( jpi, jpj     , zht, htn, tscale, tsp, hur_n, hvr_n, hu_n, hv_n  ) 
      CALL wrk_dealloc( jpi, jpj, jpk, 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( .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, ld_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 ld_option.
      !!----------------------------------------------------------------------
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )  ::  htold      ! temporary 2D workspace
      LOGICAL                     , INTENT(in   )  ::  ld_option  ! temporary 2D workspace
      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)  ::  htnew      ! temporary 2D workspace
      !
      INTEGER                           ::  ji, jj  ! dummy loop indices
      INTEGER , POINTER, DIMENSION(:,:) ::  nb, iwork
      REAL(wp), POINTER, DIMENSION(:,:) ::  work    ! temporary 2D workspace
      !!----------------------------------------------------------------------
      !
      CALL wrk_alloc( jpi, jpj, nb, iwork ) 
      CALL wrk_alloc( jpi, jpj, work      ) 
      !
      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 ld_option
      !! nb is set to the sum of the weights of the valid values used in work
      IF( ld_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
      CALL wrk_dealloc( jpi, jpj, nb, iwork ) 
      CALL wrk_dealloc( jpi, jpj, work      ) 
      !
   END SUBROUTINE dyn_nept_smooth_vel

END MODULE dynnept