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

source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90 @ 3953

Last change on this file since 3953 was 3723, checked in by acc, 11 years ago

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