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 branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/DYN – NEMO

source: branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/DYN/divcur.F90 @ 2207

Last change on this file since 2207 was 2207, checked in by acc, 14 years ago

#733 DEV_r2191_3partymerge2010. Merged in changes from devukmo2010 branch

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