source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 20 months ago

GMED 450 add flush after prints

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