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

source: branches/2011/dev_r2787_NOCS_NEPTUNE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90 @ 2795

Last change on this file since 2795 was 2795, checked in by acc, 13 years ago

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