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 @ 1058

Last change on this file since 1058 was 1058, checked in by rblod, 16 years ago

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