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

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90 @ 4372

Last change on this file since 4372 was 4372, checked in by jchanut, 10 years ago

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