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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago

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