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/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90 @ 7649

Last change on this file since 7649 was 7649, checked in by glong, 8 years ago

dyn_vrt_dia subroutine added and calls added for most processes - spg, zdf, and bfr not yet done
File size: 21.6 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   !!----------------------------------------------------------------------
21
22   !!----------------------------------------------------------------------
23   !!   div_cur    : Compute the horizontal divergence and relative
24   !!                vorticity fields
25   !!----------------------------------------------------------------------
26   USE oce             ! ocean dynamics and tracers
27   USE dom_oce         ! ocean space and time domain
28   USE sbc_oce, ONLY : ln_rnf, nn_isf ! surface boundary condition: ocean
29   USE sbcrnf          ! river runoff
30   USE sbcisf          ! ice shelf
31   USE cla             ! cross land advection             (cla_div routine)
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   USE iom             ! I/O Manager for dyn_vrt_dia
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   div_cur     ! routine called by step.F90 and istate.F90
43   PUBLIC   dyn_vrt_dia ! routine called by various modules
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), ice shelf melting (div_isf)
73      !!      and cross land flow (div_cla)
74      !!              II. vorticity :
75      !!         - save the curl computed at the previous time-step
76      !!            rotb = rotn
77      !!      (note that the Asselin time filter has not been applied to rotb)
78      !!         - compute the now curl in tensorial formalism:
79      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
80      !!         - Coastal boundary condition: 'key_noslip_accurate' defined,
81      !!      the no-slip boundary condition is computed using Schchepetkin
82      !!      and O'Brien (1996) scheme (i.e. 4th order at the coast).
83      !!      For example, along east coast, the one-sided finite difference
84      !!      approximation used for di[v] is:
85      !!         di[e2v vn] =  1/(e1f*e2f) * ( (e2v vn)(i) + (e2v vn)(i-1) + (e2v vn)(i-2) )
86      !!
87      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
88      !!              - update rotb , rotn , the before & now rel. vorticity
89      !!----------------------------------------------------------------------
90      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
91      !
92      INTEGER ::   ji, jj, jk, jl           ! dummy loop indices
93      INTEGER ::   ii, ij, ijt, iju, ierr   ! local integer
94      REAL(wp) ::  zraur, zdep              ! local scalar
95      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwu   ! specific 2D workspace
96      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwv   ! specific 2D workspace
97      !!----------------------------------------------------------------------
98      !
99      IF( nn_timing == 1 )  CALL timing_start('div_cur')
100      !
101      CALL wrk_alloc( jpi  , jpj+2, zwu  )
102      CALL wrk_alloc( jpi+2, jpj  , zwv  )
103      !
104      IF( kt == nit000 ) THEN
105         IF(lwp) WRITE(numout,*)
106         IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and relative vorticity'
107         IF(lwp) WRITE(numout,*) '~~~~~~~   NOT optimal for auto-tasking case'
108      ENDIF
109
110      !                                                ! ===============
111      DO jk = 1, jpkm1                                 ! Horizontal slab
112         !                                             ! ===============
113         !
114         hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
115         rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
116         !
117         !                                             ! --------
118         ! Horizontal divergence                       !   div
119         !                                             ! --------
120         DO jj = 2, jpjm1
121            DO ji = fs_2, fs_jpim1   ! vector opt.
122               hdivn(ji,jj,jk) =   &
123                  (  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)       &
124                   + 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)  )    &
125                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
126            END DO
127         END DO
128
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( ln_divisf .AND. (nn_isf /= 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field)
234      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (Update Hor. divergence)
235     
236      ! 4. Lateral boundary conditions on hdivn and rotn
237      ! ---------------------------------=======---======
238      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )    ! lateral boundary cond. (no sign change)
239      !
240      CALL wrk_dealloc( jpi  , jpj+2, zwu )
241      CALL wrk_dealloc( jpi+2, jpj  , zwv )
242      !
243      IF( nn_timing == 1 )  CALL timing_stop('div_cur')
244      !
245   END SUBROUTINE div_cur
246   
247#else
248   !!----------------------------------------------------------------------
249   !!   Default option                           2nd order centered schemes
250   !!----------------------------------------------------------------------
251
252   SUBROUTINE div_cur( kt )
253      !!----------------------------------------------------------------------
254      !!                  ***  ROUTINE div_cur  ***
255      !!                   
256      !! ** Purpose :   compute the horizontal divergence and the relative
257      !!      vorticity at before and now time-step
258      !!
259      !! ** Method  : - Divergence:
260      !!      - save the divergence computed at the previous time-step
261      !!      (note that the Asselin filter has not been applied on hdivb)
262      !!      - compute the now divergence given by :
263      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
264      !!      correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla)
265      !!              - Relavtive Vorticity :
266      !!      - save the curl computed at the previous time-step (rotb = rotn)
267      !!      (note that the Asselin time filter has not been applied to rotb)
268      !!      - compute the now curl in tensorial formalism:
269      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
270      !!      Note: Coastal boundary condition: lateral friction set through
271      !!      the value of fmask along the coast (see dommsk.F90) and shlat
272      !!      (namelist parameter)
273      !!
274      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
275      !!              - update rotb , rotn , the before & now rel. vorticity
276      !!----------------------------------------------------------------------
277      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
278      !
279      INTEGER  ::   ji, jj, jk    ! dummy loop indices
280      REAL(wp) ::   zraur, zdep   ! local scalars
281      !!----------------------------------------------------------------------
282      !
283      IF( nn_timing == 1 )  CALL timing_start('div_cur')
284      !
285      IF( kt == nit000 ) THEN
286         IF(lwp) WRITE(numout,*)
287         IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and'
288         IF(lwp) WRITE(numout,*) '~~~~~~~   relative vorticity'
289      ENDIF
290
291      !                                                ! ===============
292      DO jk = 1, jpkm1                                 ! Horizontal slab
293         !                                             ! ===============
294         !
295         hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
296         rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
297         !
298         !                                             ! --------
299         ! Horizontal divergence                       !   div
300         !                                             ! --------
301         DO jj = 2, jpjm1
302            DO ji = fs_2, fs_jpim1   ! vector opt.
303               hdivn(ji,jj,jk) =   &
304                  (  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)       &
305                   + 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)  )    &
306                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
307            END DO 
308         END DO 
309
310         IF( .NOT. AGRIF_Root() ) THEN
311            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
312            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
313            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
314            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
315         ENDIF
316
317         !                                             ! --------
318         ! relative vorticity                          !   rot
319         !                                             ! --------
320         DO jj = 1, jpjm1
321            DO ji = 1, fs_jpim1   ! vector opt.
322               rotn(ji,jj,jk) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
323                  &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
324                  &           * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
325            END DO
326         END DO
327         !                                             ! ===============
328      END DO                                           !   End of slab
329      !                                                ! ===============
330
331      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )                            ! runoffs (update hdivn field)
332      IF( ln_divisf .AND. (nn_isf .GT. 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field)
333      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
334      !
335      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )     ! lateral boundary cond. (no sign change)
336      !
337      IF( nn_timing == 1 )  CALL timing_stop('div_cur')
338      !
339   END SUBROUTINE div_cur
340   
341#endif
342
343
344   SUBROUTINE dyn_vrt_dia( utend, vtend, id_dia_vor_int, id_dia_vor_mn)
345
346      !!----------------------------------------------------------------------
347      !!
348      !! ** Purpose :  compute the integral and mean vorticity tendencies.
349      !!
350      !! ** Action : a) Calculate the vertical integrals of utend & of vtend
351      !!                (u_int & v_int)
352      !!             b) Calculate the vorticity tendencies for the vertical
353      !!                integrals.
354      !!             c) Calculate the vertical means, u_mn, v_mn from u_int
355      !!                and v_int by dividing by the depth
356      !!             d) Calculate the vorticity tendencies for the vertical
357      !!                means
358      !!             e) Call iom_put for the vertical integral vorticity
359      !!                tendencies (using id_dia_vor_int)
360      !!             f) Call iom_put for the vertical mean vorticity
361      !!                tendencies (using id_dia_vor_mn)
362      !!
363      !!----------------------------------------------------------------------
364      REAL :: utend(jpi,jpj,jpk) ! contribution to du/dt
365      REAL :: vtend(jpi,jpj,jpk) ! contribution to dv/dt
366      INTEGER :: id_dia_vor_int  ! identifier for the vertical integral vorticity diagnostic
367      INTEGER :: id_dia_vor_mn   ! identifier for the vertical mean vorticity diagnostic
368      !
369      !!----------------------------------------------------------------------
370      !
371      INTEGER  ::   ji, jj, jk   ! dummy loop indices
372      INTEGER  ::   ikbu, ikbv   ! dummy loop indices
373      !
374      REAL(wp), POINTER, DIMENSION(:,:) :: u_int   ! u vertical integral
375      REAL(wp), POINTER, DIMENSION(:,:) :: v_int   ! v vertical integral
376      REAL(wp), POINTER, DIMENSION(:,:) :: u_mn    ! u vertical means
377      REAL(wp), POINTER, DIMENSION(:,:) :: v_mn    ! u vertical means
378      REAL(wp), POINTER, DIMENSION(:,:) :: vor_int ! vort trend of vert integrals
379      REAL(wp), POINTER, DIMENSION(:,:) :: vor_mn  ! vort trend of vert means
380
381      CALL wrk_alloc(jpi, jpj, u_int)
382      CALL wrk_alloc(jpi, jpj, v_int)
383      CALL wrk_alloc(jpi, jpj, u_mn)
384      CALL wrk_alloc(jpi, jpj, v_mn)
385      CALL wrk_alloc(jpi, jpj, vor_int)
386      CALL wrk_alloc(jpi, jpj, vor_mn)
387
388      u_int(:,:) = 0.0_wp
389      v_int(:,:) = 0.0_wp
390
391      !
392      ! Calculate the vertical integrals of utend & of vtend
393      !
394     
395      DO jk = 1,jpk
396          DO jj = 1,jpj
397              DO ji = 1,jpi
398                  u_int(ji,jj) = u_int(ji,jj) + utend(ji,jj,jk)*fse3u(ji,jj,jk)
399                  v_int(ji,jj) = v_int(ji,jj) + vtend(ji,jj,jk)*fse3v(ji,jj,jk)
400              END DO
401          END DO
402      END DO
403
404      !
405      ! Calculate the vorticity tendencies for the vertical integrals.
406      ! 1/e1e2 * ((e2*d(vtend)/dx) - (e1*d(utend)/dy))
407      !
408
409      DO jj = 1,jpjm1
410          DO ji = 1,jpim1
411              vor_int(ji,jj) = ( v_int(ji+1,jj) * e2v(ji+1,jj)     &
412                  &            - v_int(ji,jj)   * e2v(ji,jj)       &
413                  &            + u_int(ji,jj)   * e1u(ji,jj)       &
414                  &            - u_int(ji,jj+1) * e1u(ji,jj+1) )   &
415                  &           / ( e1f(ji,jj)    * e2f(ji,jj) )
416          END DO
417      END DO
418
419
420      !
421      ! Calculate the vertical means, u_mn, v_mn from u_int & v_int by dividing
422      ! by the depth
423      ! mbku & mbkv - vertical index of the bottom last U- & W- ocean level
424      !
425
426      DO jj = 1, jpj
427          DO ji = 1, jpi
428              ikbu = mbku(ji,jj)
429              ikbv = mbkv(ji,jj)
430
431              IF (ikbu .ne. 0.0_wp) THEN      ! Don't divide by 0!
432                  u_mn(ji,jj) = u_int(ji,jj) / ikbu
433              ELSE
434                  u_mn(ji,jj) = 0.0_wp
435              END IF
436
437              IF (ikbv .ne. 0.0_wp) THEN      ! Don't divide by 0!
438                  v_mn(ji,jj) = v_int(ji,jj) / ikbv
439              ELSE
440                  v_mn(ji,jj) = 0.0_wp
441              END IF
442          END DO
443      END DO
444
445      !
446      ! Calculate the vorticity tendencies for the vertical means
447      ! 1/e1e2 * ((e2*d(v_mn)/dx) - (e1*d(u_mn)/dy))
448      !
449
450      DO jj = 1,jpjm1
451          DO ji = 1,jpim1
452              vor_mn(ji,jj) = ( v_mn(ji+1,jj) * e2v(ji+1,jj)     &
453                  &           - v_mn(ji,jj)   * e2v(ji,jj)       &
454                  &           + u_mn(ji,jj)   * e1u(ji,jj)       &
455                  &           - u_mn(ji,jj+1) * e1u(ji,jj+1) )   &
456                  &          / ( e1f(ji,jj)   * e2f(ji,jj) )
457          END DO
458      END DO
459
460
461      ! Call iom_put for the vertical integral vorticity tendencies
462      IF (id_dia_vor_int == 1) THEN
463          CALL iom_put( "dia_vor_int", vor_int(:,:))
464      ENDIF
465
466      ! Call iom_put for the vertical mean vorticity tendencies
467      IF (id_dia_vor_int == 1) THEN
468          CALL iom_put( "dia_vor_mn", vor_mn(:,:))
469      ENDIF
470
471      CALL wrk_dealloc(jpi, jpj, u_int)
472      CALL wrk_dealloc(jpi, jpj, v_int)
473      CALL wrk_dealloc(jpi, jpj, u_mn)
474      CALL wrk_dealloc(jpi, jpj, v_mn)
475      CALL wrk_dealloc(jpi, jpj, vor_int)
476      CALL wrk_dealloc(jpi, jpj, vor_mn)
477
478   END SUBROUTINE dyn_vrt_dia
479
480
481   !!======================================================================
482END MODULE divcur
Note: See TracBrowser for help on using the repository browser.