New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynnept.F90 in branches/UKMO/dev_r5107_mld_zint/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/dev_r5107_mld_zint/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90 @ 5450

Last change on this file since 5450 was 5450, checked in by davestorkey, 9 years ago

Clear SVn keywords from UKMO/dev_r5107_mld_zint branch.

File size: 27.1 KB
Line 
1MODULE dynnept
2   !!======================================================================
3   !!                       ***  MODULE  dynnept  ***
4   !! Ocean dynamics: Neptune effect as proposed by Greg Holloway,
5   !!                 recoded version of simplest case (u*, v* only)
6   !!======================================================================
7   !! History :  1.0  !  2007-06  (Michael Dunphy)  Modular form: - new namelist parameters
8   !!                                                             - horizontal diffusion for Neptune
9   !!                                                             - vertical diffusion for gm in momentum eqns
10   !!                                                             - option to use Neptune in Coriolis eqn
11   !!                    2011-08  (Jeff Blundell, NOCS) Simplified form for temporally invariant u*, v*
12   !!                                               Horizontal and vertical diffusivity formulations removed
13   !!                                               Dynamic allocation of storage added
14   !!                                               Option of ramping Neptune vel. down
15   !!                                               to zero added in shallow depths added
16   !!----------------------------------------------------------------------
17   !! dynnept_alloc        :
18   !! dyn_nept_init        :
19   !! dyn_nept_div_cur_init:
20   !! dyn_nept_cor         :
21   !! dyn_nept_vel         :
22   !! dyn_nept_smooth_vel  :
23   !!----------------------------------------------------------------------
24   USE oce              ! ocean dynamics and tracers
25   USE dom_oce          ! ocean space and time domain
26   USE in_out_manager   ! I/O manager
27   USE lib_mpp          ! distributed memory computing
28   USE prtctl           ! Print control
29   USE phycst
30   USE lbclnk
31   USE wrk_nemo        ! Memory Allocation
32
33   IMPLICIT NONE
34   PRIVATE
35
36   !! * Routine accessibility
37   PUBLIC dyn_nept_init      ! routine called by nemogcm.F90
38   PUBLIC dyn_nept_cor       ! routine called by step.F90
39   !! dynnept_alloc()         is called only by dyn_nept_init, within this module
40   !! dyn_nept_div_cur_init   is called only by dyn_nept_init, within this module
41   !! dyn_nept_vel            is called only by dyn_nept_cor,  within this module
42
43   !! * Shared module variables
44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: zunep, zvnep  ! Neptune u and v
45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: zhdivnep      ! hor. div for Neptune vel.
46   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: zmrotnep      ! curl for Neptune vel.
47
48
49   !! * Namelist namdyn_nept variables
50   LOGICAL, PUBLIC  ::  ln_neptsimp          ! yes/no simplified neptune
51
52   LOGICAL          ::  ln_smooth_neptvel    ! yes/no smooth zunep, zvnep
53   REAL(wp)         ::  rn_tslse             ! value of lengthscale L at the equator
54   REAL(wp)         ::  rn_tslsp             ! value of lengthscale L at the pole
55!! Specify whether to ramp down the Neptune velocity in shallow
56!! water, and the depth range controlling such ramping down
57   LOGICAL          ::  ln_neptramp          ! ramp down Neptune velocity in shallow water
58   REAL(wp)         ::  rn_htrmin            ! min. depth of transition range
59   REAL(wp)         ::  rn_htrmax            ! max. depth of transition range
60
61   !! * Module variables
62
63
64   !! * Substitutions
65#  include "vectopt_loop_substitute.h90"
66#  include "domzgr_substitute.h90"
67   !!----------------------------------------------------------------------
68   !!   OPA 9.0 , implemented by Bedford Institute of Oceanography
69   !!----------------------------------------------------------------------
70
71   !! $Id$
72 CONTAINS
73
74   INTEGER FUNCTION dynnept_alloc()
75      !!----------------------------------------------------------------------
76      !!                    ***  ROUTINE dynnept_alloc  ***
77      !!----------------------------------------------------------------------
78      ALLOCATE( zunep(jpi,jpj) , zvnep(jpi,jpj) ,     &
79         &      zhdivnep(jpi,jpj,jpk) , zmrotnep(jpi,jpj,jpk) , STAT=dynnept_alloc )
80         !
81      IF( dynnept_alloc /= 0 )   CALL ctl_warn('dynnept_alloc: array allocate failed.')
82   END FUNCTION dynnept_alloc
83
84
85   SUBROUTINE dyn_nept_init
86      !!----------------------------------------------------------------------
87      !!                  ***  ROUTINE dyn_nept_init  ***
88      !!
89      !! ** Purpose :   Read namelist parameters, initialise arrays
90      !!                and compute the arrays zunep and zvnep
91      !!
92      !! ** Method  :   zunep =
93      !!                zvnep =
94      !!
95      !! ** History :  1.0  !   07-05  (Zeliang Wang)   Original code for zunep, zvnep
96      !!               1.1  !   07-06  (Michael Dunphy) namelist and  initialisation
97      !!               2.0  ! 2011-07  (Jeff Blundell, NOCS)
98      !!                    ! Simplified form for temporally invariant u*, v*
99      !!                    ! Horizontal and vertical diffusivity formulations removed
100      !!                    ! Includes optional tapering-off in shallow depths
101      !!----------------------------------------------------------------------
102      USE iom
103      !!
104      INTEGER  ::   ji, jj, jk    ! dummy loop indices
105      REAL(wp) :: unemin,unemax,vnemin,vnemax   ! extrema of (u*, v*) fields
106      REAL(wp) :: zhdivmin,zhdivmax             ! extrema of horizontal divergence of (u*, v*) fields
107      REAL(wp) :: zmrotmin,zmrotmax             ! extrema of the curl of the (u*, v*) fields
108      REAL(wp) :: ustar,vstar                   ! (u*, v*) before tapering in shallow water
109      REAL(wp) :: hramp                         ! depth over which Neptune vel. is ramped down
110      !
111      REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, htn, tscale, tsp, hur_n, hvr_n, hu_n, hv_n     
112      REAL(wp), POINTER, DIMENSION(:,:,:) :: znmask
113      !!
114      NAMELIST/namdyn_nept/ ln_neptsimp, ln_smooth_neptvel, rn_tslse, rn_tslsp,      &
115                            ln_neptramp, rn_htrmin, rn_htrmax
116      INTEGER  ::   ios
117      !!----------------------------------------------------------------------
118      ! Define the (simplified) Neptune parameters
119      ! ==========================================
120
121      REWIND( numnam_ref )              ! Namelist namdyn_nept in reference namelist : Simplified Neptune
122      READ  ( numnam_ref, namdyn_nept, IOSTAT = ios, ERR = 901)
123901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_nept in reference namelist', lwp )
124
125      REWIND( numnam_cfg )              ! Namelist namdyn_nept in reference namelist : Simplified Neptune
126      READ  ( numnam_cfg, namdyn_nept, IOSTAT = ios, ERR = 902 )
127902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_nept in configuration namelist', lwp )
128      IF(lwm) WRITE ( numond, namdyn_nept )
129
130      IF(lwp) THEN                      ! Control print
131         WRITE(numout,*)
132         WRITE(numout,*) 'dyn_nept_init : Simplified Neptune module'
133         WRITE(numout,*) '~~~~~~~~~~~~~'
134         WRITE(numout,*) ' -->   Reading namelist namdyn_nept parameters:'
135         WRITE(numout,*) '       ln_neptsimp          = ', ln_neptsimp
136         WRITE(numout,*)
137         IF( ln_neptsimp ) THEN
138            WRITE(numout,*) '       ln_smooth_neptvel    = ', ln_smooth_neptvel
139            WRITE(numout,*) '       rn_tslse             = ', rn_tslse
140            WRITE(numout,*) '       rn_tslsp             = ', rn_tslsp
141            WRITE(numout,*)
142            WRITE(numout,*) '       ln_neptramp          = ', ln_neptramp
143            WRITE(numout,*) '       rn_htrmin            = ', rn_htrmin
144            WRITE(numout,*) '       rn_htrmax            = ', rn_htrmax
145            WRITE(numout,*)
146         ENDIF
147      ENDIF
148      !
149      IF( .NOT. ln_neptsimp ) RETURN
150      !                                 ! Dynamically allocate local work arrays
151      CALL wrk_alloc( jpi, jpj     , zht, htn, tscale, tsp, hur_n, hvr_n, hu_n, hv_n  ) 
152      CALL wrk_alloc( jpi, jpj, jpk, znmask                                          ) 
153
154      IF( ln_smooth_neptvel ) THEN
155         IF(lwp) WRITE(numout,*) ' -->   neptune velocities will be smoothed'
156      ELSE
157         IF(lwp) WRITE(numout,*) ' -->   neptune velocities will not be smoothed'
158      ENDIF
159
160      IF( ln_neptramp ) THEN
161          IF(lwp) WRITE(numout,*) ' -->   ln_neptramp enabled, ramp down Neptune'
162          IF(lwp) WRITE(numout,*) ' -->   velocity components in shallow water'
163      ELSE
164          IF(lwp) WRITE(numout,*) ' -->   ln_neptramp disabled'
165      ENDIF
166
167
168!!    Perform dynamic allocation of shared module variables
169      IF( dynnept_alloc() /= 0 )   CALL ctl_warn('dynnept_alloc: array allocate failed.')
170
171      IF( .not. ln_rstart ) THEN      ! If restarting, these arrays are read from the restart file
172         zhdivnep(:,:,:) = 0.0_wp
173         zmrotnep(:,:,:) = 0.0_wp
174      END IF
175
176      ! Computation of nmask: same as fmask, but fmask cannot be used
177      ! because it is modified after it is computed in dom_msk
178      ! (this can be optimised to save memory, such as merge into next loop)
179      DO jk = 1, jpk
180         DO jj = 1, jpjm1
181            DO ji = 1, fs_jpim1   ! vector loop
182               znmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
183                   &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
184            END DO
185         END DO
186      END DO
187
188      CALL lbc_lnk( znmask, 'F', 1.0_wp )
189
190
191      ! now compute zunep, zvnep (renamed from earlier versions)
192
193      zunep(:,:) = 0.0_wp
194      zvnep(:,:) = 0.0_wp
195
196      htn(:,:) = 0.0_wp            ! ocean depth at F-point
197      DO jk = 1, jpk
198         htn(:,:) = htn(:,:) + fse3f(:,:,jk) * znmask(:,:,jk)
199      END DO
200
201      IF( ln_smooth_neptvel ) THEN
202         CALL dyn_nept_smooth_vel( htn, zht, .TRUE. )
203      !! overwrites zht with a smoothed version of htn
204      ELSE
205         zht(:,:) = htn(:,:)
206      !! use unsmoothed version of htn
207      ENDIF
208      CALL lbc_lnk( zht, 'F', 1.0_wp )
209
210      !! Compute tsp, a stream function for the Neptune velocity,
211      !! with the usual geophysical sign convention
212      !! Then zunep = -latitudinal derivative "-(1/H)*d(tsp)/dy"
213      !!      zvnep = longitudinal derivative " (1/H)*d(tsp)/dx"
214
215      tsp(:,:)    = 0.0_wp
216      tscale(:,:) = 0.0_wp
217
218      tscale(:,:) = rn_tslsp + (rn_tslse - rn_tslsp) *   &
219                   ( 0.5_wp + 0.5_wp * COS( 2.0_wp * rad * gphif(:,:) )  )
220      tsp   (:,:) = -2.0_wp * omega * SIN( rad * gphif(:,:) ) * tscale(:,:) * tscale(:,:) * zht(:,:)
221
222
223      IF( ln_smooth_neptvel ) THEN
224         CALL dyn_nept_smooth_vel( hu, hu_n, .TRUE. )
225      !! overwrites hu_n with a smoothed version of hu
226      ELSE
227         hu_n(:,:) = hu(:,:)
228      !! use unsmoothed version of hu
229      ENDIF
230      CALL lbc_lnk( hu_n, 'U', 1.0_wp )
231      hu_n(:,:) = hu_n(:,:) * umask(:,:,1)
232
233      WHERE( hu_n(:,:) == 0.0_wp )
234         hur_n(:,:) = 0.0_wp
235      ELSEWHERE
236         hur_n(:,:) = 1.0_wp / hu_n(:,:)
237      END WHERE
238
239
240      IF( ln_smooth_neptvel ) THEN
241         CALL dyn_nept_smooth_vel( hv, hv_n, .TRUE. )
242      !! overwrites hv_n with a smoothed version of hv
243      ELSE
244         hv_n(:,:) = hv(:,:)
245      !! use unsmoothed version of hv
246      ENDIF
247      CALL lbc_lnk( hv_n, 'V', 1.0_wp )
248      hv_n(:,:) = hv_n(:,:) * vmask(:,:,1)
249
250      WHERE( hv_n == 0.0_wp )
251         hvr_n(:,:) = 0.0_wp
252      ELSEWHERE
253         hvr_n(:,:) = 1.0_wp / hv_n(:,:)
254      END WHERE
255
256
257      unemin =  1.0e35
258      unemax = -1.0e35
259      vnemin =  1.0e35
260      vnemax = -1.0e35
261      hramp = rn_htrmax - rn_htrmin
262      DO jj = 2, jpj-1
263         DO ji = 2, jpi-1
264            if ( umask(ji,jj,1) /= 0.0_wp ) then
265               ustar =-1.0_wp/e2u(ji,jj) * hur_n(ji,jj) * ( tsp(ji,jj)-tsp(ji,jj-1) ) * umask(ji,jj,1)
266               if ( ln_neptramp ) then
267!!                Apply ramp down to velocity component
268                  if ( hu_n(ji,jj) <= rn_htrmin ) then
269                    zunep(ji,jj) = 0.0_wp
270                   else if ( hu_n(ji,jj) >= rn_htrmax ) then
271                    zunep(ji,jj) = ustar
272                   else if ( hramp > 0.0_wp ) then
273                    zunep(ji,jj) = ( hu_n(ji,jj) - rn_htrmin) * ustar/hramp
274                  endif
275                else
276                 zunep(ji,jj) = ustar
277               endif
278             else
279              zunep(ji,jj) = 0.0_wp
280            endif
281            if ( vmask(ji,jj,1) /= 0.0_wp ) then
282               vstar = 1.0_wp/e1v(ji,jj) * hvr_n(ji,jj) * ( tsp(ji,jj)-tsp(ji-1,jj) ) * vmask(ji,jj,1)
283               if ( ln_neptramp ) then
284!!                Apply ramp down to velocity component
285                  if ( hv_n(ji,jj) <= rn_htrmin ) then
286                    zvnep(ji,jj) = 0.0_wp
287                   else if ( hv_n(ji,jj) >= rn_htrmax ) then
288                    zvnep(ji,jj) = vstar
289                   else if ( hramp > 0.0_wp ) then
290                    zvnep(ji,jj) = ( hv_n(ji,jj) - rn_htrmin) * vstar/hramp
291                  endif
292                else
293                  zvnep(ji,jj) = vstar
294               endif
295             else
296              zvnep(ji,jj) = 0.0_wp
297            endif
298            unemin = min( unemin, zunep(ji,jj) )
299            unemax = max( unemax, zunep(ji,jj) )
300            vnemin = min( vnemin, zvnep(ji,jj) )
301            vnemax = max( vnemax, zvnep(ji,jj) )
302         END DO
303      END DO
304      CALL lbc_lnk( zunep, 'U', -1.0_wp )
305      CALL lbc_lnk( zvnep, 'V', -1.0_wp )
306      WRITE(numout,*) '      zunep: min, max       = ', unemin,unemax
307      WRITE(numout,*) '      zvnep: min, max       = ', vnemin,vnemax
308      WRITE(numout,*)
309
310      !!  Compute, once and for all, the horizontal divergence (zhdivnep)
311      !!  and the curl (zmrotnep) of the Neptune velocity field (zunep, zvnep)
312      CALL dyn_nept_div_cur_init
313
314      !! Check the ranges of the computed divergence & vorticity
315      zhdivmin =  1.0e35
316      zhdivmax = -1.0e35
317      zmrotmin =  1.0e35
318      zmrotmax = -1.0e35
319      hramp = rn_htrmax - rn_htrmin
320      DO jk = 1, jpkm1                                 ! Horizontal slab
321         DO jj = 2, jpj-1
322            DO ji = 2, jpi-1
323               zhdivmin = min( zhdivmin, zhdivnep(ji,jj,jk) )
324               zhdivmax = max( zhdivmax, zhdivnep(ji,jj,jk) )
325               zmrotmin = min( zmrotmin, zmrotnep(ji,jj,jk) )
326               zmrotmax = max( zmrotmax, zmrotnep(ji,jj,jk) )
327            END DO
328         END DO
329      END DO
330      WRITE(numout,*) '   zhdivnep: min, max       = ', zhdivmin,zhdivmax
331      WRITE(numout,*) '   zmrotnep: min, max       = ', zmrotmin,zmrotmax
332      WRITE(numout,*)
333
334!!    Deallocate temporary workspace arrays, which are all local to
335!!    this routine, except where passed as arguments to other routines
336      CALL wrk_dealloc( jpi, jpj     , zht, htn, tscale, tsp, hur_n, hvr_n, hu_n, hv_n  ) 
337      CALL wrk_dealloc( jpi, jpj, jpk, znmask                                          ) 
338      !
339   END SUBROUTINE dyn_nept_init
340
341
342   SUBROUTINE dyn_nept_div_cur_init
343      !!----------------------------------------------------------------------
344      !!             ***  ROUTINE dyn_nept_div_cur_init  ***
345      !!
346      !! ** Purpose :   compute the horizontal divergence and the relative
347      !!                vorticity of the time-invariant u* and v* Neptune
348      !!                effect velocities (called zunep, zvnep)
349      !!
350      !! ** Method  : - Divergence:
351      !!      - compute the divergence given by :
352      !!         zhdivnep = 1/(e1t*e2t*e3t) ( di[e2u*e3u zunep] + dj[e1v*e3v zvnep] )
353      !!      - compute the curl in tensorial formalism:
354      !!         zmrotnep = 1/(e1f*e2f) ( di[e2v zvnep] - dj[e1u zunep] )
355      !!      Note: Coastal boundary condition: lateral friction set through
356      !!      the value of fmask along the coast (see dommsk.F90) and shlat
357      !!      (namelist parameter)
358      !!
359      !! ** Action  : - compute zhdivnep, the hor. divergence of (u*, v*)
360      !!              - compute zmrotnep, the rel. vorticity  of (u*, v*)
361      !!
362      !! History :  OPA  ! 1987-06  (P. Andrich, D. L Hostis)  Original code
363      !!            4.0  ! 1991-11  (G. Madec)
364      !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions
365      !!            7.0  ! 1996-01  (G. Madec)  s-coordinates
366      !!            8.0  ! 1997-06  (G. Madec)  lateral boundary cond., lbc
367      !!            8.1  ! 1997-08  (J.M. Molines)  Open boundaries
368      !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
369      !!  NEMO      1.0  ! 2002-09  (G. Madec, E. Durand)  Free form, F90
370      !!             -   ! 2005-01  (J. Chanut) Unstructured open boundaries
371      !!             -   ! 2003-08  (G. Madec)  merged of cur and div, free form, F90
372      !!             -   ! 2005-01  (J. Chanut, A. Sellar) unstructured open boundaries
373      !!            3.3  ! 2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
374      !!                 ! 2011-06  (Jeff Blundell, NOCS) Adapt code from divcur.F90
375      !!                 !           to compute Neptune effect fields only
376      !!----------------------------------------------------------------------
377      !
378      INTEGER  ::   ji, jj, jk          ! dummy loop indices
379      !!----------------------------------------------------------------------
380      !   
381      IF(lwp) WRITE(numout,*)
382      IF(lwp) WRITE(numout,*) 'dyn_nept_div_cur_init :'
383      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~'
384      IF(lwp) WRITE(numout,*) 'horizontal velocity divergence and'
385      IF(lwp) WRITE(numout,*) 'relative vorticity of Neptune flow'
386#if defined key_noslip_accurate
387   !!----------------------------------------------------------------------
388   !!   'key_noslip_accurate'                     2nd order centered scheme
389   !!                                                4th order at the coast
390   !!----------------------------------------------------------------------
391      IF(lwp) WRITE(numout,*)
392      IF(lwp) WRITE(numout,*) 'WARNING: key_noslip_accurate option'
393      IF(lwp) WRITE(numout,*) 'not implemented in simplified Neptune'
394      CALL ctl_warn( ' noslip_accurate option not implemented' )
395#endif
396
397   !!----------------------------------------------------------------------
398   !!   Default option                           2nd order centered schemes
399   !!----------------------------------------------------------------------
400
401      ! Apply the div and curl operators to the depth-dependent velocity
402      ! field produced by multiplying (zunep, zvnep) by (umask, vmask), exactly
403      ! equivalent to the equivalent calculation in the unsimplified code
404      !                                                ! ===============
405      DO jk = 1, jpkm1                                 ! Horizontal slab
406         !                                             ! ===============
407         !                                             ! --------
408         ! Horizontal divergence                       !   div
409         !                                             ! --------
410         DO jj = 2, jpjm1
411            DO ji = fs_2, fs_jpim1   ! vector opt.
412               zhdivnep(ji,jj,jk) =   &
413               &   (  e2u(ji  ,jj  )*fse3u(ji  ,jj  ,jk) * zunep(ji  ,jj  ) * umask(ji  ,jj  ,jk)    &
414               &    - e2u(ji-1,jj  )*fse3u(ji-1,jj  ,jk) * zunep(ji-1,jj  ) * umask(ji-1,jj  ,jk)    &
415               &    + e1v(ji  ,jj  )*fse3v(ji  ,jj  ,jk) * zvnep(ji  ,jj  ) * vmask(ji  ,jj  ,jk)    &
416               &    - e1v(ji  ,jj-1)*fse3v(ji  ,jj-1,jk) * zvnep(ji  ,jj-1) * vmask(ji  ,jj-1,jk) )  &
417               &  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
418            END DO
419         END DO
420
421         IF( .NOT. AGRIF_Root() ) THEN
422            IF ((nbondi ==  1).OR.(nbondi == 2))  zhdivnep(nlci-1 , :     ,jk) = 0.0_wp   ! east
423            IF ((nbondi == -1).OR.(nbondi == 2))  zhdivnep(2      , :     ,jk) = 0.0_wp   ! west
424            IF ((nbondj ==  1).OR.(nbondj == 2))  zhdivnep(:      ,nlcj-1 ,jk) = 0.0_wp   ! north
425            IF ((nbondj == -1).OR.(nbondj == 2))  zhdivnep(:      ,2      ,jk) = 0.0_wp   ! south
426         ENDIF
427
428         !                                             ! --------
429         ! relative vorticity                          !   rot
430         !                                             ! --------
431         DO jj = 1, jpjm1
432            DO ji = 1, fs_jpim1   ! vector opt.
433               zmrotnep(ji,jj,jk) =   &
434                  &       (  e2v(ji+1,jj  ) * zvnep(ji+1,jj  ) * vmask(ji+1,jj  ,jk)     &
435                  &        - e2v(ji  ,jj  ) * zvnep(ji  ,jj  ) * vmask(ji  ,jj  ,jk)     &
436                  &        - e1u(ji  ,jj+1) * zunep(ji  ,jj+1) * umask(ji  ,jj+1,jk)     &
437                  &        + e1u(ji  ,jj  ) * zunep(ji  ,jj  ) * umask(ji  ,jj  ,jk)  )  &
438                  &       * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
439            END DO
440         END DO
441         !                                             ! ===============
442      END DO                                           !   End of slab
443      !                                                ! ===============
444
445      ! 4. Lateral boundary conditions on zhdivnep and zmrotnep
446      ! ----------------------------------=======-----=======
447      CALL lbc_lnk( zhdivnep, 'T', 1. )   ;   CALL lbc_lnk( zmrotnep , 'F', 1. )     ! lateral boundary cond. (no sign change)
448      !
449   END SUBROUTINE dyn_nept_div_cur_init
450
451
452   SUBROUTINE dyn_nept_cor( kt )
453      !!----------------------------------------------------------------------
454      !!                  ***  ROUTINE dyn_nept_cor  ***
455      !!
456      !! ** Purpose :  Add or subtract the Neptune velocity from the now velocities
457      !!
458      !! ** Method  :  First call : kt not equal to lastkt -> subtract zunep, zvnep
459      !!               Second call: kt     equal to lastkt -> add zunep, zvnep
460      !!----------------------------------------------------------------------
461      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
462      !!
463      INTEGER, SAVE :: lastkt             ! store previous kt
464      DATA lastkt/-1/                     ! initialise previous kt
465      !!----------------------------------------------------------------------
466      !
467      IF( ln_neptsimp ) THEN
468         !
469         IF( lastkt /= kt ) THEN          ! 1st call for this kt: subtract the Neptune velocities zunep, zvnep from un, vn
470            CALL dyn_nept_vel( -1 )      ! -1 = subtract
471            !
472         ELSE                              ! 2nd call for this kt: add the Neptune velocities zunep, zvnep to un, vn
473            CALL dyn_nept_vel(  1 )      !  1 = add
474            !
475         ENDIF
476         !
477         lastkt = kt     ! Store kt
478         !
479      ENDIF
480      !
481   END SUBROUTINE dyn_nept_cor
482
483
484   SUBROUTINE dyn_nept_vel( ksign )
485      !!----------------------------------------------------------------------
486      !!                  ***  ROUTINE dyn_nept_vel  ***
487      !!
488      !! ** Purpose :  Add or subtract the Neptune velocity from the now
489      !!               velocities based on ksign
490      !!----------------------------------------------------------------------
491      INTEGER, INTENT( in ) ::   ksign    ! 1 or -1 to add or subtract neptune velocities
492      !!
493      INTEGER :: jk                       ! dummy loop index
494      !!----------------------------------------------------------------------
495      !
496      ! Adjust the current velocity un, vn by adding or subtracting the
497      ! Neptune velocities zunep, zvnep, as determined by argument ksign
498      DO jk=1, jpk
499         un(:,:,jk) = un(:,:,jk) + ksign * zunep(:,:) * umask(:,:,jk)
500         vn(:,:,jk) = vn(:,:,jk) + ksign * zvnep(:,:) * vmask(:,:,jk)
501      END DO
502      !
503   END SUBROUTINE dyn_nept_vel
504
505
506   SUBROUTINE dyn_nept_smooth_vel( htold, htnew, ld_option )
507
508      !!----------------------------------------------------------------------
509      !!                  ***  ROUTINE dyn_nept_smooth_vel  ***
510      !!
511      !! ** Purpose :  Compute smoothed topography field.
512      !!
513      !! ** Action : - Updates the array htnew (output) with a smoothed
514      !!               version of the (input) array htold. Form of smoothing
515      !!               algorithm is controlled by the (logical) argument ld_option.
516      !!----------------------------------------------------------------------
517      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )  ::  htold      ! temporary 2D workspace
518      LOGICAL                     , INTENT(in   )  ::  ld_option  ! temporary 2D workspace
519      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)  ::  htnew      ! temporary 2D workspace
520      !
521      INTEGER                           ::  ji, jj  ! dummy loop indices
522      INTEGER , POINTER, DIMENSION(:,:) ::  nb, iwork
523      REAL(wp), POINTER, DIMENSION(:,:) ::  work    ! temporary 2D workspace
524      !!----------------------------------------------------------------------
525      !
526      CALL wrk_alloc( jpi, jpj, nb, iwork ) 
527      CALL wrk_alloc( jpi, jpj, work      ) 
528      !
529      iwork(:,:) = 0
530
531      !! iwork is a mask of gridpoints: iwork = 1 => ocean, iwork = 0 => land
532      WHERE( htold(:,:) > 0 )
533         iwork(:,:) = 1
534         htnew(:,:) = htold(:,:)
535      ELSEWHERE
536         iwork(:,:) = 0
537         htnew(:,:) = 0.0_wp
538      END WHERE
539      !! htnew contains valid ocean depths from htold, or zero
540
541      !! set work to a smoothed/averaged version of htnew; choice controlled by ld_option
542      !! nb is set to the sum of the weights of the valid values used in work
543      IF( ld_option ) THEN
544
545         !! Apply scale-selective smoothing in determining work from htnew
546         DO jj=2,jpj-1
547            DO ji=2,jpi-1
548               work(ji,jj) = 4.0*htnew( ji , jj ) +                        &
549                           & 2.0*htnew(ji+1, jj ) + 2.0*htnew(ji-1, jj ) + &
550                           & 2.0*htnew( ji ,jj+1) + 2.0*htnew( ji ,jj-1) + &
551                           &     htnew(ji+1,jj+1) +     htnew(ji+1,jj-1) + &
552                           &     htnew(ji-1,jj+1) +     htnew(ji-1,jj-1)
553
554               nb(ji,jj)   = 4 * iwork( ji , jj ) +                        &
555                           & 2 * iwork(ji+1, jj ) + 2 * iwork(ji-1, jj ) + &
556                           & 2 * iwork( ji ,jj+1) + 2 * iwork( ji ,jj-1) + &
557                           &     iwork(ji+1,jj+1) +     iwork(ji+1,jj-1) + &
558                           &     iwork(ji-1,jj+1) +     iwork(ji-1,jj-1)
559            END DO
560         END DO
561
562      ELSE
563
564         !! Apply simple 9-point averaging in determining work from htnew
565         DO jj=2,jpj-1
566            DO ji=2,jpi-1
567               work(ji,jj) =  htnew( ji , jj ) +                    &
568                           &  htnew(ji+1, jj ) + htnew(ji-1, jj ) + &
569                           &  htnew( ji ,jj+1) + htnew( ji ,jj-1) + &
570                           &  htnew(ji+1,jj+1) + htnew(ji+1,jj-1) + &
571                           &  htnew(ji-1,jj+1) + htnew(ji-1,jj-1)
572
573               nb(ji,jj) =    iwork( ji , jj ) +                    &
574                           &  iwork(ji+1, jj ) + iwork(ji-1, jj ) + &
575                           &  iwork( ji ,jj+1) + iwork( ji ,jj-1) + &
576                           &  iwork(ji+1,jj+1) + iwork(ji+1,jj-1) + &
577                           &  iwork(ji-1,jj+1) + iwork(ji-1,jj-1)
578            END DO
579         END DO
580
581      ENDIF
582
583      !! write averaged value of work into htnew,
584      !! if average is valid and point is unmasked
585      WHERE( (htold(:,:) /= 0.0_wp ) .AND. ( nb(:,:) /= 0 ) )
586         htnew(:,:) = work(:,:)/real(nb(:,:))
587      ELSEWHERE
588         htnew(:,:) = 0.0_wp
589      END WHERE
590
591      !!    Deallocate temporary workspace arrays, all local to this routine
592      CALL wrk_dealloc( jpi, jpj, nb, iwork ) 
593      CALL wrk_dealloc( jpi, jpj, work      ) 
594      !
595   END SUBROUTINE dyn_nept_smooth_vel
596
597END MODULE dynnept
Note: See TracBrowser for help on using the repository browser.