source: branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90 @ 3876

Last change on this file since 3876 was 3876, checked in by gm, 8 years ago

dev_r3858_CNRS3_Ediag: #927 phasing with 2011/dev_r3309_LOCEAN12_Ediag branche + mxl diag update

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