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

Last change on this file since 4417 was 4417, checked in by trackstand2, 10 years ago

Added debug SUM outputs to div_cur
  • Property svn:keywords set to Id
File size: 19.9 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
309#if defined key_obc
310      IF( Agrif_Root() ) THEN
311         ! open boundaries (div must be zero behind the open boundary)
312         !  mpp remark: The zeroing of hdivn can probably be extended to 1->jpi/jpj for the correct row/column
313         IF( lp_obc_east  )   hdivn(nie0p1:nie1p1,nje0  :nje1  ,1:jpkm1) = 0.e0      ! east
314         IF( lp_obc_west  )   hdivn(niw0  :niw1  ,njw0  :njw1  ,1:jpkm1) = 0.e0      ! west
315         IF( lp_obc_north )   hdivn(nin0  :nin1  ,njn0p1:njn1p1,1:jpkm1) = 0.e0      ! north
316         IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,1:jpkm1) = 0.e0      ! south
317      ENDIF
318#endif         
319      IF( .NOT. AGRIF_Root() ) THEN
320         IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,1:jpkm1) = 0.e0      ! east
321         IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,1:jpkm1) = 0.e0      ! west
322         IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,1:jpkm1) = 0.e0      ! north
323         IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,1:jpkm1) = 0.e0      ! south
324      ENDIF
325
326         !                                             ! --------
327         ! relative vorticity                          !   rot
328         !                                             ! --------
329      rotb (:,:,1:jpkm1) = rotn (:,:,1:jpkm1)    ! time swap of rot arrays
330      DO jj = 1, jpjm1
331         DO ji = 1, jpim1
332            DO jk = 1, jpkm1
333               rotn(ji,jj,jk) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
334                  &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
335                  &           * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
336            END DO
337         END DO
338      END DO
339#else
340
341      !                                                ! ===============
342      DO jk = 1, jpkm1                                 ! Horizontal slab
343         !                                             ! ===============
344         !
345         hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
346         rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
347         !
348         !                                             ! --------
349         ! Horizontal divergence                       !   div
350         !                                             ! --------
351         DO jj = 2, jpjm1
352            DO ji = fs_2, fs_jpim1   ! vector opt.
353               hdivn(ji,jj,jk) =   &
354                  (  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)       &
355                   + 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)  )    &
356                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
357            END DO 
358         END DO 
359
360#if defined key_obc
361         IF( Agrif_Root() ) THEN
362            ! open boundaries (div must be zero behind the open boundary)
363            !  mpp remark: The zeroing of hdivn can probably be extended to 1->jpi/jpj for the correct row/column
364            IF( lp_obc_east  )   hdivn(nie0p1:nie1p1,nje0  :nje1  ,jk) = 0.e0      ! east
365            IF( lp_obc_west  )   hdivn(niw0  :niw1  ,njw0  :njw1  ,jk) = 0.e0      ! west
366            IF( lp_obc_north )   hdivn(nin0  :nin1  ,njn0p1:njn1p1,jk) = 0.e0      ! north
367            IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.e0      ! south
368         ENDIF
369#endif         
370         IF( .NOT. AGRIF_Root() ) THEN
371            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
372            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
373            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
374            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
375         ENDIF
376
377         !                                             ! --------
378         ! relative vorticity                          !   rot
379         !                                             ! --------
380         DO jj = 1, jpjm1
381            DO ji = 1, fs_jpim1   ! vector opt.
382               rotn(ji,jj,jk) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
383                  &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
384                  &           * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
385            END DO
386         END DO
387         !                                             ! ===============
388      END DO                                           !   End of slab
389      !                                                ! ===============
390      WRITE(*,*) 'ARPDBG, div_cur: sum hdivn=',SUM(hdivn),'at step=',kt
391      WRITE(*,*) 'ARPDBG, div_cur: sum un=',SUM(un),'at step=',kt
392      WRITE(*,*) 'ARPDBG, div_cur: sum vn=',SUM(vn),'at step=',kt
393
394#endif
395
396      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
397      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
398      !
399      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )     ! lateral boundary cond. (no sign change)
400      !
401   END SUBROUTINE div_cur
402   
403#endif
404   !!======================================================================
405END MODULE divcur
Note: See TracBrowser for help on using the repository browser.