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.
divcur.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/divcur.F90 @ 234

Last change on this file since 234 was 234, checked in by opalod, 19 years ago

CT : BUGFIX165 : no sign change on the vorticity field (rotn(:,:) field) at the North fold boundary
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.3 KB
Line 
1MODULE divcur
2   !!==============================================================================
3   !!                       ***  MODULE  divcur  ***
4   !! Ocean diagnostic variable : horizontal divergence and relative vorticity
5   !!==============================================================================
6
7   !!----------------------------------------------------------------------
8   !!   div_cur    : Compute the horizontal divergence and relative
9   !!                vorticity fields
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE oce             ! ocean dynamics and tracers
13   USE dom_oce         ! ocean space and time domain
14   USE in_out_manager  ! I/O manager
15   USE obc_oce         ! ocean lateral open boundary condition
16   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
17
18   IMPLICIT NONE
19   PRIVATE
20
21   !! * Accessibility
22   PUBLIC div_cur    ! routine called by step.F90 and istate.F90
23
24   !! * Substitutions
25#  include "domzgr_substitute.h90"
26#  include "vectopt_loop_substitute.h90"
27   !!----------------------------------------------------------------------
28   !!   OPA 9.0 , LODYC-IPSL  (2003)
29   !!----------------------------------------------------------------------
30
31CONTAINS
32
33#if defined key_noslip_accurate
34   !!----------------------------------------------------------------------
35   !!   'key_noslip_accurate'                     2nd order centered scheme
36   !!                                                4th order at the coast
37   !!----------------------------------------------------------------------
38
39   SUBROUTINE div_cur( kt )
40      !!----------------------------------------------------------------------
41      !!                  ***  ROUTINE div_cur  ***
42      !!
43      !! ** Purpose :   compute the horizontal divergence and the relative
44      !!      vorticity at before and now time-step
45      !!
46      !! ** Method  :
47      !!      I.  divergence :
48      !!         - save the divergence computed at the previous time-step
49      !!      (note that the Asselin filter has not been applied on hdivb)
50      !!         - compute the now divergence given by :
51      !!            * s-coordinate ('key_s_coord' defined)
52      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
53      !!         * z-coordinate (default key)
54      !!         hdivn = 1/(e1t*e2t) [ di(e2u  un) + dj(e1v  vn) ]
55      !!         - apply lateral boundary conditions on hdivn
56      !!      II. vorticity :
57      !!         - save the curl computed at the previous time-step
58      !!            rotb = rotn
59      !!      (note that the Asselin time filter has not been applied to rotb)
60      !!         - compute the now curl in tensorial formalism:
61      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
62      !!         - apply lateral boundary conditions on rotn through a call
63      !!      of lbc_lnk routine.
64      !!         - Coastal boundary condition: 'key_noslip_accurate' defined,
65      !!      the no-slip boundary condition is computed using Schchepetkin
66      !!      and O'Brien (1996) scheme (i.e. 4th order at the coast).
67      !!      For example, along east coast, the one-sided finite difference
68      !!      approximation used for di[v] is:
69      !!         di[e2v vn] =  1/(e1f*e2f)
70      !!                    * ( (e2v vn)(i) + (e2v vn)(i-1) + (e2v vn)(i-2) )
71      !!
72      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
73      !!              - update rotb , rotn , the before & now rel. vorticity
74      !!
75      !! History :
76      !!   8.2  !  00-03  (G. Madec)  no slip accurate
77      !!   9.0  !  03-08  (G. Madec)  merged of cur and div, free form, F90
78      !!----------------------------------------------------------------------
79      !! * Arguments
80      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
81     
82      !! * Local declarations
83      INTEGER ::   ji, jj, jk     ! dummy loop indices
84      INTEGER ::   ii, ij, jl     ! temporary integer
85      INTEGER ::   ijt, iju       ! temporary integer
86      REAL(wp) ::   zdiv, zdju
87      REAL(wp), DIMENSION(   jpi  ,1:jpj+2) ::   zwu   ! workspace
88      REAL(wp), DIMENSION(-1:jpi+2,  jpj  ) ::   zwv   ! workspace
89      !!----------------------------------------------------------------------
90
91      IF( kt == nit000 ) THEN
92         IF(lwp) WRITE(numout,*)
93         IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and relative vorticity'
94         IF(lwp) WRITE(numout,*) '~~~~~~~   NOT optimal for auto-tasking case'
95      ENDIF
96
97      !                                                ! ===============
98      DO jk = 1, jpkm1                                 ! Horizontal slab
99         !                                             ! ===============
100
101         hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
102         rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
103
104         !                                             ! --------
105         ! Horizontal divergence                       !   div
106         !                                             ! --------
107         DO jj = 2, jpjm1
108            DO ji = fs_2, fs_jpim1   ! vector opt.
109#if defined key_s_coord || defined key_partial_steps
110               hdivn(ji,jj,jk) =   &
111                  (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  )*fse3u(ji-1,jj  ,jk)  * un(ji-1,jj  ,jk)       &
112                   + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1)*fse3v(ji  ,jj-1,jk)  * vn(ji  ,jj-1,jk)  )    &
113                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
114#else
115               hdivn(ji,jj,jk) = (  e2u(ji,jj) * un(ji,jj,jk) - e2u(ji-1,jj  ) * un(ji-1,jj  ,jk)      &
116                  &               + e1v(ji,jj) * vn(ji,jj,jk) - e1v(ji  ,jj-1) * vn(ji  ,jj-1,jk)  )   &
117                  &            / ( e1t(ji,jj) * e2t(ji,jj) )
118#endif
119            END DO
120         END DO
121
122#if defined key_obc
123         ! open boundaries (div must be zero behind the open boundary)
124         !  mpp remark: The zeroing of hdivn can probably be extended to 1->jpi/jpj for the correct row/column
125         IF( lp_obc_east  )   hdivn(nie0p1:nie1p1,nje0  :nje1  ,jk) = 0.e0      ! east
126         IF( lp_obc_west  )   hdivn(niw0  :niw1  ,njw0  :njw1  ,jk) = 0.e0      ! west
127         IF( lp_obc_north )   hdivn(nin0  :nin1  ,njn0p1:njn1p1,jk) = 0.e0      ! north
128         IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.e0      ! south
129#endif         
130
131         !                                             ! --------
132         ! relative vorticity                          !   rot
133         !                                             ! --------
134         ! contravariant velocity (extended for lateral b.c.)
135         ! inside the model domain
136         DO jj = 1, jpj
137            DO ji = 1, jpi
138               zwu(ji,jj) = e1u(ji,jj) * un(ji,jj,jk)
139               zwv(ji,jj) = e2v(ji,jj) * vn(ji,jj,jk)
140            END DO 
141         END DO 
142 
143         ! East-West boundary conditions
144         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN
145            zwv(  0  ,:) = zwv(jpi-2,:)
146            zwv( -1  ,:) = zwv(jpi-3,:)
147            zwv(jpi+1,:) = zwv(  3  ,:)
148            zwv(jpi+2,:) = zwv(  4  ,:)
149         ELSE
150            zwv(  0  ,:) = 0.e0
151            zwv( -1  ,:) = 0.e0
152            zwv(jpi+1,:) = 0.e0
153            zwv(jpi+2,:) = 0.e0
154         ENDIF
155
156         ! North-South boundary conditions
157         IF( nperio == 3 .OR. nperio == 4 ) THEN
158            ! north fold ( Grid defined with a T-point pivot) ORCA 2 degre
159            zwu(jpi,jpj+1) = 0.e0
160            zwu(jpi,jpj+2) = 0.e0
161            DO ji = 1, jpi-1
162               iju = jpi - ji + 1
163               zwu(ji,jpj+1) = - zwu(iju,jpj-3)
164               zwu(ji,jpj+2) = - zwu(iju,jpj-4)
165            END DO
166         ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN
167            ! north fold ( Grid defined with a F-point pivot) ORCA 0.5 degre\
168            zwu(jpi,jpj+1) = 0.e0
169            zwu(jpi,jpj+2) = 0.e0
170            DO ji = 1, jpi-1
171               iju = jpi - ji
172               zwu(ji,jpj  ) = - zwu(iju,jpj-1)
173               zwu(ji,jpj+1) = - zwu(iju,jpj-2)
174               zwu(ji,jpj+2) = - zwu(iju,jpj-3)
175            END DO
176            DO ji = -1, jpi+2
177               ijt = jpi - ji + 1
178               zwv(ji,jpj) = - zwv(ijt,jpj-2)
179            END DO
180            DO ji = jpi/2+1, jpi+2
181               ijt = jpi - ji + 1
182               zwv(ji,jpjm1) = - zwv(ijt,jpjm1)
183            END DO
184         ELSE
185            ! closed
186            zwu(:,jpj+1) = 0.e0
187            zwu(:,jpj+2) = 0.e0
188         ENDIF
189
190         ! relative vorticity (vertical component of the velocity curl)
191         DO jj = 1, jpjm1
192            DO ji = 1, fs_jpim1   ! vector opt.
193               rotn(ji,jj,jk) = (  zwv(ji+1,jj  ) - zwv(ji,jj)      &
194                                 - zwu(ji  ,jj+1) + zwu(ji,jj)  )   &
195                              * fmask(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
196            END DO
197         END DO
198
199         ! second order accurate scheme along straight coast
200         DO jl = 1, npcoa(1,jk)
201            ii = nicoa(jl,1,jk)
202            ij = njcoa(jl,1,jk)
203            rotn(ii,ij,jk) = 1. / ( e1f(ii,ij) * e2f(ii,ij) )   &
204                           * ( + 4. * zwv(ii+1,ij) - zwv(ii+2,ij) + 0.2 * zwv(ii+3,ij) )
205         END DO
206         DO jl = 1, npcoa(2,jk)
207            ii = nicoa(jl,2,jk)
208            ij = njcoa(jl,2,jk)
209            rotn(ii,ij,jk) = 1./(e1f(ii,ij)*e2f(ii,ij))   &
210               *(-4.*zwv(ii,ij)+zwv(ii-1,ij)-0.2*zwv(ii-2,ij))
211         END DO
212         DO jl = 1, npcoa(3,jk)
213            ii = nicoa(jl,3,jk)
214            ij = njcoa(jl,3,jk)
215            rotn(ii,ij,jk) = -1. / ( e1f(ii,ij)*e2f(ii,ij) )   &
216               * ( +4. * zwu(ii,ij+1) - zwu(ii,ij+2) + 0.2 * zwu(ii,ij+3) )
217         END DO
218         DO jl = 1, npcoa(4,jk)
219            ii = nicoa(jl,4,jk)
220            ij = njcoa(jl,4,jk)
221            rotn(ii,ij,jk) = -1. / ( e1f(ii,ij)*e2f(ii,ij) )   &
222               * ( -4. * zwu(ii,ij) + zwu(ii,ij-1) - 0.2 * zwu(ii,ij-2) )
223         END DO
224
225         !                                             ! ===============
226      END DO                                           !   End of slab
227      !                                                ! ===============
228     
229      ! 4. Lateral boundary conditions on hdivn and rotn
230      ! ---------------------------------=======---======
231      CALL lbc_lnk( hdivn, 'T', 1. )     ! T-point, no sign change
232      CALL lbc_lnk( rotn , 'F', 1. )     ! F-point, no sign change
233
234   END SUBROUTINE div_cur
235   
236#else
237   !!----------------------------------------------------------------------
238   !!   Default option                           2nd order centered schemes
239   !!----------------------------------------------------------------------
240
241   SUBROUTINE div_cur( kt )
242      !!----------------------------------------------------------------------
243      !!                  ***  ROUTINE div_cur  ***
244      !!                   
245      !! ** Purpose :   compute the horizontal divergence and the relative
246      !!      vorticity at before and now time-step
247      !!
248      !! ** Method  : - Divergence:
249      !!      - save the divergence computed at the previous time-step
250      !!      (note that the Asselin filter has not been applied on hdivb)
251      !!      - compute the now divergence given by :
252      !!         * s-coordinate ('key_s_coord' defined)
253      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
254      !!         * z-coordinate (default key)
255      !!         hdivn = 1/(e1t*e2t) [ di(e2u  un) + dj(e1v  vn) ]
256      !!      - apply lateral boundary conditions on hdivn
257      !!              - Relavtive Vorticity :
258      !!      - save the curl computed at the previous time-step (rotb = rotn)
259      !!      (note that the Asselin time filter has not been applied to rotb)
260      !!      - compute the now curl in tensorial formalism:
261      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
262      !!      - apply lateral boundary conditions on rotn through a call to
263      !!      routine lbc_lnk routine.
264      !!      Note: Coastal boundary condition: lateral friction set through
265      !!      the value of fmask along the coast (see dommsk.F90) and shlat
266      !!      (namelist parameter)
267      !!
268      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
269      !!              - update rotb , rotn , the before & now rel. vorticity
270      !!
271      !! History :
272      !!   1.0  !  87-06  (P. Andrich, D. L Hostis)  Original code
273      !!   4.0  !  91-11  (G. Madec)
274      !!   6.0  !  93-03  (M. Guyon)  symetrical conditions
275      !!   7.0  !  96-01  (G. Madec)  s-coordinates
276      !!   8.0  !  97-06  (G. Madec)  lateral boundary cond., lbc
277      !!   8.1  !  97-08  (J.M. Molines)  Open boundaries
278      !!   9.0  !  02-09  (G. Madec, E. Durand)  Free form, F90
279      !!----------------------------------------------------------------------
280      !! * Arguments
281      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
282     
283      !! * Local declarations
284      INTEGER  ::   ji, jj, jk          ! dummy loop indices
285      !!----------------------------------------------------------------------
286
287      IF( kt == nit000 ) THEN
288         IF(lwp) WRITE(numout,*)
289         IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and'
290         IF(lwp) WRITE(numout,*) '~~~~~~~   relative vorticity'
291      ENDIF
292
293      !                                                ! ===============
294      DO jk = 1, jpkm1                                 ! Horizontal slab
295         !                                             ! ===============
296
297         hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
298         rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
299
300         !                                             ! --------
301         ! Horizontal divergence                       !   div
302         !                                             ! --------
303         DO jj = 2, jpjm1
304            DO ji = fs_2, fs_jpim1   ! vector opt.
305#if defined key_s_coord || defined key_partial_steps
306               hdivn(ji,jj,jk) =   &
307                  (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  )*fse3u(ji-1,jj  ,jk)  * un(ji-1,jj  ,jk)       &
308                   + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1)*fse3v(ji  ,jj-1,jk)  * vn(ji  ,jj-1,jk)  )    &
309                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
310#else
311               hdivn(ji,jj,jk) = (  e2u(ji,jj) * un(ji,jj,jk) - e2u(ji-1,jj  ) * un(ji-1,jj  ,jk)      &
312                  &               + e1v(ji,jj) * vn(ji,jj,jk) - e1v(ji  ,jj-1) * vn(ji  ,jj-1,jk)  )   & 
313                  / ( e1t(ji,jj) * e2t(ji,jj) )
314#endif
315            END DO 
316         END DO 
317
318#if defined key_obc
319         ! open boundaries (div must be zero behind the open boundary)
320         !  mpp remark: The zeroing of hdivn can probably be extended to 1->jpi/jpj for the correct row/column
321         IF( lp_obc_east  )   hdivn(nie0p1:nie1p1,nje0  :nje1  ,jk) = 0.e0      ! east
322         IF( lp_obc_west  )   hdivn(niw0  :niw1  ,njw0  :njw1  ,jk) = 0.e0      ! west
323         IF( lp_obc_north )   hdivn(nin0  :nin1  ,njn0p1:njn1p1,jk) = 0.e0      ! north
324         IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.e0      ! south
325#endif         
326         !                                             ! --------
327         ! relative vorticity                          !   rot
328         !                                             ! --------
329         DO jj = 1, jpjm1
330            DO ji = 1, fs_jpim1   ! vector opt.
331               rotn(ji,jj,jk) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
332                  &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
333                  &           * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
334            END DO
335         END DO
336         !                                             ! ===============
337      END DO                                           !   End of slab
338      !                                                ! ===============
339     
340      ! 4. Lateral boundary conditions on hdivn and rotn
341      ! ---------------------------------=======---======
342      CALL lbc_lnk( hdivn, 'T', 1. )       ! T-point, no sign change
343      CALL lbc_lnk( rotn , 'F', 1. )       ! F-point, no sign change
344
345   END SUBROUTINE div_cur
346   
347#endif
348   !!======================================================================
349END MODULE divcur
Note: See TracBrowser for help on using the repository browser.