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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 19.7 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   USE lib_mpp         ! MPP library
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   div_cur    ! routine called by step.F90 and istate.F90
39
40   !! * Control permutation of array indices
41#  include "oce_ftrans.h90"
42#  include "dom_oce_ftrans.h90"
43#  include "obc_oce_ftrans.h90"
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
50   !! $Id$
51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55#if defined key_noslip_accurate
56   !!----------------------------------------------------------------------
57   !!   'key_noslip_accurate'   2nd order interior + 4th order at the coast
58   !!----------------------------------------------------------------------
59
60   SUBROUTINE div_cur( kt )
61      !!----------------------------------------------------------------------
62      !!                  ***  ROUTINE div_cur  ***
63      !!
64      !! ** Purpose :   compute the horizontal divergence and the relative
65      !!              vorticity at before and now time-step
66      !!
67      !! ** Method  : I.  divergence :
68      !!         - save the divergence computed at the previous time-step
69      !!      (note that the Asselin filter has not been applied on hdivb)
70      !!         - compute the now divergence given by :
71      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
72      !!      correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla)
73      !!              II. vorticity :
74      !!         - save the curl computed at the previous time-step
75      !!            rotb = rotn
76      !!      (note that the Asselin time filter has not been applied to rotb)
77      !!         - compute the now curl in tensorial formalism:
78      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
79      !!         - Coastal boundary condition: 'key_noslip_accurate' defined,
80      !!      the no-slip boundary condition is computed using Schchepetkin
81      !!      and O'Brien (1996) scheme (i.e. 4th order at the coast).
82      !!      For example, along east coast, the one-sided finite difference
83      !!      approximation used for di[v] is:
84      !!         di[e2v vn] =  1/(e1f*e2f) * ( (e2v vn)(i) + (e2v vn)(i-1) + (e2v vn)(i-2) )
85      !!
86      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
87      !!              - update rotb , rotn , the before & now rel. vorticity
88      !!----------------------------------------------------------------------
89      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
90      !
91      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zwu   ! specific 2D workspace
92      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zwv   ! specific 2D workspace
93      !
94      INTEGER ::   ji, jj, jk, jl           ! dummy loop indices
95      INTEGER ::   ii, ij, ijt, iju, ierr   ! local integer
96      REAL(wp) ::  zraur, zdep              ! local scalar
97      !!----------------------------------------------------------------------
98
99      IF( kt == nit000 ) THEN
100         IF(lwp) WRITE(numout,*)
101         IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and relative vorticity'
102         IF(lwp) WRITE(numout,*) '~~~~~~~   NOT optimal for auto-tasking case'
103         !
104         ALLOCATE( zwu( jpi, 1:jpj+2) , zwv(-1:jpi+2, jpj) , STAT=ierr )
105         IF( lk_mpp    )   CALL mpp_sum( ierr )
106         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'div_cur : unable to allocate arrays' )
107      ENDIF
108
109      !                                                ! ===============
110      DO jk = 1, jpkm1                                 ! Horizontal slab
111         !                                             ! ===============
112         !
113         hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
114         rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
115         !
116         !                                             ! --------
117         ! Horizontal divergence                       !   div
118         !                                             ! --------
119         DO jj = 2, jpjm1
120            DO ji = fs_2, fs_jpim1   ! vector opt.
121               hdivn(ji,jj,jk) =   &
122                  (  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)       &
123                   + 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)  )    &
124                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
125            END DO
126         END DO
127
128#if defined key_obc
129         IF( Agrif_Root() ) THEN
130            ! open boundaries (div must be zero behind the open boundary)
131            !  mpp remark: The zeroing of hdivn can probably be extended to 1->jpi/jpj for the correct row/column
132            IF( lp_obc_east  )   hdivn(nie0p1:nie1p1,nje0  :nje1  ,jk) = 0.e0      ! east
133            IF( lp_obc_west  )   hdivn(niw0  :niw1  ,njw0  :njw1  ,jk) = 0.e0      ! west
134            IF( lp_obc_north )   hdivn(nin0  :nin1  ,njn0p1:njn1p1,jk) = 0.e0      ! north
135            IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.e0      ! south
136         ENDIF
137#endif         
138         IF( .NOT. AGRIF_Root() ) THEN
139            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
140            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
141            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
142            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
143         ENDIF
144
145         !                                             ! --------
146         ! relative vorticity                          !   rot
147         !                                             ! --------
148         ! contravariant velocity (extended for lateral b.c.)
149         ! inside the model domain
150         DO jj = 1, jpj
151            DO ji = 1, jpi
152               zwu(ji,jj) = e1u(ji,jj) * un(ji,jj,jk)
153               zwv(ji,jj) = e2v(ji,jj) * vn(ji,jj,jk)
154            END DO 
155         END DO 
156 
157         ! East-West boundary conditions
158         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN
159            zwv(  0  ,:) = zwv(jpi-2,:)
160            zwv( -1  ,:) = zwv(jpi-3,:)
161            zwv(jpi+1,:) = zwv(  3  ,:)
162            zwv(jpi+2,:) = zwv(  4  ,:)
163         ELSE
164            zwv(  0  ,:) = 0.e0
165            zwv( -1  ,:) = 0.e0
166            zwv(jpi+1,:) = 0.e0
167            zwv(jpi+2,:) = 0.e0
168         ENDIF
169
170         ! North-South boundary conditions
171         IF( nperio == 3 .OR. nperio == 4 ) THEN
172            ! north fold ( Grid defined with a T-point pivot) ORCA 2 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 + 1
177               zwu(ji,jpj+1) = - zwu(iju,jpj-3)
178               zwu(ji,jpj+2) = - zwu(iju,jpj-4)
179            END DO
180         ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN
181            ! north fold ( Grid defined with a F-point pivot) ORCA 0.5 degre\
182            zwu(jpi,jpj+1) = 0.e0
183            zwu(jpi,jpj+2) = 0.e0
184            DO ji = 1, jpi-1
185               iju = jpi - ji
186               zwu(ji,jpj  ) = - zwu(iju,jpj-1)
187               zwu(ji,jpj+1) = - zwu(iju,jpj-2)
188               zwu(ji,jpj+2) = - zwu(iju,jpj-3)
189            END DO
190            DO ji = -1, jpi+2
191               ijt = jpi - ji + 1
192               zwv(ji,jpj) = - zwv(ijt,jpj-2)
193            END DO
194            DO ji = jpi/2+1, jpi+2
195               ijt = jpi - ji + 1
196               zwv(ji,jpjm1) = - zwv(ijt,jpjm1)
197            END DO
198         ELSE
199            ! closed
200            zwu(:,jpj+1) = 0.e0
201            zwu(:,jpj+2) = 0.e0
202         ENDIF
203
204         ! relative vorticity (vertical component of the velocity curl)
205         DO jj = 1, jpjm1
206            DO ji = 1, fs_jpim1   ! vector opt.
207               rotn(ji,jj,jk) = (  zwv(ji+1,jj  ) - zwv(ji,jj)      &
208                  &              - zwu(ji  ,jj+1) + zwu(ji,jj)  ) * fmask(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
209            END DO
210         END DO
211
212         ! second order accurate scheme along straight coast
213         DO jl = 1, npcoa(1,jk)
214            ii = nicoa(jl,1,jk)
215            ij = njcoa(jl,1,jk)
216            rotn(ii,ij,jk) = 1. / ( e1f(ii,ij) * e2f(ii,ij) )   &
217                           * ( + 4. * zwv(ii+1,ij) - zwv(ii+2,ij) + 0.2 * zwv(ii+3,ij) )
218         END DO
219         DO jl = 1, npcoa(2,jk)
220            ii = nicoa(jl,2,jk)
221            ij = njcoa(jl,2,jk)
222            rotn(ii,ij,jk) = 1./(e1f(ii,ij)*e2f(ii,ij))   &
223               *(-4.*zwv(ii,ij)+zwv(ii-1,ij)-0.2*zwv(ii-2,ij))
224         END DO
225         DO jl = 1, npcoa(3,jk)
226            ii = nicoa(jl,3,jk)
227            ij = njcoa(jl,3,jk)
228            rotn(ii,ij,jk) = -1. / ( e1f(ii,ij)*e2f(ii,ij) )   &
229               * ( +4. * zwu(ii,ij+1) - zwu(ii,ij+2) + 0.2 * zwu(ii,ij+3) )
230         END DO
231         DO jl = 1, npcoa(4,jk)
232            ii = nicoa(jl,4,jk)
233            ij = njcoa(jl,4,jk)
234            rotn(ii,ij,jk) = -1. / ( e1f(ii,ij)*e2f(ii,ij) )   &
235               * ( -4. * zwu(ii,ij) + zwu(ii,ij-1) - 0.2 * zwu(ii,ij-2) )
236         END DO
237         !                                             ! ===============
238      END DO                                           !   End of slab
239      !                                                ! ===============
240
241      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
242      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (Update Hor. divergence)
243     
244      ! 4. Lateral boundary conditions on hdivn and rotn
245      ! ---------------------------------=======---======
246      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )    ! lateral boundary cond. (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      !!      correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla)
268      !!              - Relavtive Vorticity :
269      !!      - save the curl computed at the previous time-step (rotb = rotn)
270      !!      (note that the Asselin time filter has not been applied to rotb)
271      !!      - compute the now curl in tensorial formalism:
272      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
273      !!      Note: Coastal boundary condition: lateral friction set through
274      !!      the value of fmask along the coast (see dommsk.F90) and shlat
275      !!      (namelist parameter)
276      !!
277      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
278      !!              - update rotb , rotn , the before & now rel. vorticity
279      !!----------------------------------------------------------------------
280      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
281      !
282      INTEGER  ::   ji, jj, jk    ! dummy loop indices
283      REAL(wp) ::   zraur, zdep   ! local scalars
284      !!----------------------------------------------------------------------
285
286      IF( kt == nit000 ) THEN
287         IF(lwp) WRITE(numout,*)
288         IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and'
289         IF(lwp) WRITE(numout,*) '~~~~~~~   relative vorticity'
290      ENDIF
291
292#if defined key_z_first
293      !                                             ! --------
294      ! Horizontal divergence                       !   div
295      !                                             ! --------
296      hdivb(:,:,1:jpkm1) = hdivn(:,:,1:jpkm1)    ! time swap of div arrays
297      DO jj = 2, jpjm1
298         DO ji = 2, jpim1
299            DO jk = 1, jpkm1
300               hdivn(ji,jj,jk) =   &
301                  (  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)       &
302                   + 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)  )    &
303                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
304            END DO 
305         END DO 
306      END DO
307
308#if defined key_obc
309      IF( Agrif_Root() ) THEN
310         ! open boundaries (div must be zero behind the open boundary)
311         !  mpp remark: The zeroing of hdivn can probably be extended to 1->jpi/jpj for the correct row/column
312         IF( lp_obc_east  )   hdivn(nie0p1:nie1p1,nje0  :nje1  ,1:jpkm1) = 0.e0      ! east
313         IF( lp_obc_west  )   hdivn(niw0  :niw1  ,njw0  :njw1  ,1:jpkm1) = 0.e0      ! west
314         IF( lp_obc_north )   hdivn(nin0  :nin1  ,njn0p1:njn1p1,1:jpkm1) = 0.e0      ! north
315         IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,1:jpkm1) = 0.e0      ! south
316      ENDIF
317#endif         
318      IF( .NOT. AGRIF_Root() ) THEN
319         IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,1:jpkm1) = 0.e0      ! east
320         IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,1:jpkm1) = 0.e0      ! west
321         IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,1:jpkm1) = 0.e0      ! north
322         IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,1:jpkm1) = 0.e0      ! south
323      ENDIF
324
325         !                                             ! --------
326         ! relative vorticity                          !   rot
327         !                                             ! --------
328      rotb (:,:,1:jpkm1) = rotn (:,:,1:jpkm1)    ! time swap of rot arrays
329      DO jj = 1, jpjm1
330         DO ji = 1, jpim1
331            DO jk = 1, jpkm1
332               rotn(ji,jj,jk) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
333                  &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
334                  &           * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
335            END DO
336         END DO
337      END DO
338#else
339
340      !                                                ! ===============
341      DO jk = 1, jpkm1                                 ! Horizontal slab
342         !                                             ! ===============
343         !
344         hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
345         rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
346         !
347         !                                             ! --------
348         ! Horizontal divergence                       !   div
349         !                                             ! --------
350         DO jj = 2, jpjm1
351            DO ji = fs_2, fs_jpim1   ! vector opt.
352               hdivn(ji,jj,jk) =   &
353                  (  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)       &
354                   + 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)  )    &
355                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
356            END DO 
357         END DO 
358
359#if defined key_obc
360         IF( Agrif_Root() ) THEN
361            ! open boundaries (div must be zero behind the open boundary)
362            !  mpp remark: The zeroing of hdivn can probably be extended to 1->jpi/jpj for the correct row/column
363            IF( lp_obc_east  )   hdivn(nie0p1:nie1p1,nje0  :nje1  ,jk) = 0.e0      ! east
364            IF( lp_obc_west  )   hdivn(niw0  :niw1  ,njw0  :njw1  ,jk) = 0.e0      ! west
365            IF( lp_obc_north )   hdivn(nin0  :nin1  ,njn0p1:njn1p1,jk) = 0.e0      ! north
366            IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.e0      ! south
367         ENDIF
368#endif         
369         IF( .NOT. AGRIF_Root() ) THEN
370            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
371            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
372            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
373            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
374         ENDIF
375
376         !                                             ! --------
377         ! relative vorticity                          !   rot
378         !                                             ! --------
379         DO jj = 1, jpjm1
380            DO ji = 1, fs_jpim1   ! vector opt.
381               rotn(ji,jj,jk) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
382                  &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
383                  &           * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
384            END DO
385         END DO
386         !                                             ! ===============
387      END DO                                           !   End of slab
388      !                                                ! ===============
389#endif
390
391      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
392      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
393      !
394      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )     ! lateral boundary cond. (no sign change)
395      !
396   END SUBROUTINE div_cur
397   
398#endif
399   !!======================================================================
400END MODULE divcur
Note: See TracBrowser for help on using the repository browser.