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/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90 @ 4147

Last change on this file since 4147 was 4147, checked in by cetlod, 10 years ago

merge in dev_LOCEAN_2013, the 1st development branch dev_r3853_CNRS9_Confsetting, from its starting point ( r3853 ) on the trunk: see ticket #1169
File size: 27.8 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 obc_oce          ! ocean lateral open boundary condition
27   USE in_out_manager   ! I/O manager
28   USE lib_mpp          ! distributed memory computing
29   USE prtctl           ! Print control
30   USE phycst
31   USE lbclnk
32   USE wrk_nemo        ! Memory Allocation
33
34   IMPLICIT NONE
35   PRIVATE
36
37   !! * Routine accessibility
38   PUBLIC dyn_nept_init      ! routine called by nemogcm.F90
39   PUBLIC dyn_nept_cor       ! routine called by step.F90
40   !! dynnept_alloc()         is called only by dyn_nept_init, within this module
41   !! dyn_nept_div_cur_init   is called only by dyn_nept_init, within this module
42   !! dyn_nept_vel            is called only by dyn_nept_cor,  within this module
43
44   !! * Shared module variables
45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: zunep, zvnep  ! Neptune u and v
46   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: zhdivnep      ! hor. div for Neptune vel.
47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: zmrotnep      ! curl for Neptune vel.
48
49
50   !! * Namelist namdyn_nept variables
51   LOGICAL, PUBLIC  ::  ln_neptsimp          ! yes/no simplified neptune
52
53   LOGICAL          ::  ln_smooth_neptvel    ! yes/no smooth zunep, zvnep
54   REAL(wp)         ::  rn_tslse             ! value of lengthscale L at the equator
55   REAL(wp)         ::  rn_tslsp             ! value of lengthscale L at the pole
56!! Specify whether to ramp down the Neptune velocity in shallow
57!! water, and the depth range controlling such ramping down
58   LOGICAL          ::  ln_neptramp          ! ramp down Neptune velocity in shallow water
59   REAL(wp)         ::  rn_htrmin            ! min. depth of transition range
60   REAL(wp)         ::  rn_htrmax            ! max. depth of transition range
61
62   !! * Module variables
63
64
65   !! * Substitutions
66#  include "vectopt_loop_substitute.h90"
67#  include "domzgr_substitute.h90"
68   !!----------------------------------------------------------------------
69   !!   OPA 9.0 , implemented by Bedford Institute of Oceanography
70   !!----------------------------------------------------------------------
71
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(:,:  ) :: ht, 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      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     , ht, 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, ht, .TRUE. )
203      !! overwrites ht with a smoothed version of htn
204      ELSE
205         ht(:,:) = htn(:,:)
206      !! use unsmoothed version of htn
207      ENDIF
208      CALL lbc_lnk( ht, '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(:,:) * ht(:,:)
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     , ht, 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 defined key_obc
422         IF( Agrif_Root() ) THEN
423            ! open boundaries (div must be zero behind the open boundary)
424            !  mpp remark: The zeroing of zhdivnep can probably be extended to 1->jpi/jpj for the correct row/column
425            IF( lp_obc_east  )  zhdivnep(nie0p1:nie1p1,nje0  :nje1  ,jk) = 0.0_wp      ! east
426            IF( lp_obc_west  )  zhdivnep(niw0  :niw1  ,njw0  :njw1  ,jk) = 0.0_wp      ! west
427            IF( lp_obc_north )  zhdivnep(nin0  :nin1  ,njn0p1:njn1p1,jk) = 0.0_wp      ! north
428            IF( lp_obc_south )  zhdivnep(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.0_wp      ! south
429         ENDIF
430#endif
431         IF( .NOT. AGRIF_Root() ) THEN
432            IF ((nbondi ==  1).OR.(nbondi == 2))  zhdivnep(nlci-1 , :     ,jk) = 0.0_wp   ! east
433            IF ((nbondi == -1).OR.(nbondi == 2))  zhdivnep(2      , :     ,jk) = 0.0_wp   ! west
434            IF ((nbondj ==  1).OR.(nbondj == 2))  zhdivnep(:      ,nlcj-1 ,jk) = 0.0_wp   ! north
435            IF ((nbondj == -1).OR.(nbondj == 2))  zhdivnep(:      ,2      ,jk) = 0.0_wp   ! south
436         ENDIF
437
438         !                                             ! --------
439         ! relative vorticity                          !   rot
440         !                                             ! --------
441         DO jj = 1, jpjm1
442            DO ji = 1, fs_jpim1   ! vector opt.
443               zmrotnep(ji,jj,jk) =   &
444                  &       (  e2v(ji+1,jj  ) * zvnep(ji+1,jj  ) * vmask(ji+1,jj  ,jk)     &
445                  &        - e2v(ji  ,jj  ) * zvnep(ji  ,jj  ) * vmask(ji  ,jj  ,jk)     &
446                  &        - e1u(ji  ,jj+1) * zunep(ji  ,jj+1) * umask(ji  ,jj+1,jk)     &
447                  &        + e1u(ji  ,jj  ) * zunep(ji  ,jj  ) * umask(ji  ,jj  ,jk)  )  &
448                  &       * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
449            END DO
450         END DO
451         !                                             ! ===============
452      END DO                                           !   End of slab
453      !                                                ! ===============
454
455      ! 4. Lateral boundary conditions on zhdivnep and zmrotnep
456      ! ----------------------------------=======-----=======
457      CALL lbc_lnk( zhdivnep, 'T', 1. )   ;   CALL lbc_lnk( zmrotnep , 'F', 1. )     ! lateral boundary cond. (no sign change)
458      !
459   END SUBROUTINE dyn_nept_div_cur_init
460
461
462   SUBROUTINE dyn_nept_cor( kt )
463      !!----------------------------------------------------------------------
464      !!                  ***  ROUTINE dyn_nept_cor  ***
465      !!
466      !! ** Purpose :  Add or subtract the Neptune velocity from the now velocities
467      !!
468      !! ** Method  :  First call : kt not equal to lastkt -> subtract zunep, zvnep
469      !!               Second call: kt     equal to lastkt -> add zunep, zvnep
470      !!----------------------------------------------------------------------
471      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
472      !!
473      INTEGER, SAVE :: lastkt             ! store previous kt
474      DATA lastkt/-1/                     ! initialise previous kt
475      !!----------------------------------------------------------------------
476      !
477      IF( ln_neptsimp ) THEN
478         !
479         IF( lastkt /= kt ) THEN          ! 1st call for this kt: subtract the Neptune velocities zunep, zvnep from un, vn
480            CALL dyn_nept_vel( -1 )      ! -1 = subtract
481            !
482         ELSE                              ! 2nd call for this kt: add the Neptune velocities zunep, zvnep to un, vn
483            CALL dyn_nept_vel(  1 )      !  1 = add
484            !
485         ENDIF
486         !
487         lastkt = kt     ! Store kt
488         !
489      ENDIF
490      !
491   END SUBROUTINE dyn_nept_cor
492
493
494   SUBROUTINE dyn_nept_vel( ksign )
495      !!----------------------------------------------------------------------
496      !!                  ***  ROUTINE dyn_nept_vel  ***
497      !!
498      !! ** Purpose :  Add or subtract the Neptune velocity from the now
499      !!               velocities based on ksign
500      !!----------------------------------------------------------------------
501      INTEGER, INTENT( in ) ::   ksign    ! 1 or -1 to add or subtract neptune velocities
502      !!
503      INTEGER :: jk                       ! dummy loop index
504      !!----------------------------------------------------------------------
505      !
506      ! Adjust the current velocity un, vn by adding or subtracting the
507      ! Neptune velocities zunep, zvnep, as determined by argument ksign
508      DO jk=1, jpk
509         un(:,:,jk) = un(:,:,jk) + ksign * zunep(:,:) * umask(:,:,jk)
510         vn(:,:,jk) = vn(:,:,jk) + ksign * zvnep(:,:) * vmask(:,:,jk)
511      END DO
512      !
513   END SUBROUTINE dyn_nept_vel
514
515
516   SUBROUTINE dyn_nept_smooth_vel( htold, htnew, ld_option )
517
518      !!----------------------------------------------------------------------
519      !!                  ***  ROUTINE dyn_nept_smooth_vel  ***
520      !!
521      !! ** Purpose :  Compute smoothed topography field.
522      !!
523      !! ** Action : - Updates the array htnew (output) with a smoothed
524      !!               version of the (input) array htold. Form of smoothing
525      !!               algorithm is controlled by the (logical) argument ld_option.
526      !!----------------------------------------------------------------------
527      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )  ::  htold      ! temporary 2D workspace
528      LOGICAL                     , INTENT(in   )  ::  ld_option  ! temporary 2D workspace
529      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)  ::  htnew      ! temporary 2D workspace
530      !
531      INTEGER                           ::  ji, jj  ! dummy loop indices
532      INTEGER , POINTER, DIMENSION(:,:) ::  nb, iwork
533      REAL(wp), POINTER, DIMENSION(:,:) ::  work    ! temporary 2D workspace
534      !!----------------------------------------------------------------------
535      !
536      CALL wrk_alloc( jpi, jpj, nb, iwork ) 
537      CALL wrk_alloc( jpi, jpj, work      ) 
538      !
539      iwork(:,:) = 0
540
541      !! iwork is a mask of gridpoints: iwork = 1 => ocean, iwork = 0 => land
542      WHERE( htold(:,:) > 0 )
543         iwork(:,:) = 1
544         htnew(:,:) = htold(:,:)
545      ELSEWHERE
546         iwork(:,:) = 0
547         htnew(:,:) = 0.0_wp
548      END WHERE
549      !! htnew contains valid ocean depths from htold, or zero
550
551      !! set work to a smoothed/averaged version of htnew; choice controlled by ld_option
552      !! nb is set to the sum of the weights of the valid values used in work
553      IF( ld_option ) THEN
554
555         !! Apply scale-selective smoothing in determining work from htnew
556         DO jj=2,jpj-1
557            DO ji=2,jpi-1
558               work(ji,jj) = 4.0*htnew( ji , jj ) +                        &
559                           & 2.0*htnew(ji+1, jj ) + 2.0*htnew(ji-1, jj ) + &
560                           & 2.0*htnew( ji ,jj+1) + 2.0*htnew( ji ,jj-1) + &
561                           &     htnew(ji+1,jj+1) +     htnew(ji+1,jj-1) + &
562                           &     htnew(ji-1,jj+1) +     htnew(ji-1,jj-1)
563
564               nb(ji,jj)   = 4 * iwork( ji , jj ) +                        &
565                           & 2 * iwork(ji+1, jj ) + 2 * iwork(ji-1, jj ) + &
566                           & 2 * iwork( ji ,jj+1) + 2 * iwork( ji ,jj-1) + &
567                           &     iwork(ji+1,jj+1) +     iwork(ji+1,jj-1) + &
568                           &     iwork(ji-1,jj+1) +     iwork(ji-1,jj-1)
569            END DO
570         END DO
571
572      ELSE
573
574         !! Apply simple 9-point averaging in determining work from htnew
575         DO jj=2,jpj-1
576            DO ji=2,jpi-1
577               work(ji,jj) =  htnew( ji , jj ) +                    &
578                           &  htnew(ji+1, jj ) + htnew(ji-1, jj ) + &
579                           &  htnew( ji ,jj+1) + htnew( ji ,jj-1) + &
580                           &  htnew(ji+1,jj+1) + htnew(ji+1,jj-1) + &
581                           &  htnew(ji-1,jj+1) + htnew(ji-1,jj-1)
582
583               nb(ji,jj) =    iwork( ji , jj ) +                    &
584                           &  iwork(ji+1, jj ) + iwork(ji-1, jj ) + &
585                           &  iwork( ji ,jj+1) + iwork( ji ,jj-1) + &
586                           &  iwork(ji+1,jj+1) + iwork(ji+1,jj-1) + &
587                           &  iwork(ji-1,jj+1) + iwork(ji-1,jj-1)
588            END DO
589         END DO
590
591      ENDIF
592
593      !! write averaged value of work into htnew,
594      !! if average is valid and point is unmasked
595      WHERE( (htold(:,:) /= 0.0_wp ) .AND. ( nb(:,:) /= 0 ) )
596         htnew(:,:) = work(:,:)/real(nb(:,:))
597      ELSEWHERE
598         htnew(:,:) = 0.0_wp
599      END WHERE
600
601      !!    Deallocate temporary workspace arrays, all local to this routine
602      CALL wrk_dealloc( jpi, jpj, nb, iwork ) 
603      CALL wrk_dealloc( jpi, jpj, work      ) 
604      !
605   END SUBROUTINE dyn_nept_smooth_vel
606
607END MODULE dynnept
Note: See TracBrowser for help on using the repository browser.