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

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90 @ 3318

Last change on this file since 3318 was 3318, checked in by gm, 12 years ago

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