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

Last change on this file since 8197 was 8197, checked in by glong, 7 years ago

Changed id's to be chars e.g. hpg to more easily identify output (and updated field_def.xml accordingly). Also rearranged scaling factors in dyn_vrt_dia subroutines in divcur.F90

File size: 27.4 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_2d
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   div_cur        ! routine called by step.F90 and istate.F90
43   PUBLIC   dyn_vrt_dia_3d ! routine called by various modules
44   PUBLIC   dyn_vrt_dia_2d ! routine called by various modules
45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
51   !! $Id$
52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56#if defined key_noslip_accurate
57   !!----------------------------------------------------------------------
58   !!   'key_noslip_accurate'   2nd order interior + 4th order at the coast
59   !!----------------------------------------------------------------------
60
61   SUBROUTINE div_cur( kt )
62      !!----------------------------------------------------------------------
63      !!                  ***  ROUTINE div_cur  ***
64      !!
65      !! ** Purpose :   compute the horizontal divergence and the relative
66      !!              vorticity at before and now time-step
67      !!
68      !! ** Method  : I.  divergence :
69      !!         - save the divergence computed at the previous time-step
70      !!      (note that the Asselin filter has not been applied on hdivb)
71      !!         - compute the now divergence given by :
72      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
73      !!      correct hdiv with runoff inflow (div_rnf), ice shelf melting (div_isf)
74      !!      and cross land flow (div_cla)
75      !!              II. vorticity :
76      !!         - save the curl computed at the previous time-step
77      !!            rotb = rotn
78      !!      (note that the Asselin time filter has not been applied to rotb)
79      !!         - compute the now curl in tensorial formalism:
80      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
81      !!         - Coastal boundary condition: 'key_noslip_accurate' defined,
82      !!      the no-slip boundary condition is computed using Schchepetkin
83      !!      and O'Brien (1996) scheme (i.e. 4th order at the coast).
84      !!      For example, along east coast, the one-sided finite difference
85      !!      approximation used for di[v] is:
86      !!         di[e2v vn] =  1/(e1f*e2f) * ( (e2v vn)(i) + (e2v vn)(i-1) + (e2v vn)(i-2) )
87      !!
88      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
89      !!              - update rotb , rotn , the before & now rel. vorticity
90      !!----------------------------------------------------------------------
91      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
92      !
93      INTEGER ::   ji, jj, jk, jl           ! dummy loop indices
94      INTEGER ::   ii, ij, ijt, iju, ierr   ! local integer
95      REAL(wp) ::  zraur, zdep              ! local scalar
96      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwu   ! specific 2D workspace
97      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwv   ! specific 2D workspace
98      !!----------------------------------------------------------------------
99      !
100      IF( nn_timing == 1 )  CALL timing_start('div_cur')
101      !
102      CALL wrk_alloc( jpi  , jpj+2, zwu  )
103      CALL wrk_alloc( jpi+2, jpj  , zwv  )
104      !
105      IF( kt == nit000 ) THEN
106         IF(lwp) WRITE(numout,*)
107         IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and relative vorticity'
108         IF(lwp) WRITE(numout,*) '~~~~~~~   NOT optimal for auto-tasking case'
109      ENDIF
110
111      !                                                ! ===============
112      DO jk = 1, jpkm1                                 ! Horizontal slab
113         !                                             ! ===============
114         !
115         hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
116         rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
117         !
118         !                                             ! --------
119         ! Horizontal divergence                       !   div
120         !                                             ! --------
121         DO jj = 2, jpjm1
122            DO ji = fs_2, fs_jpim1   ! vector opt.
123               hdivn(ji,jj,jk) =   &
124                  (  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)       &
125                   + 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)  )    &
126                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
127            END DO
128         END DO
129
130         IF( .NOT. AGRIF_Root() ) THEN
131            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
132            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
133            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
134            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
135         ENDIF
136
137         !                                             ! --------
138         ! relative vorticity                          !   rot
139         !                                             ! --------
140         ! contravariant velocity (extended for lateral b.c.)
141         ! inside the model domain
142         DO jj = 1, jpj
143            DO ji = 1, jpi
144               zwu(ji,jj) = e1u(ji,jj) * un(ji,jj,jk)
145               zwv(ji,jj) = e2v(ji,jj) * vn(ji,jj,jk)
146            END DO 
147         END DO 
148 
149         ! East-West boundary conditions
150         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN
151            zwv(  0  ,:) = zwv(jpi-2,:)
152            zwv( -1  ,:) = zwv(jpi-3,:)
153            zwv(jpi+1,:) = zwv(  3  ,:)
154            zwv(jpi+2,:) = zwv(  4  ,:)
155         ELSE
156            zwv(  0  ,:) = 0.e0
157            zwv( -1  ,:) = 0.e0
158            zwv(jpi+1,:) = 0.e0
159            zwv(jpi+2,:) = 0.e0
160         ENDIF
161
162         ! North-South boundary conditions
163         IF( nperio == 3 .OR. nperio == 4 ) THEN
164            ! north fold ( Grid defined with a T-point pivot) ORCA 2 degre
165            zwu(jpi,jpj+1) = 0.e0
166            zwu(jpi,jpj+2) = 0.e0
167            DO ji = 1, jpi-1
168               iju = jpi - ji + 1
169               zwu(ji,jpj+1) = - zwu(iju,jpj-3)
170               zwu(ji,jpj+2) = - zwu(iju,jpj-4)
171            END DO
172         ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN
173            ! north fold ( Grid defined with a F-point pivot) ORCA 0.5 degre\
174            zwu(jpi,jpj+1) = 0.e0
175            zwu(jpi,jpj+2) = 0.e0
176            DO ji = 1, jpi-1
177               iju = jpi - ji
178               zwu(ji,jpj  ) = - zwu(iju,jpj-1)
179               zwu(ji,jpj+1) = - zwu(iju,jpj-2)
180               zwu(ji,jpj+2) = - zwu(iju,jpj-3)
181            END DO
182            DO ji = -1, jpi+2
183               ijt = jpi - ji + 1
184               zwv(ji,jpj) = - zwv(ijt,jpj-2)
185            END DO
186            DO ji = jpi/2+1, jpi+2
187               ijt = jpi - ji + 1
188               zwv(ji,jpjm1) = - zwv(ijt,jpjm1)
189            END DO
190         ELSE
191            ! closed
192            zwu(:,jpj+1) = 0.e0
193            zwu(:,jpj+2) = 0.e0
194         ENDIF
195
196         ! relative vorticity (vertical component of the velocity curl)
197         DO jj = 1, jpjm1
198            DO ji = 1, fs_jpim1   ! vector opt.
199               rotn(ji,jj,jk) = (  zwv(ji+1,jj  ) - zwv(ji,jj)      &
200                  &              - zwu(ji  ,jj+1) + zwu(ji,jj)  ) * fmask(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
201            END DO
202         END DO
203
204         ! second order accurate scheme along straight coast
205         DO jl = 1, npcoa(1,jk)
206            ii = nicoa(jl,1,jk)
207            ij = njcoa(jl,1,jk)
208            rotn(ii,ij,jk) = 1. / ( e1f(ii,ij) * e2f(ii,ij) )   &
209                           * ( + 4. * zwv(ii+1,ij) - zwv(ii+2,ij) + 0.2 * zwv(ii+3,ij) )
210         END DO
211         DO jl = 1, npcoa(2,jk)
212            ii = nicoa(jl,2,jk)
213            ij = njcoa(jl,2,jk)
214            rotn(ii,ij,jk) = 1./(e1f(ii,ij)*e2f(ii,ij))   &
215               *(-4.*zwv(ii,ij)+zwv(ii-1,ij)-0.2*zwv(ii-2,ij))
216         END DO
217         DO jl = 1, npcoa(3,jk)
218            ii = nicoa(jl,3,jk)
219            ij = njcoa(jl,3,jk)
220            rotn(ii,ij,jk) = -1. / ( e1f(ii,ij)*e2f(ii,ij) )   &
221               * ( +4. * zwu(ii,ij+1) - zwu(ii,ij+2) + 0.2 * zwu(ii,ij+3) )
222         END DO
223         DO jl = 1, npcoa(4,jk)
224            ii = nicoa(jl,4,jk)
225            ij = njcoa(jl,4,jk)
226            rotn(ii,ij,jk) = -1. / ( e1f(ii,ij)*e2f(ii,ij) )   &
227               * ( -4. * zwu(ii,ij) + zwu(ii,ij-1) - 0.2 * zwu(ii,ij-2) )
228         END DO
229         !                                             ! ===============
230      END DO                                           !   End of slab
231      !                                                ! ===============
232
233      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs   (update hdivn field)
234      IF( ln_divisf .AND. (nn_isf /= 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field)
235      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (Update Hor. divergence)
236     
237      ! 4. Lateral boundary conditions on hdivn and rotn
238      ! ---------------------------------=======---======
239      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )    ! lateral boundary cond. (no sign change)
240      !
241      CALL wrk_dealloc( jpi  , jpj+2, zwu )
242      CALL wrk_dealloc( jpi+2, jpj  , zwv )
243      !
244      IF( nn_timing == 1 )  CALL timing_stop('div_cur')
245      !
246   END SUBROUTINE div_cur
247   
248#else
249   !!----------------------------------------------------------------------
250   !!   Default option                           2nd order centered schemes
251   !!----------------------------------------------------------------------
252
253   SUBROUTINE div_cur( kt )
254      !!----------------------------------------------------------------------
255      !!                  ***  ROUTINE div_cur  ***
256      !!                   
257      !! ** Purpose :   compute the horizontal divergence and the relative
258      !!      vorticity at before and now time-step
259      !!
260      !! ** Method  : - Divergence:
261      !!      - save the divergence computed at the previous time-step
262      !!      (note that the Asselin filter has not been applied on hdivb)
263      !!      - compute the now divergence given by :
264      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
265      !!      correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla)
266      !!              - Relavtive Vorticity :
267      !!      - save the curl computed at the previous time-step (rotb = rotn)
268      !!      (note that the Asselin time filter has not been applied to rotb)
269      !!      - compute the now curl in tensorial formalism:
270      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
271      !!      Note: Coastal boundary condition: lateral friction set through
272      !!      the value of fmask along the coast (see dommsk.F90) and shlat
273      !!      (namelist parameter)
274      !!
275      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
276      !!              - update rotb , rotn , the before & now rel. vorticity
277      !!----------------------------------------------------------------------
278      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
279      !
280      INTEGER  ::   ji, jj, jk    ! dummy loop indices
281      REAL(wp) ::   zraur, zdep   ! local scalars
282      !!----------------------------------------------------------------------
283      !
284      IF( nn_timing == 1 )  CALL timing_start('div_cur')
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      !                                                ! ===============
293      DO jk = 1, jpkm1                                 ! Horizontal slab
294         !                                             ! ===============
295         !
296         hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
297         rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
298         !
299         !                                             ! --------
300         ! Horizontal divergence                       !   div
301         !                                             ! --------
302         DO jj = 2, jpjm1
303            DO ji = fs_2, fs_jpim1   ! vector opt.
304               hdivn(ji,jj,jk) =   &
305                  (  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)       &
306                   + 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)  )    &
307                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
308            END DO 
309         END DO 
310
311         IF( .NOT. AGRIF_Root() ) THEN
312            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
313            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
314            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
315            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
316         ENDIF
317
318         !                                             ! --------
319         ! relative vorticity                          !   rot
320         !                                             ! --------
321         DO jj = 1, jpjm1
322            DO ji = 1, fs_jpim1   ! vector opt.
323               rotn(ji,jj,jk) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
324                  &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
325                  &           * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
326            END DO
327         END DO
328         !                                             ! ===============
329      END DO                                           !   End of slab
330      !                                                ! ===============
331
332      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )                            ! runoffs (update hdivn field)
333      IF( ln_divisf .AND. (nn_isf .GT. 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field)
334      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
335      !
336      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )     ! lateral boundary cond. (no sign change)
337      !
338      IF( nn_timing == 1 )  CALL timing_stop('div_cur')
339      !
340   END SUBROUTINE div_cur
341   
342#endif
343
344
345   ! TODO - remove kt only used for validation
346   SUBROUTINE dyn_vrt_dia_3d( utend, vtend, id_dia_vor, kt)
347
348      !!----------------------------------------------------------------------
349      !!
350      !! ** Purpose :  compute the vertical integrals of utend and vtend, and
351      !!               then pass to dyn_vrt_dia_2d to calculate vorticity
352      !!               tendencies.
353      !!
354      !! ** Action : a) Calculate the vertical integrals of utend & of vtend
355      !!                (u_int & v_int)
356      !!             b) Call dyn_vrt_dia_2d with vertical integrals
357      !!
358      !!----------------------------------------------------------------------
359      REAL             :: utend(jpi,jpj,jpk) ! contribution to du/dt
360      REAL             :: vtend(jpi,jpj,jpk) ! contribution to dv/dt
361      CHARACTER(len=3) :: id_dia_vor         ! identifier for the diagnostic
362      INTEGER          :: kt                 ! ocean time-step index TODO remove after validation
363      !
364      !!----------------------------------------------------------------------
365      !
366      INTEGER  ::   ji, jj, jk     ! dummy loop indices
367      !
368      REAL(wp), POINTER, DIMENSION(:,:) :: u_int   ! u vertical integral
369      REAL(wp), POINTER, DIMENSION(:,:) :: v_int   ! v vertical integral
370
371      CALL wrk_alloc(jpi, jpj, u_int)
372      CALL wrk_alloc(jpi, jpj, v_int)
373
374      u_int(:,:) = 0.0_wp
375      v_int(:,:) = 0.0_wp
376
377      !
378      ! Calculate the vertical integrals of utend & of vtend
379      !
380      ! TODO remove - for validation only
381      IF ( kt == 1 ) THEN
382          WRITE(numout, *) 'dyn_vrt_dia_3d id:', id_dia_vor,       &
383                 &         ':bathy:', bathy(17,12),                &
384                 &         ':mbathy:', mbathy(17,12)
385      END IF
386     
387      DO jk = 1, jpkm1
388          DO jj = 2, jpjm1
389              DO ji = fs_2, fs_jpim1
390                  u_int(ji,jj) = u_int(ji,jj) + ( utend(ji,jj,jk) * fse3u(ji,jj,jk) &
391                                 &                 * e1u(ji,jj) * umask(ji,jj,jk) )
392                  v_int(ji,jj) = v_int(ji,jj) + ( vtend(ji,jj,jk) * fse3v(ji,jj,jk) &
393                                 &                 * e2v(ji,jj) * vmask(ji,jj,jk) )
394
395                  ! TODO remove - for validation only
396                  IF ( kt == 1 .AND. ji == 17 .AND. jj == 12 ) THEN
397                      WRITE(numout, *) 'dyn_vrt_dia_3d id:', id_dia_vor,       &
398                             &         ' :ji: ',     ji,                       &
399                             &         ' :jj: ',     jj,                       &
400                             &         ' :jk: ',     jk,                       &
401                             &         ' :u_int:',   u_int(ji,jj),             &
402                             &         ' :u_tend: ', utend(ji,jj,jk),          &
403                             &         ' :e1u: ',    e1u(ji,jj),               &
404                             &         ' :umask: ',  umask(ji,jj,jk),          &
405                             &         ' :fse3u: ',  fse3u(ji,jj,jk)
406                      WRITE(numout, *) 'dyn_vrt_dia_3d id:', id_dia_vor,       &
407                             &         ' :ji: ',     ji,                       &
408                             &         ' :jj: ',     jj,                       &
409                             &         ' :jk: ',     jk,                       &
410                             &         ' :v_int:',   v_int(ji,jj),             &
411                             &         ' :v_tend: ', vtend(ji,jj,jk),          &
412                             &         ' :e2v: ',    e2v(ji,jj),               &
413                             &         ' :vmask: ',  vmask(ji,jj,jk),          &
414                             &         ' :fse3v: ',  fse3v(ji,jj,jk)
415                  END IF
416              END DO
417          END DO
418      END DO
419
420      CALL dyn_vrt_dia_2d(u_int, v_int, id_dia_vor, kt)
421
422      CALL wrk_dealloc(jpi, jpj, u_int)
423      CALL wrk_dealloc(jpi, jpj, v_int)
424
425   END SUBROUTINE dyn_vrt_dia_3d
426
427
428   ! TODO - remove kt only used for validation
429   SUBROUTINE dyn_vrt_dia_2d( u_int, v_int, id_dia_vor, kt)
430
431      !!----------------------------------------------------------------------
432      !!
433      !! ** Purpose :  compute the integral and mean vorticity tendencies.
434      !!
435      !! ** Action : a) Calculate the vorticity tendencies for the vertical
436      !!                integrals.
437      !!             b) Calculate the vertical means, u_mn, v_mn from u_int
438      !!                and v_int by dividing by the depth
439      !!             c) Calculate the vorticity tendencies for the vertical
440      !!                means
441      !!             d) Call iom_put for the vertical integral vorticity
442      !!                tendencies (using cid_dia_vor_int)
443      !!             e) Call iom_put for the vertical mean vorticity
444      !!                tendencies (using cid_dia_vor_mn)
445      !!
446      !!----------------------------------------------------------------------
447      REAL             :: u_int(jpi,jpj)  ! u vertical integral
448      REAL             :: v_int(jpi,jpj)  ! v vertical integral
449      CHARACTER(len=3) :: id_dia_vor      ! identifier for the vorticity diagnostic
450      INTEGER          :: kt              ! ocean time-step index TODO remove after validation
451      !
452      !!----------------------------------------------------------------------
453      !
454      INTEGER  ::   ji, jj       ! dummy loop indices
455      INTEGER  ::   ikbu, ikbv   ! dummy loop indices
456      !
457      REAL(wp), POINTER, DIMENSION(:,:) :: u_mn    ! u vertical means
458      REAL(wp), POINTER, DIMENSION(:,:) :: v_mn    ! u vertical means
459      REAL(wp), POINTER, DIMENSION(:,:) :: vor_int ! vort trend of vert integrals
460      REAL(wp), POINTER, DIMENSION(:,:) :: vor_mn  ! vort trend of vert means
461      CHARACTER(len=20) ::   cid_dia_vor_int       ! Identifier for iom_put
462      CHARACTER(len=20) ::   cid_dia_vor_mn        ! Identifier for iom_put
463
464      CALL wrk_alloc(jpi, jpj, u_mn)
465      CALL wrk_alloc(jpi, jpj, v_mn)
466      CALL wrk_alloc(jpi, jpj, vor_int)
467      CALL wrk_alloc(jpi, jpj, vor_mn)
468
469      CALL lbc_lnk( u_int, 'U', 1. )
470      CALL lbc_lnk( v_int, 'V', 1. )
471
472      WRITE ( cid_dia_vor_int, "(A16,A3)" ) "dia_vor_int-",  id_dia_vor
473      WRITE ( cid_dia_vor_mn,  "(A17,A3)" ) "dia_vor_mean-", id_dia_vor
474
475      !
476      ! Calculate the vorticity tendencies for the vertical integrals.
477      ! 1/e1e2 * ((d(vtend)/dx) - (d(utend)/dy))
478      !
479
480      DO jj = 2, jpjm1
481          DO ji = fs_2, fs_jpim1
482              vor_int(ji,jj) = (   ( v_int(ji+1,jj) - v_int(ji,jj) )     &
483                  &              - ( u_int(ji,jj+1) - u_int(ji,jj) ) )   &
484                  &            / ( e1f(ji,jj)    * e2f(ji,jj) )
485
486              ! TODO remove - for validation only
487              IF ( kt == 1 .AND. ji == 17 .AND. jj == 12 ) THEN
488                  WRITE(numout, *) 'dyn_vrt_dia_2d id:', id_dia_vor,     &
489                       &           ':ji:', ji,                           &
490                       &           ':jj:', jj,                           &
491                       &           ':vor_int:', vor_int(ji,jj),          &
492                       &           ':v_int(i+1):', v_int(ji+1,jj),       &
493                       &           ':v_int:', v_int(ji,jj),              &
494                       &           ':u_int:', u_int(ji,jj),              &
495                       &           ':u_int(j+1):', u_int(ji,jj+1),       &
496                       &           ':e1f:', e1f(ji,jj),                  &
497                       &           ':e2f:', e2f(ji,jj)
498              END IF
499          END DO
500      END DO
501
502      ! Multiply by the surface mask
503      vor_int(:,:) = vor_int(:,:) * fmask(:,:,1)
504
505      !
506      ! Calculate the vertical means, u_mn, v_mn from u_int & v_int by dividing
507      ! by the depth
508      ! mbku & mbkv - vertical index of the bottom last U- & W- ocean level
509      !
510
511      DO jj = 2, jpjm1
512          DO ji = fs_2, fs_jpim1
513              ikbu = mbku(ji,jj)
514              ikbv = mbkv(ji,jj)
515
516              IF (ikbu .ne. 0.0_wp) THEN      ! Don't divide by 0!
517                  u_mn(ji,jj) = u_int(ji,jj) / gdepw_n(ji,jj,ikbu+1)
518              ELSE
519                  u_mn(ji,jj) = 0.0_wp
520              END IF
521
522              IF (ikbv .ne. 0.0_wp) THEN      ! Don't divide by 0!
523                  v_mn(ji,jj) = v_int(ji,jj) / gdepw_n(ji,jj,ikbv+1)
524              ELSE
525                  v_mn(ji,jj) = 0.0_wp
526              END IF
527
528              ! TODO remove - for validation only
529              IF ( kt == 1 .AND. ji == 17 .AND. jj == 12 ) THEN
530                  WRITE(numout, *) 'dyn_vrt_dia_2d id:', id_dia_vor,          &
531                       &           ':gdepw_n(ikbu):', gdepw_n(ji,jj,ikbu+1),  &
532                       &           ':u_int:', u_int(ji,jj),                   &
533                       &           ':u_mn:', u_mn(ji,jj),                     &
534                       &           ':gdepw_n(ikbv):', gdepw_n(ji,jj,ikbv+1),  &
535                       &           ':v_int:', v_int(ji,jj),                   &
536                       &           ':v_mn:', v_mn(ji,jj)
537              END IF
538          END DO
539      END DO
540
541      !
542      ! Calculate the vorticity tendencies for the vertical means
543      ! 1/e1e2 * ((d(v_mn)/dx) - (d(u_mn)/dy))
544      !
545
546      DO jj = 2, jpjm1
547          DO ji = fs_2, fs_jpim1
548              vor_mn(ji,jj) = (   ( v_mn(ji+1,jj) - v_mn(ji,jj) )        &
549                  &             - ( u_mn(ji,jj+1) - u_mn(ji,jj) ) )      &
550                  &          / ( e1f(ji,jj)   * e2f(ji,jj) )
551
552              ! TODO remove - for validation only
553              IF ( kt == 1 .AND. ji == 17 .AND. jj == 12 ) THEN
554                  WRITE(numout, *) 'dyn_vrt_dia_2d id:', id_dia_vor,     &
555                       &           ':ji:', ji,                           &
556                       &           ':jj:', jj,                           &
557                       &           ':vor_mn:', vor_mn(ji,jj),            &
558                       &           ':v_mn(i+1):', v_mn(ji+1,jj),         &
559                       &           ':v_mn:', v_mn(ji,jj),                &
560                       &           ':u_mn:', u_mn(ji,jj),                &
561                       &           ':u_mn(j+1):', u_mn(ji,jj+1),         &
562                       &           ':e1f:', e1f(ji,jj),                  &
563                       &           ':e2f:', e2f(ji,jj)
564              END IF
565          END DO
566      END DO
567
568      ! Multiply by the surface mask
569      vor_mn(:,:) = vor_mn(:,:) * fmask(:,:,1)
570
571
572      ! Call iom_put for the vertical integral vorticity tendencies
573      CALL iom_put( cid_dia_vor_int, vor_int(:,:))
574
575      ! Call iom_put for the vertical mean vorticity tendencies
576      CALL iom_put( cid_dia_vor_mn, vor_mn(:,:))
577
578      CALL wrk_dealloc(jpi, jpj, u_mn)
579      CALL wrk_dealloc(jpi, jpj, v_mn)
580      CALL wrk_dealloc(jpi, jpj, vor_int)
581      CALL wrk_dealloc(jpi, jpj, vor_mn)
582
583   END SUBROUTINE dyn_vrt_dia_2d
584
585
586   !!======================================================================
587END MODULE divcur
Note: See TracBrowser for help on using the repository browser.