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/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90 @ 3027

Last change on this file since 3027 was 3027, checked in by acc, 13 years ago

Branch dev_NOC_2011_MERGE. #874. Step 8: minor corrections for single processor running as uncovered during SETTE testing
File size: 28.2 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
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        = .FALSE.  ! yes/no simplified neptune
52
53   LOGICAL          ::  ln_smooth_neptvel  = .FALSE.  ! yes/no smooth zunep, zvnep
54   REAL(wp)         ::  rn_tslse           =  1.2e4   ! value of lengthscale L at the equator
55   REAL(wp)         ::  rn_tslsp           =  3.0e3   ! 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        = .FALSE.  ! ramp down Neptune velocity in shallow water
59   REAL(wp)         ::  rn_htrmin          =  100.0   ! min. depth of transition range
60   REAL(wp)         ::  rn_htrmax          =  200.0   ! 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!!    Local dynamic allocation
105      INTEGER  ::   ji, jj, jk    ! dummy loop indices
106      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   ht       ! temporary 2D workspace
107      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   htn      ! temporary 2D workspace
108      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   tscale   ! temporary 2D workspace
109      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   tsp      ! temporary 2D workspace
110      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   hur_n    ! temporary 2D workspace
111      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   hvr_n    ! temporary 2D workspace
112      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   hu_n     ! temporary 2D workspace
113      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   hv_n     ! temporary 2D workspace
114      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   znmask   ! temporary 3D array for nmask
115      REAL(wp) :: unemin,unemax,vnemin,vnemax   ! extrema of (u*, v*) fields
116      REAL(wp) :: zhdivmin,zhdivmax             ! extrema of horizontal divergence of (u*, v*) fields
117      REAL(wp) :: zmrotmin,zmrotmax             ! extrema of the curl of the (u*, v*) fields
118      REAL(wp) :: ustar,vstar                   ! (u*, v*) before tapering in shallow water
119      REAL(wp) :: hramp                         ! depth over which Neptune vel. is ramped down
120      !!
121      NAMELIST/namdyn_nept/ ln_neptsimp,      &
122                            ln_smooth_neptvel,&
123             rn_tslse,        &
124             rn_tslsp,        &
125                            ln_neptramp,      &
126                            rn_htrmin,        &
127             rn_htrmax
128      !!----------------------------------------------------------------------
129
130      ! Define the (simplified) Neptune parameters
131      ! ==========================================
132
133!!    WRITE(numout,*) ' start dynnept namelist'
134!!    CALL FLUSH(numout)
135      REWIND( numnam )                  ! Read Namelist namdyn_nept:  Simplified Neptune
136      READ  ( numnam, namdyn_nept )
137!!    WRITE(numout,*) ' dynnept namelist done'
138!!    CALL FLUSH(numout)
139
140      IF(lwp) THEN                      ! Control print
141         WRITE(numout,*)
142         WRITE(numout,*) 'dyn_nept_init : Simplified Neptune module enabled'
143         WRITE(numout,*) '~~~~~~~~~~~~~'
144         WRITE(numout,*) ' -->   Reading namelist namdyn_nept parameters:'
145         WRITE(numout,*) '       ln_neptsimp          = ', ln_neptsimp
146         WRITE(numout,*)
147         WRITE(numout,*) '       ln_smooth_neptvel    = ', ln_smooth_neptvel
148         WRITE(numout,*) '       rn_tslse             = ', rn_tslse
149         WRITE(numout,*) '       rn_tslsp             = ', rn_tslsp
150         WRITE(numout,*)
151         WRITE(numout,*) '       ln_neptramp          = ', ln_neptramp
152         WRITE(numout,*) '       rn_htrmin            = ', rn_htrmin
153         WRITE(numout,*) '       rn_htrmax            = ', rn_htrmax
154         WRITE(numout,*)
155         CALL FLUSH(numout)
156      ENDIF
157
158      IF( ln_smooth_neptvel ) THEN
159         IF(lwp) WRITE(numout,*) ' -->   neptune velocities will be smoothed'
160      ELSE
161         IF(lwp) WRITE(numout,*) ' -->   neptune velocities will not be smoothed'
162      ENDIF
163
164      IF( ln_neptsimp ) THEN
165          IF(lwp) WRITE(numout,*) ' -->   ln_neptsimp enabled, solving for U-UN'
166      ELSE
167          IF(lwp) WRITE(numout,*) ' -->   ln_neptsimp disabled'
168          RETURN
169      ENDIF
170
171      IF( ln_neptramp ) THEN
172          IF(lwp) WRITE(numout,*) ' -->   ln_neptramp enabled, ramp down Neptune'
173          IF(lwp) WRITE(numout,*) ' -->   velocity components in shallow water'
174      ELSE
175          IF(lwp) WRITE(numout,*) ' -->   ln_neptramp disabled'
176      ENDIF
177
178
179!!    Perform dynamic allocation of shared module variables
180      IF( dynnept_alloc() /= 0 )   CALL ctl_warn('dynnept_alloc: array allocate failed.')
181
182!!    Dynamically allocate local work arrays
183      ALLOCATE( ht(jpi,jpj), htn(jpi,jpj), tscale(jpi,jpj), tsp(jpi,jpj),      &
184         &      hur_n(jpi,jpj), hvr_n(jpi,jpj), hu_n(jpi,jpj), hv_n(jpi,jpj),  &
185         &      znmask(jpi,jpj,jpk) )
186
187      IF( .not. ln_rstart ) THEN      ! If restarting, these arrays are read from the restart file
188         zhdivnep(:,:,:) = 0.0_wp
189         zmrotnep(:,:,:) = 0.0_wp
190      END IF
191
192      ! Computation of nmask: same as fmask, but fmask cannot be used
193      ! because it is modified after it is computed in dom_msk
194      ! (this can be optimised to save memory, such as merge into next loop)
195      DO jk = 1, jpk
196         DO jj = 1, jpjm1
197            DO ji = 1, fs_jpim1   ! vector loop
198               znmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
199                   &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
200            END DO
201         END DO
202      END DO
203
204      CALL lbc_lnk( znmask, 'F', 1.0_wp )
205
206
207      ! now compute zunep, zvnep (renamed from earlier versions)
208
209      zunep(:,:) = 0.0_wp
210      zvnep(:,:) = 0.0_wp
211
212      htn(:,:) = 0.0_wp            ! ocean depth at F-point
213      DO jk = 1, jpk
214         htn(:,:) = htn(:,:) + fse3f(:,:,jk) * znmask(:,:,jk)
215      END DO
216
217      IF( ln_smooth_neptvel ) THEN
218         CALL dyn_nept_smooth_vel( htn, ht, .TRUE. )
219      !! overwrites ht with a smoothed version of htn
220      ELSE
221         ht(:,:) = htn(:,:)
222      !! use unsmoothed version of htn
223      ENDIF
224      CALL lbc_lnk( ht, 'F', 1.0_wp )
225
226      !! Compute tsp, a stream function for the Neptune velocity,
227      !! with the usual geophysical sign convention
228      !! Then zunep = -latitudinal derivative "-(1/H)*d(tsp)/dy"
229      !!      zvnep = longitudinal derivative " (1/H)*d(tsp)/dx"
230
231      tsp(:,:)    = 0.0_wp
232      tscale(:,:) = 0.0_wp
233
234      tscale(:,:) = rn_tslsp + (rn_tslse - rn_tslsp) *   &
235                   ( 0.5_wp + 0.5_wp * COS( 2.0_wp * rad * gphif(:,:) )  )
236      tsp   (:,:) = -2.0_wp * omega * SIN( rad * gphif(:,:) ) * tscale(:,:) * tscale(:,:) * ht(:,:)
237
238
239      IF( ln_smooth_neptvel ) THEN
240         CALL dyn_nept_smooth_vel( hu, hu_n, .TRUE. )
241      !! overwrites hu_n with a smoothed version of hu
242      ELSE
243         hu_n(:,:) = hu(:,:)
244      !! use unsmoothed version of hu
245      ENDIF
246      CALL lbc_lnk( hu_n, 'U', 1.0_wp )
247      hu_n(:,:) = hu_n(:,:) * umask(:,:,1)
248
249      WHERE( hu_n(:,:) == 0.0_wp )
250         hur_n(:,:) = 0.0_wp
251      ELSEWHERE
252         hur_n(:,:) = 1.0_wp / hu_n(:,:)
253      END WHERE
254
255
256      IF( ln_smooth_neptvel ) THEN
257         CALL dyn_nept_smooth_vel( hv, hv_n, .TRUE. )
258      !! overwrites hv_n with a smoothed version of hv
259      ELSE
260         hv_n(:,:) = hv(:,:)
261      !! use unsmoothed version of hv
262      ENDIF
263      CALL lbc_lnk( hv_n, 'V', 1.0_wp )
264      hv_n(:,:) = hv_n(:,:) * vmask(:,:,1)
265
266      WHERE( hv_n == 0.0_wp )
267         hvr_n(:,:) = 0.0_wp
268      ELSEWHERE
269         hvr_n(:,:) = 1.0_wp / hv_n(:,:)
270      END WHERE
271
272
273      unemin =  1.0e35
274      unemax = -1.0e35
275      vnemin =  1.0e35
276      vnemax = -1.0e35
277      hramp = rn_htrmax - rn_htrmin
278      DO jj = 2, jpj-1
279         DO ji = 2, jpi-1
280            if ( umask(ji,jj,1) /= 0.0_wp ) then
281               ustar =-1.0_wp/e2u(ji,jj) * hur_n(ji,jj) * ( tsp(ji,jj)-tsp(ji,jj-1) ) * umask(ji,jj,1)
282               if ( ln_neptramp ) then
283!!                Apply ramp down to velocity component
284                  if ( hu_n(ji,jj) <= rn_htrmin ) then
285                    zunep(ji,jj) = 0.0_wp
286                   else if ( hu_n(ji,jj) >= rn_htrmax ) then
287                    zunep(ji,jj) = ustar
288                   else if ( hramp > 0.0_wp ) then
289                    zunep(ji,jj) = ( hu_n(ji,jj) - rn_htrmin) * ustar/hramp
290                  endif
291                else
292                 zunep(ji,jj) = ustar
293               endif
294             else
295              zunep(ji,jj) = 0.0_wp
296            endif
297            if ( vmask(ji,jj,1) /= 0.0_wp ) then
298               vstar = 1.0_wp/e1v(ji,jj) * hvr_n(ji,jj) * ( tsp(ji,jj)-tsp(ji-1,jj) ) * vmask(ji,jj,1)
299               if ( ln_neptramp ) then
300!!                Apply ramp down to velocity component
301                  if ( hv_n(ji,jj) <= rn_htrmin ) then
302                    zvnep(ji,jj) = 0.0_wp
303                   else if ( hv_n(ji,jj) >= rn_htrmax ) then
304                    zvnep(ji,jj) = vstar
305                   else if ( hramp > 0.0_wp ) then
306                    zvnep(ji,jj) = ( hv_n(ji,jj) - rn_htrmin) * vstar/hramp
307                  endif
308                else
309                  zvnep(ji,jj) = vstar
310               endif
311             else
312              zvnep(ji,jj) = 0.0_wp
313            endif
314            unemin = min( unemin, zunep(ji,jj) )
315            unemax = max( unemax, zunep(ji,jj) )
316            vnemin = min( vnemin, zvnep(ji,jj) )
317            vnemax = max( vnemax, zvnep(ji,jj) )
318         END DO
319      END DO
320      CALL lbc_lnk( zunep, 'U', -1.0_wp )
321      CALL lbc_lnk( zvnep, 'V', -1.0_wp )
322      WRITE(numout,*) '      zunep: min, max       = ', unemin,unemax
323      WRITE(numout,*) '      zvnep: min, max       = ', vnemin,vnemax
324      WRITE(numout,*)
325
326      !!  Compute, once and for all, the horizontal divergence (zhdivnep)
327      !!  and the curl (zmrotnep) of the Neptune velocity field (zunep, zvnep)
328      CALL dyn_nept_div_cur_init
329
330      !! Check the ranges of the computed divergence & vorticity
331      zhdivmin =  1.0e35
332      zhdivmax = -1.0e35
333      zmrotmin =  1.0e35
334      zmrotmax = -1.0e35
335      hramp = rn_htrmax - rn_htrmin
336      DO jk = 1, jpkm1                                 ! Horizontal slab
337         DO jj = 2, jpj-1
338            DO ji = 2, jpi-1
339               zhdivmin = min( zhdivmin, zhdivnep(ji,jj,jk) )
340               zhdivmax = max( zhdivmax, zhdivnep(ji,jj,jk) )
341               zmrotmin = min( zmrotmin, zmrotnep(ji,jj,jk) )
342               zmrotmax = max( zmrotmax, zmrotnep(ji,jj,jk) )
343            END DO
344         END DO
345      END DO
346      WRITE(numout,*) '   zhdivnep: min, max       = ', zhdivmin,zhdivmax
347      WRITE(numout,*) '   zmrotnep: min, max       = ', zmrotmin,zmrotmax
348      WRITE(numout,*)
349
350!!    Deallocate temporary workspace arrays, which are all local to
351!!    this routine, except where passed as arguments to other routines
352      DEALLOCATE( ht, htn, tscale, tsp, hur_n, hvr_n, hu_n, hv_n, znmask )
353
354   END SUBROUTINE dyn_nept_init
355
356
357   SUBROUTINE dyn_nept_div_cur_init
358      !!----------------------------------------------------------------------
359      !!             ***  ROUTINE dyn_nept_div_cur_init  ***
360      !!
361      !! ** Purpose :   compute the horizontal divergence and the relative
362      !!                vorticity of the time-invariant u* and v* Neptune
363      !!                effect velocities (called zunep, zvnep)
364      !!
365      !! ** Method  : - Divergence:
366      !!      - compute the divergence given by :
367      !!         zhdivnep = 1/(e1t*e2t*e3t) ( di[e2u*e3u zunep] + dj[e1v*e3v zvnep] )
368      !!      - compute the curl in tensorial formalism:
369      !!         zmrotnep = 1/(e1f*e2f) ( di[e2v zvnep] - dj[e1u zunep] )
370      !!      Note: Coastal boundary condition: lateral friction set through
371      !!      the value of fmask along the coast (see dommsk.F90) and shlat
372      !!      (namelist parameter)
373      !!
374      !! ** Action  : - compute zhdivnep, the hor. divergence of (u*, v*)
375      !!              - compute zmrotnep, the rel. vorticity  of (u*, v*)
376      !!
377      !! History :  OPA  ! 1987-06  (P. Andrich, D. L Hostis)  Original code
378      !!            4.0  ! 1991-11  (G. Madec)
379      !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions
380      !!            7.0  ! 1996-01  (G. Madec)  s-coordinates
381      !!            8.0  ! 1997-06  (G. Madec)  lateral boundary cond., lbc
382      !!            8.1  ! 1997-08  (J.M. Molines)  Open boundaries
383      !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
384      !!  NEMO      1.0  ! 2002-09  (G. Madec, E. Durand)  Free form, F90
385      !!             -   ! 2005-01  (J. Chanut) Unstructured open boundaries
386      !!             -   ! 2003-08  (G. Madec)  merged of cur and div, free form, F90
387      !!             -   ! 2005-01  (J. Chanut, A. Sellar) unstructured open boundaries
388      !!            3.3  ! 2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
389      !!                 ! 2011-06  (Jeff Blundell, NOCS) Adapt code from divcur.F90
390      !!                 !           to compute Neptune effect fields only
391      !!----------------------------------------------------------------------
392      !
393      INTEGER  ::   ji, jj, jk          ! dummy loop indices
394      !!----------------------------------------------------------------------
395
396
397      IF(lwp) WRITE(numout,*)
398      IF(lwp) WRITE(numout,*) 'dyn_nept_div_cur_init :'
399      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~'
400      IF(lwp) WRITE(numout,*) 'horizontal velocity divergence and'
401      IF(lwp) WRITE(numout,*) 'relative vorticity of Neptune flow'
402#if defined key_noslip_accurate
403   !!----------------------------------------------------------------------
404   !!   'key_noslip_accurate'                     2nd order centered scheme
405   !!                                                4th order at the coast
406   !!----------------------------------------------------------------------
407      IF(lwp) WRITE(numout,*)
408      IF(lwp) WRITE(numout,*) 'WARNING: key_noslip_accurate option'
409      IF(lwp) WRITE(numout,*) 'not implemented in simplified Neptune'
410      CALL ctl_warn( ' noslip_accurate option not implemented' )
411#endif
412
413   !!----------------------------------------------------------------------
414   !!   Default option                           2nd order centered schemes
415   !!----------------------------------------------------------------------
416
417      ! Apply the div and curl operators to the depth-dependent velocity
418      ! field produced by multiplying (zunep, zvnep) by (umask, vmask), exactly
419      ! equivalent to the equivalent calculation in the unsimplified code
420      !                                                ! ===============
421      DO jk = 1, jpkm1                                 ! Horizontal slab
422         !                                             ! ===============
423         !                                             ! --------
424         ! Horizontal divergence                       !   div
425         !                                             ! --------
426         DO jj = 2, jpjm1
427            DO ji = fs_2, fs_jpim1   ! vector opt.
428               zhdivnep(ji,jj,jk) =   &
429               &   (  e2u(ji  ,jj  )*fse3u(ji  ,jj  ,jk) * zunep(ji  ,jj  ) * umask(ji  ,jj  ,jk)    &
430               &    - e2u(ji-1,jj  )*fse3u(ji-1,jj  ,jk) * zunep(ji-1,jj  ) * umask(ji-1,jj  ,jk)    &
431               &    + e1v(ji  ,jj  )*fse3v(ji  ,jj  ,jk) * zvnep(ji  ,jj  ) * vmask(ji  ,jj  ,jk)    &
432               &    - e1v(ji  ,jj-1)*fse3v(ji  ,jj-1,jk) * zvnep(ji  ,jj-1) * vmask(ji  ,jj-1,jk) )  &
433               &  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
434            END DO
435         END DO
436
437#if defined key_obc
438         IF( Agrif_Root() ) THEN
439            ! open boundaries (div must be zero behind the open boundary)
440            !  mpp remark: The zeroing of zhdivnep can probably be extended to 1->jpi/jpj for the correct row/column
441            IF( lp_obc_east  )  zhdivnep(nie0p1:nie1p1,nje0  :nje1  ,jk) = 0.0_wp      ! east
442            IF( lp_obc_west  )  zhdivnep(niw0  :niw1  ,njw0  :njw1  ,jk) = 0.0_wp      ! west
443            IF( lp_obc_north )  zhdivnep(nin0  :nin1  ,njn0p1:njn1p1,jk) = 0.0_wp      ! north
444            IF( lp_obc_south )  zhdivnep(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.0_wp      ! south
445         ENDIF
446#endif
447         IF( .NOT. AGRIF_Root() ) THEN
448            IF ((nbondi ==  1).OR.(nbondi == 2))  zhdivnep(nlci-1 , :     ,jk) = 0.0_wp   ! east
449            IF ((nbondi == -1).OR.(nbondi == 2))  zhdivnep(2      , :     ,jk) = 0.0_wp   ! west
450            IF ((nbondj ==  1).OR.(nbondj == 2))  zhdivnep(:      ,nlcj-1 ,jk) = 0.0_wp   ! north
451            IF ((nbondj == -1).OR.(nbondj == 2))  zhdivnep(:      ,2      ,jk) = 0.0_wp   ! south
452         ENDIF
453
454         !                                             ! --------
455         ! relative vorticity                          !   rot
456         !                                             ! --------
457         DO jj = 1, jpjm1
458            DO ji = 1, fs_jpim1   ! vector opt.
459               zmrotnep(ji,jj,jk) =   &
460                  &       (  e2v(ji+1,jj  ) * zvnep(ji+1,jj  ) * vmask(ji+1,jj  ,jk)     &
461                  &        - e2v(ji  ,jj  ) * zvnep(ji  ,jj  ) * vmask(ji  ,jj  ,jk)     &
462                  &        - e1u(ji  ,jj+1) * zunep(ji  ,jj+1) * umask(ji  ,jj+1,jk)     &
463                  &        + e1u(ji  ,jj  ) * zunep(ji  ,jj  ) * umask(ji  ,jj  ,jk)  )  &
464                  &       * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
465            END DO
466         END DO
467         !                                             ! ===============
468      END DO                                           !   End of slab
469      !                                                ! ===============
470
471      ! 4. Lateral boundary conditions on zhdivnep and zmrotnep
472      ! ----------------------------------=======-----=======
473      CALL lbc_lnk( zhdivnep, 'T', 1. )   ;   CALL lbc_lnk( zmrotnep , 'F', 1. )     ! lateral boundary cond. (no sign change)
474      !
475   END SUBROUTINE dyn_nept_div_cur_init
476
477
478   SUBROUTINE dyn_nept_cor( kt )
479      !!----------------------------------------------------------------------
480      !!                  ***  ROUTINE dyn_nept_cor  ***
481      !!
482      !! ** Purpose :  Add or subtract the Neptune velocity from the now velocities
483      !!
484      !! ** Method  :  First call : kt not equal to lastkt -> subtract zunep, zvnep
485      !!               Second call: kt     equal to lastkt -> add zunep, zvnep
486      !!----------------------------------------------------------------------
487      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
488      !!
489      INTEGER, SAVE :: lastkt             ! store previous kt
490      DATA lastkt/-1/                     ! initialise previous kt
491      !!----------------------------------------------------------------------
492      !
493      IF( ln_neptsimp ) THEN
494         !
495         IF( lastkt /= kt ) THEN          ! 1st call for this kt: subtract the Neptune velocities zunep, zvnep from un, vn
496            CALL dyn_nept_vel( -1 )      ! -1 = subtract
497            !
498         ELSE                              ! 2nd call for this kt: add the Neptune velocities zunep, zvnep to un, vn
499            CALL dyn_nept_vel(  1 )      !  1 = add
500            !
501         ENDIF
502         !
503    lastkt = kt        ! Store kt
504    !
505      ENDIF
506      !
507   END SUBROUTINE dyn_nept_cor
508
509
510   SUBROUTINE dyn_nept_vel( ksign )
511      !!----------------------------------------------------------------------
512      !!                  ***  ROUTINE dyn_nept_vel  ***
513      !!
514      !! ** Purpose :  Add or subtract the Neptune velocity from the now
515      !!               velocities based on ksign
516      !!----------------------------------------------------------------------
517      INTEGER, INTENT( in ) ::   ksign    ! 1 or -1 to add or subtract neptune velocities
518      !!
519      INTEGER :: jk                       ! dummy loop index
520      !!----------------------------------------------------------------------
521      !
522      ! Adjust the current velocity un, vn by adding or subtracting the
523      ! Neptune velocities zunep, zvnep, as determined by argument ksign
524      DO jk=1, jpk
525         un(:,:,jk) = un(:,:,jk) + ksign * zunep(:,:) * umask(:,:,jk)
526         vn(:,:,jk) = vn(:,:,jk) + ksign * zvnep(:,:) * vmask(:,:,jk)
527      END DO
528      !
529   END SUBROUTINE dyn_nept_vel
530
531
532   SUBROUTINE dyn_nept_smooth_vel( htold, htnew, option )
533
534      !!----------------------------------------------------------------------
535      !!                  ***  ROUTINE dyn_nept_smooth_vel  ***
536      !!
537      !! ** Purpose :  Compute smoothed topography field.
538      !!
539      !! ** Action : - Updates the array htnew (output) with a smoothed
540      !!               version of the (input) array htold. Form of smoothing
541      !!               algorithm is controlled by the (logical) argument option.
542      !!----------------------------------------------------------------------
543
544      INTEGER                                   ::   ji, jj  ! dummy loop indices
545      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN)  ::   htold   ! temporary 2D workspace
546      LOGICAL, INTENT(IN)                       ::   option  ! temporary 2D workspace
547      REAL(wp), DIMENSION(jpi,jpj)              ::   htnew   ! temporary 2D workspace
548      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   work    ! temporary 2D workspace
549      INTEGER,  ALLOCATABLE, DIMENSION(:,:)     ::   nb      ! temporary 2D workspace
550      INTEGER,  ALLOCATABLE, DIMENSION(:,:)     ::   iwork   ! temporary 2D workspace
551
552!!    Dynamically allocate local work arrays
553      ALLOCATE( work(jpi,jpj), nb(jpi,jpj), iwork(jpi,jpj) )
554
555      iwork(:,:) = 0
556
557      !! iwork is a mask of gridpoints: iwork = 1 => ocean, iwork = 0 => land
558      WHERE( htold(:,:) > 0 )
559         iwork(:,:) = 1
560         htnew(:,:) = htold(:,:)
561      ELSEWHERE
562         iwork(:,:) = 0
563         htnew(:,:) = 0.0_wp
564      END WHERE
565      !! htnew contains valid ocean depths from htold, or zero
566
567      !! set work to a smoothed/averaged version of htnew; choice controlled by option
568      !! nb is set to the sum of the weights of the valid values used in work
569      IF( option ) THEN
570
571         !! Apply scale-selective smoothing in determining work from htnew
572         DO jj=2,jpj-1
573            DO ji=2,jpi-1
574               work(ji,jj) = 4.0*htnew( ji , jj ) +                        &
575                           & 2.0*htnew(ji+1, jj ) + 2.0*htnew(ji-1, jj ) + &
576                           & 2.0*htnew( ji ,jj+1) + 2.0*htnew( ji ,jj-1) + &
577                           &     htnew(ji+1,jj+1) +     htnew(ji+1,jj-1) + &
578                           &     htnew(ji-1,jj+1) +     htnew(ji-1,jj-1)
579
580               nb(ji,jj)   = 4 * iwork( ji , jj ) +                        &
581                           & 2 * iwork(ji+1, jj ) + 2 * iwork(ji-1, jj ) + &
582                           & 2 * iwork( ji ,jj+1) + 2 * iwork( ji ,jj-1) + &
583                           &     iwork(ji+1,jj+1) +     iwork(ji+1,jj-1) + &
584                           &     iwork(ji-1,jj+1) +     iwork(ji-1,jj-1)
585            END DO
586         END DO
587
588      ELSE
589
590         !! Apply simple 9-point averaging in determining work from htnew
591         DO jj=2,jpj-1
592            DO ji=2,jpi-1
593               work(ji,jj) =  htnew( ji , jj ) +                    &
594                           &  htnew(ji+1, jj ) + htnew(ji-1, jj ) + &
595                           &  htnew( ji ,jj+1) + htnew( ji ,jj-1) + &
596                           &  htnew(ji+1,jj+1) + htnew(ji+1,jj-1) + &
597                           &  htnew(ji-1,jj+1) + htnew(ji-1,jj-1)
598
599               nb(ji,jj) =    iwork( ji , jj ) +                    &
600                           &  iwork(ji+1, jj ) + iwork(ji-1, jj ) + &
601                           &  iwork( ji ,jj+1) + iwork( ji ,jj-1) + &
602                           &  iwork(ji+1,jj+1) + iwork(ji+1,jj-1) + &
603                           &  iwork(ji-1,jj+1) + iwork(ji-1,jj-1)
604            END DO
605         END DO
606
607      ENDIF
608
609      !! write averaged value of work into htnew,
610      !! if average is valid and point is unmasked
611      WHERE( (htold(:,:) /= 0.0_wp ) .AND. ( nb(:,:) /= 0 ) )
612         htnew(:,:) = work(:,:)/real(nb(:,:))
613      ELSEWHERE
614         htnew(:,:) = 0.0_wp
615      END WHERE
616
617!!    Deallocate temporary workspace arrays, all local to this routine
618      DEALLOCATE( work, nb, iwork )
619
620   END SUBROUTINE dyn_nept_smooth_vel
621
622END MODULE dynnept
Note: See TracBrowser for help on using the repository browser.