source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90 @ 2392

Last change on this file since 2392 was 2392, checked in by gm, 10 years ago

v3.3beta: Cross Land Advection (ticket #127) full rewriting + MPP bug corrections

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