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

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90 @ 9366

Last change on this file since 9366 was 9366, checked in by andmirek, 6 years ago

#2050 first version. Compiled OK in moci test suite

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