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/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90 @ 5770

Last change on this file since 5770 was 5770, checked in by gm, 9 years ago

#1593: LDF-ADV, step II.2: phasing the improvements/simplifications of advective tracer trend (see wiki page)
  • Property svn:keywords set to Id
File size: 16.1 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   !!            3.6  ! 2014-11  (P. Mathiot)          isf            added directly here
20   !!            3.7  ! 2015-10  (G. Madec) remove cross-land advection
21   !!----------------------------------------------------------------------
22
23   !!----------------------------------------------------------------------
24   !!   div_cur    : Compute the horizontal divergence and relative
25   !!                vorticity fields
26   !!----------------------------------------------------------------------
27   USE oce             ! ocean dynamics and tracers
28   USE dom_oce         ! ocean space and time domain
29   USE sbc_oce, ONLY : ln_rnf, nn_isf ! surface boundary condition: ocean
30   USE sbcrnf          ! river runoff
31   USE sbcisf          ! ice shelf
32   USE in_out_manager  ! I/O manager
33   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
34   USE lib_mpp         ! MPP library
35   USE wrk_nemo        ! Memory Allocation
36   USE timing          ! Timing
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   div_cur    ! routine called by step.F90 and istate.F90
42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
48   !! $Id$
49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53#if defined key_noslip_accurate
54   !!----------------------------------------------------------------------
55   !!   'key_noslip_accurate'   2nd order interior + 4th order at the coast
56   !!----------------------------------------------------------------------
57
58   SUBROUTINE div_cur( kt )
59      !!----------------------------------------------------------------------
60      !!                  ***  ROUTINE div_cur  ***
61      !!
62      !! ** Purpose :   compute the horizontal divergence and the relative
63      !!              vorticity at before and now time-step
64      !!
65      !! ** Method  : I.  divergence :
66      !!         - save the divergence computed at the previous time-step
67      !!      (note that the Asselin filter has not been applied on hdivb)
68      !!         - compute the now divergence given by :
69      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
70      !!      correct hdiv with runoff inflow (div_rnf) and ice shelf melting (div_isf)
71      !!              II. vorticity :
72      !!         - save the curl computed at the previous time-step
73      !!            rotb = rotn
74      !!      (note that the Asselin time filter has not been applied to rotb)
75      !!         - compute the now curl in tensorial formalism:
76      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
77      !!         - Coastal boundary condition: 'key_noslip_accurate' defined,
78      !!      the no-slip boundary condition is computed using Schchepetkin
79      !!      and O'Brien (1996) scheme (i.e. 4th order at the coast).
80      !!      For example, along east coast, the one-sided finite difference
81      !!      approximation used for di[v] is:
82      !!         di[e2v vn] =  1/(e1f*e2f) * ( (e2v vn)(i) + (e2v vn)(i-1) + (e2v vn)(i-2) )
83      !!
84      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
85      !!              - update rotb , rotn , the before & now rel. vorticity
86      !!----------------------------------------------------------------------
87      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
88      !
89      INTEGER ::   ji, jj, jk, jl           ! dummy loop indices
90      INTEGER ::   ii, ij, ijt, iju, ierr   ! local integer
91      REAL(wp) ::  zraur, zdep              ! local scalar
92      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwu   ! specific 2D workspace
93      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwv   ! specific 2D workspace
94      !!----------------------------------------------------------------------
95      !
96      IF( nn_timing == 1 )  CALL timing_start('div_cur')
97      !
98      CALL wrk_alloc( jpi  , jpj+2, zwu               )
99      CALL wrk_alloc( jpi+4, jpj  , zwv, kistart = -1 )
100      !
101      IF( kt == nit000 ) THEN
102         IF(lwp) WRITE(numout,*)
103         IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and relative vorticity'
104         IF(lwp) WRITE(numout,*) '~~~~~~~   NOT optimal for auto-tasking case'
105      ENDIF
106
107      !                                                ! ===============
108      DO jk = 1, jpkm1                                 ! Horizontal slab
109         !                                             ! ===============
110         !
111         hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
112         rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
113         !
114         !                                             ! --------
115         ! Horizontal divergence                       !   div
116         !                                             ! --------
117         DO jj = 2, jpjm1
118            DO ji = fs_2, fs_jpim1   ! vector opt.
119               hdivn(ji,jj,jk) =   &
120                  (  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)       &
121                   + 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)  )    &
122                  / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
123            END DO
124         END DO
125
126         IF( .NOT. AGRIF_Root() ) THEN
127            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
128            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
129            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
130            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
131         ENDIF
132
133         !                                             ! --------
134         ! relative vorticity                          !   rot
135         !                                             ! --------
136         ! contravariant velocity (extended for lateral b.c.)
137         ! inside the model domain
138         DO jj = 1, jpj
139            DO ji = 1, jpi
140               zwu(ji,jj) = e1u(ji,jj) * un(ji,jj,jk)
141               zwv(ji,jj) = e2v(ji,jj) * vn(ji,jj,jk)
142            END DO 
143         END DO 
144 
145         ! East-West boundary conditions
146         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN
147            zwv(  0  ,:) = zwv(jpi-2,:)
148            zwv( -1  ,:) = zwv(jpi-3,:)
149            zwv(jpi+1,:) = zwv(  3  ,:)
150            zwv(jpi+2,:) = zwv(  4  ,:)
151         ELSE
152            zwv(  0  ,:) = 0.e0
153            zwv( -1  ,:) = 0.e0
154            zwv(jpi+1,:) = 0.e0
155            zwv(jpi+2,:) = 0.e0
156         ENDIF
157
158         ! North-South boundary conditions
159         IF( nperio == 3 .OR. nperio == 4 ) THEN
160            ! north fold ( Grid defined with a T-point pivot) ORCA 2 degre
161            zwu(jpi,jpj+1) = 0.e0
162            zwu(jpi,jpj+2) = 0.e0
163            DO ji = 1, jpi-1
164               iju = jpi - ji + 1
165               zwu(ji,jpj+1) = - zwu(iju,jpj-3)
166               zwu(ji,jpj+2) = - zwu(iju,jpj-4)
167            END DO
168         ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN
169            ! north fold ( Grid defined with a F-point pivot) ORCA 0.5 degre\
170            zwu(jpi,jpj+1) = 0.e0
171            zwu(jpi,jpj+2) = 0.e0
172            DO ji = 1, jpi-1
173               iju = jpi - ji
174               zwu(ji,jpj  ) = - zwu(iju,jpj-1)
175               zwu(ji,jpj+1) = - zwu(iju,jpj-2)
176               zwu(ji,jpj+2) = - zwu(iju,jpj-3)
177            END DO
178            DO ji = -1, jpi+2
179               ijt = jpi - ji + 1
180               zwv(ji,jpj) = - zwv(ijt,jpj-2)
181            END DO
182            DO ji = jpi/2+1, jpi+2
183               ijt = jpi - ji + 1
184               zwv(ji,jpjm1) = - zwv(ijt,jpjm1)
185            END DO
186         ELSE
187            ! closed
188            zwu(:,jpj+1) = 0.e0
189            zwu(:,jpj+2) = 0.e0
190         ENDIF
191
192         ! relative vorticity (vertical component of the velocity curl)
193         DO jj = 1, jpjm1
194            DO ji = 1, fs_jpim1   ! vector opt.
195               rotn(ji,jj,jk) = (  zwv(ji+1,jj  ) - zwv(ji,jj)      &
196                  &              - zwu(ji  ,jj+1) + zwu(ji,jj)  ) * fmask(ji,jj,jk) * r1_e1e2f(ji,jj)
197            END DO
198         END DO
199
200         ! second order accurate scheme along straight coast
201         DO jl = 1, npcoa(1,jk)
202            ii = nicoa(jl,1,jk)
203            ij = njcoa(jl,1,jk)
204            rotn(ii,ij,jk) = r1_e1e2f(ji,jj) * ( + 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) = r1_e1e2f(ji,jj) * (-4.*zwv(ii,ij)+zwv(ii-1,ij)-0.2*zwv(ii-2,ij))
210         END DO
211         DO jl = 1, npcoa(3,jk)
212            ii = nicoa(jl,3,jk)
213            ij = njcoa(jl,3,jk)
214            rotn(ii,ij,jk) = -r1_e1e2f(ji,jj) * ( +4. * zwu(ii,ij+1) - zwu(ii,ij+2) + 0.2 * zwu(ii,ij+3) )
215         END DO
216         DO jl = 1, npcoa(4,jk)
217            ii = nicoa(jl,4,jk)
218            ij = njcoa(jl,4,jk)
219            rotn(ii,ij,jk) = -r1_e1e2f(ji,jj) * ( -4. * zwu(ii,ij) + zwu(ii,ij-1) - 0.2 * zwu(ii,ij-2) )
220         END DO
221         !                                             ! ===============
222      END DO                                           !   End of slab
223      !                                                ! ===============
224
225      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs   (update hdivn field)
226      IF( ln_divisf .AND. (nn_isf /= 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field)
227     
228      ! 4. Lateral boundary conditions on hdivn and rotn
229      ! ---------------------------------=======---======
230      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )    ! lateral boundary cond. (no sign change)
231      !
232      CALL wrk_dealloc( jpi  , jpj+2, zwu               )
233      CALL wrk_dealloc( jpi+4, jpj  , zwv, kistart = -1 )
234      !
235      IF( nn_timing == 1 )  CALL timing_stop('div_cur')
236      !
237   END SUBROUTINE div_cur
238   
239#else
240   !!----------------------------------------------------------------------
241   !!   Default option                           2nd order centered schemes
242   !!----------------------------------------------------------------------
243
244   SUBROUTINE div_cur( kt )
245      !!----------------------------------------------------------------------
246      !!                  ***  ROUTINE div_cur  ***
247      !!                   
248      !! ** Purpose :   compute the horizontal divergence and the relative
249      !!      vorticity at before and now time-step
250      !!
251      !! ** Method  : - Divergence:
252      !!      - save the divergence computed at the previous time-step
253      !!      (note that the Asselin filter has not been applied on hdivb)
254      !!      - compute the now divergence given by :
255      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
256      !!      correct hdiv with runoff inflow (div_rnf)
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      !!      Note: Coastal boundary condition: lateral friction set through
263      !!      the value of fmask along the coast (see dommsk.F90) and shlat
264      !!      (namelist parameter)
265      !!
266      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
267      !!              - update rotb , rotn , the before & now rel. vorticity
268      !!----------------------------------------------------------------------
269      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
270      !
271      INTEGER  ::   ji, jj, jk    ! dummy loop indices
272      REAL(wp) ::   zraur, zdep   ! local scalars
273      !!----------------------------------------------------------------------
274      !
275      IF( nn_timing == 1 )  CALL timing_start('div_cur')
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                  / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
299            END DO 
300         END DO 
301
302         IF( .NOT. AGRIF_Root() ) THEN
303            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
304            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
305            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
306            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
307         ENDIF
308
309         !                                             ! --------
310         ! relative vorticity                          !   rot
311         !                                             ! --------
312         DO jj = 1, jpjm1
313            DO ji = 1, fs_jpim1   ! vector opt.
314               rotn(ji,jj,jk) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
315                  &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
316                  &           * fmask(ji,jj,jk) * r1_e1e2f(ji,jj)
317            END DO
318         END DO
319         !                                             ! ===============
320      END DO                                           !   End of slab
321      !                                                ! ===============
322
323      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )                            ! runoffs (update hdivn field)
324      IF( ln_divisf .AND. (nn_isf .GT. 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field)
325      !
326      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )     ! lateral boundary cond. (no sign change)
327      !
328      IF( nn_timing == 1 )  CALL timing_stop('div_cur')
329      !
330   END SUBROUTINE div_cur
331   
332#endif
333   !!======================================================================
334END MODULE divcur
Note: See TracBrowser for help on using the repository browser.