source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90 @ 8627

Last change on this file since 8627 was 8627, checked in by gm, 3 years ago

#1963 : Correct implementation of CMIP6 KE dissipation + move EKE⇔PE diag in traadv_eiv (in both cases, no more need of ln_trd_KE=T)

  • Property svn:keywords set to Id
File size: 18.1 KB
Line 
1MODULE traadv_eiv
2   !!======================================================================
3   !!                    ***  MODULE  traadv_eiv  ***
4   !! Ocean tracers:  advection trend - eddy induced velocity
5   !!======================================================================
6   !! History :  1.0  !  2005-11 (G. Madec)  Original code, from traldf and zdf _iso
7   !!            3.3  !  2010-05 (C. Ethe, G. Madec)  merge TRC-TRA
8   !!----------------------------------------------------------------------
9#if defined key_traldf_eiv   ||   defined key_esopa
10   !!----------------------------------------------------------------------
11   !!   'key_traldf_eiv'                  rotation of the lateral mixing tensor
12   !!----------------------------------------------------------------------
13   !!   tra_ldf_iso : update the tracer trend with the horizontal component
14   !!                 of iso neutral laplacian operator or horizontal
15   !!                 laplacian operator in s-coordinate
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers variables
18   USE dom_oce         ! ocean space and time domain variables
19   USE ldftra_oce      ! ocean active tracers: lateral physics
20   USE ldfslp          ! iso-neutral slopes
21   USE in_out_manager  ! I/O manager
22   USE iom
23   USE trc_oce         ! share passive tracers/Ocean variables
24# if defined key_diaeiv
25   USE phycst          ! physical constants
26   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
27# endif 
28   USE wrk_nemo        ! Memory Allocation
29   USE timing          ! Timing
30   USE diaptr         ! Heat/Salt transport diagnostics
31   USE trddyn
32   USE trd_oce
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   tra_adv_eiv   ! routine called by step.F90
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "ldftra_substitute.h90"
42#  include "ldfeiv_substitute.h90"
43#  include "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
46   !! $Id$
47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   SUBROUTINE tra_adv_eiv( kt, kit000, pun, pvn, pwn, cdtype )
52      !!----------------------------------------------------------------------
53      !!                  ***  ROUTINE tra_adv_eiv  ***
54      !!
55      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
56      !!      trend and add it to the general trend of tracer equation.
57      !!
58      !! ** Method  :   The eddy induced advection is computed from the slope
59      !!      of iso-neutral surfaces computed in routine ldf_slp as follows:
60      !!         zu_eiv =  1/(e2u e3u)   dk[ aeiu e2u mi(wslpi) ]
61      !!         zv_eiv =  1/(e1v e3v)   dk[ aeiv e1v mj(wslpj)
62      !!         zw_eiv = -1/(e1t e2t) { di[ aeiu e2u mi(wslpi) ]
63      !!                               + dj[ aeiv e1v mj(wslpj) ] }
64      !!      add the eiv component to the model velocity:
65      !!         p.n = p.n + z._eiv
66      !!
67      !! ** Action  : - add to p.n the eiv component
68      !!----------------------------------------------------------------------
69      INTEGER                         , INTENT(in   ) ::   kt       ! ocean time-step index
70      INTEGER                         , INTENT(in   ) ::   kit000   ! first time step index
71      CHARACTER(len=3)                , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
72      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun      ! in : 3 ocean velocity components
73      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pvn      ! out: 3 ocean velocity components
74      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pwn      ! increased by the eiv
75      !!
76      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
77      REAL(wp) ::   zuwk, zuwk1, zuwi, zuwi1   ! local scalars
78      REAL(wp) ::   zvwk, zvwk1, zvwj, zvwj1   !   -      -
79# if defined key_diaeiv 
80      REAL(wp) ::   zztmp                      ! local scalar
81# endif 
82      REAL(wp), POINTER, DIMENSION(:,:) :: zu_eiv, zv_eiv, zw_eiv, z2d
83      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, z3d_T
84      !!----------------------------------------------------------------------
85      !
86      IF( nn_timing == 1 )  CALL timing_start( 'tra_adv_eiv')
87      !
88# if defined key_diaeiv 
89      CALL wrk_alloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv, z2d )
90      CALL wrk_alloc( jpi, jpj, jpk, z3d, z3d_T )
91# else
92      CALL wrk_alloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv )
93# endif
94
95      IF( kt == kit000 )  THEN
96         IF(lwp) WRITE(numout,*)
97         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection on ', cdtype,' :'
98         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component'
99# if defined key_diaeiv 
100         IF( cdtype == 'TRA') THEN
101            u_eiv(:,:,:) = 0.e0
102            v_eiv(:,:,:) = 0.e0
103            w_eiv(:,:,:) = 0.e0
104         END IF
105# endif
106      ENDIF
107
108      zu_eiv(:,:) = 0.e0   ;   zv_eiv(:,:) = 0.e0   ;    zw_eiv(:,:) = 0.e0 
109     
110                                                    ! =================
111      DO jk = 1, jpkm1                              !  Horizontal slab
112         !                                          ! =================
113         DO jj = 1, jpjm1
114            DO ji = 1, fs_jpim1   ! vector opt.
115               zuwk = ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk  ) ) * fsaeiu(ji,jj,jk  ) * umask(ji,jj,jk  )
116               zuwk1= ( wslpi(ji,jj,jk+1) + wslpi(ji+1,jj,jk+1) ) * fsaeiu(ji,jj,jk+1) * umask(ji,jj,jk+1)
117               zvwk = ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk  ) ) * fsaeiv(ji,jj,jk  ) * vmask(ji,jj,jk  )
118               zvwk1= ( wslpj(ji,jj,jk+1) + wslpj(ji,jj+1,jk+1) ) * fsaeiv(ji,jj,jk+1) * vmask(ji,jj,jk+1)
119
120               zu_eiv(ji,jj) = 0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 ) 
121               zv_eiv(ji,jj) = 0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 ) 
122   
123               pun(ji,jj,jk) = pun(ji,jj,jk) + e2u(ji,jj) * zu_eiv(ji,jj)
124               pvn(ji,jj,jk) = pvn(ji,jj,jk) + e1v(ji,jj) * zv_eiv(ji,jj)
125            END DO
126         END DO
127# if defined key_diaeiv 
128         IF( cdtype == 'TRA') THEN
129            u_eiv(:,:,jk) = zu_eiv(:,:) / fse3u(:,:,jk)
130            v_eiv(:,:,jk) = zv_eiv(:,:) / fse3v(:,:,jk)
131         END IF
132# endif
133         IF( jk >=2 ) THEN                             ! jk=1 zw_eiv=0, not computed
134            DO jj = 2, jpjm1
135               DO ji = fs_2, fs_jpim1   ! vector opt.
136# if defined key_traldf_c2d || defined key_traldf_c3d
137                  zuwi  = ( wslpi(ji,jj,jk)+wslpi(ji-1,jj,jk) ) * fsaeiu(ji-1,jj,jk) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
138                  zuwi1 = ( wslpi(ji,jj,jk)+wslpi(ji+1,jj,jk) ) * fsaeiu(ji  ,jj,jk) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
139                  zvwj  = ( wslpj(ji,jj,jk)+wslpj(ji,jj-1,jk) ) * fsaeiv(ji,jj-1,jk) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
140                  zvwj1 = ( wslpj(ji,jj,jk)+wslpj(ji,jj+1,jk) ) * fsaeiv(ji,jj  ,jk) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
141 
142                  zw_eiv(ji,jj) = - 0.5 * tmask(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj ) 
143# else
144                  zuwi  = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) ) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
145                  zuwi1 = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) ) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
146                  zvwj  = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) ) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
147                  zvwj1 = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) ) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
148
149                  zw_eiv(ji,jj) = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj )
150# endif
151                  pwn(ji,jj,jk) = pwn(ji,jj,jk) + zw_eiv(ji,jj)
152               END DO
153            END DO
154# if defined key_diaeiv 
155            IF( cdtype == 'TRA')  w_eiv(:,:,jk) = zw_eiv(:,:) / ( e1t(:,:) * e2t(:,:) )
156# endif
157         ENDIF
158         !                                          ! =================
159      END DO                                        !    End of slab 
160      !                                             ! =================
161
162# if defined key_diaeiv 
163      IF( cdtype == 'TRA') THEN
164         CALL iom_put( "uoce_eiv", u_eiv )    ! i-eiv current
165         CALL iom_put( "voce_eiv", v_eiv )    ! j-eiv current
166         CALL iom_put( "woce_eiv", w_eiv )    ! vert. eiv current
167         !
168         IF( iom_use('weiv_masstr') ) THEN   ! vertical mass transport & its square value
169            z2d(:,:) = rau0 * e12t(:,:)
170            DO jk = 1, jpk
171               z3d(:,:,jk) = w_eiv(:,:,jk) * z2d(:,:)
172            END DO
173            CALL iom_put( "weiv_masstr" , z3d ) 
174         ENDIF
175         !
176         IF( iom_use("ueiv_masstr") .OR. iom_use("ueiv_heattr") .OR. iom_use('ueiv_heattr3d')        &
177            &                       .OR. iom_use("ueiv_salttr") .OR. iom_use('ueiv_salttr3d') ) THEN
178            z3d(:,:,jpk) = 0.e0
179            z2d(:,:) = 0.e0
180            DO jk = 1, jpkm1
181               z3d(:,:,jk) = rau0 * u_eiv(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)
182               z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
183            END DO
184            CALL iom_put( "ueiv_masstr", z3d )                  ! mass transport in i-direction
185         ENDIF
186         !
187         IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN
188            zztmp = 0.5 * rcp 
189            z2d(:,:) = 0.e0 
190            z3d_T(:,:,:) = 0.e0 
191            DO jk = 1, jpkm1
192               DO jj = 2, jpjm1
193                  DO ji = fs_2, fs_jpim1   ! vector opt.
194                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
195                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk) 
196                  END DO
197               END DO
198            END DO
199            IF (iom_use('ueiv_heattr') ) THEN
200               CALL lbc_lnk( z2d, 'U', -1. )
201               CALL iom_put( "ueiv_heattr", zztmp * z2d )                  ! 2D heat transport in i-direction
202            ENDIF
203            IF (iom_use('ueiv_heattr3d') ) THEN
204               CALL lbc_lnk( z3d_T, 'U', -1. )
205               CALL iom_put( "ueiv_heattr3d", zztmp * z3d_T )              ! 3D heat transport in i-direction
206            ENDIF
207         ENDIF
208         !
209         IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d') ) THEN
210            zztmp = 0.5 * 0.001
211            z2d(:,:) = 0.e0 
212            z3d_T(:,:,:) = 0.e0 
213            DO jk = 1, jpkm1
214               DO jj = 2, jpjm1
215                  DO ji = fs_2, fs_jpim1   ! vector opt.
216                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
217                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk) 
218                  END DO
219               END DO
220            END DO
221            IF (iom_use('ueiv_salttr') ) THEN
222               CALL lbc_lnk( z2d, 'U', -1. )
223               CALL iom_put( "ueiv_salttr", zztmp * z2d )                  ! 2D salt transport in i-direction
224            ENDIF
225            IF (iom_use('ueiv_salttr3d') ) THEN
226               CALL lbc_lnk( z3d_T, 'U', -1. )
227               CALL iom_put( "ueiv_salttr3d", zztmp * z3d_T )              ! 3D salt transport in i-direction
228            ENDIF
229         ENDIF
230         !
231         IF( iom_use("veiv_masstr") .OR. iom_use("veiv_heattr") .OR. iom_use('veiv_heattr3d')       &
232                                    .OR. iom_use("veiv_salttr") .OR. iom_use('veiv_salttr3d') ) THEN
233            z3d(:,:,jpk) = 0.e0
234            DO jk = 1, jpkm1
235               z3d(:,:,jk) = rau0 * v_eiv(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
236            END DO
237            CALL iom_put( "veiv_masstr", z3d )                  ! mass transport in j-direction
238         ENDIF
239         !   
240         IF( iom_use('veiv_heattr') .OR. iom_use('veiv_heattr3d') ) THEN
241            zztmp = 0.5 * rcp 
242            z2d(:,:) = 0.e0 
243            z3d_T(:,:,:) = 0.e0 
244            DO jk = 1, jpkm1
245               DO jj = 2, jpjm1
246                  DO ji = fs_2, fs_jpim1   ! vector opt.
247                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
248                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk) 
249                  END DO
250               END DO
251            END DO
252            IF (iom_use('veiv_heattr') ) THEN
253               CALL lbc_lnk( z2d, 'V', -1. )
254               CALL iom_put( "veiv_heattr", zztmp * z2d )                  ! 2D heat transport in j-direction
255            ENDIF
256            IF (iom_use('veiv_heattr3d') ) THEN
257               CALL lbc_lnk( z3d_T, 'V', -1. )
258               CALL iom_put( "veiv_heattr3d", zztmp * z3d_T )              ! 3D heat transport in j-direction
259            ENDIF
260         ENDIF
261         !
262         IF( iom_use('veiv_salttr') .OR. iom_use('veiv_salttr3d') ) THEN
263            zztmp = 0.5 * 0.001
264            z2d(:,:) = 0.e0 
265            z3d_T(:,:,:) = 0.e0 
266            DO jk = 1, jpkm1
267               DO jj = 2, jpjm1
268                  DO ji = fs_2, fs_jpim1   ! vector opt.
269                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
270                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk)
271                  END DO
272               END DO
273            END DO
274            IF (iom_use('veiv_salttr') ) THEN
275               CALL lbc_lnk( z2d, 'V', -1. )
276               CALL iom_put( "veiv_salttr", zztmp * z2d )                  ! 2D salt transport in i-direction
277            ENDIF
278            IF (iom_use('veiv_salttr3d') ) THEN
279               CALL lbc_lnk( z3d_T, 'V', -1. )
280               CALL iom_put( "veiv_salttr3d", zztmp * z3d_T )              ! 3D salt transport in i-direction
281            ENDIF
282         ENDIF
283         !
284         IF( iom_use('weiv_masstr') .OR. iom_use('weiv_heattr3d') .OR. iom_use('weiv_salttr3d')) THEN   ! vertical mass transport & its square value
285           z2d(:,:) = rau0 * e12t(:,:)
286           DO jk = 1, jpk
287              z3d(:,:,jk) = w_eiv(:,:,jk) * z2d(:,:)
288           END DO
289           CALL iom_put( "weiv_masstr" , z3d )                  ! mass transport in k-direction
290         ENDIF
291         !
292         IF( iom_use('weiv_heattr3d') ) THEN
293            zztmp = 0.5 * rcp 
294            DO jk = 1, jpkm1
295               DO jj = 2, jpjm1
296                  DO ji = fs_2, fs_jpim1   ! vector opt.
297                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj,jk+1,jp_tem) )
298                  END DO
299               END DO
300            END DO
301            CALL lbc_lnk( z3d_T, 'T', 1. )
302            CALL iom_put( "weiv_heattr3d", zztmp * z3d_T )                 ! 3D heat transport in k-direction
303         ENDIF
304         !
305         IF( iom_use('weiv_salttr3d') ) THEN
306            zztmp = 0.5 * 0.001 
307            DO jk = 1, jpkm1
308               DO jj = 2, jpjm1
309                  DO ji = fs_2, fs_jpim1   ! vector opt.
310                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj,jk+1,jp_sal) )
311                  END DO
312               END DO
313            END DO
314            CALL lbc_lnk( z3d_T, 'T', 1. )
315            CALL iom_put( "weiv_salttr3d", zztmp * z3d_T )                 ! 3D salt transport in k-direction
316         ENDIF
317         !
318         IF( ln_diaptr ) THEN
319            z3d(:,:,:) = 0._wp
320            DO jk = 1, jpkm1
321               DO jj = 2, jpjm1
322                  DO ji = fs_2, fs_jpim1   ! vector opt.
323                     z3d(ji,jj,jk) = v_eiv(ji,jj,jk) * 0.5 * (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) &
324                        &             * e1v(ji,jj) * fse3v(ji,jj,jk)
325                  END DO
326               END DO
327            END DO
328            CALL dia_ptr_ohst_components( jp_tem, 'eiv', z3d )
329            z3d(:,:,:) = 0._wp
330            DO jk = 1, jpkm1
331               DO jj = 2, jpjm1
332                  DO ji = fs_2, fs_jpim1   ! vector opt.
333                     z3d(ji,jj,jk) = v_eiv(ji,jj,jk) * 0.5 * (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) &
334                     &             * e1v(ji,jj) * fse3v(ji,jj,jk)
335                  END DO
336               END DO
337            END DO
338            CALL dia_ptr_ohst_components( jp_sal, 'eiv', z3d )
339         ENDIF
340         !
341!!gm add CMIP6 diag here instead of been done in trdken.F90
342         !
343         IF( iom_use('eketrd_eiv') ) THEN     ! tendency of EKE from parameterized eddy advection
344            ! CMIP6 diagnostic tknebto = tendency of EKE from parameterized mesoscale eddy advection
345            ! = vertical_integral( k (N S)^2 ) rho dz   where rho = rau0 and S = isoneutral slope.
346            z2d(:,:) = 0._wp
347            DO jk = 1, jpkm1
348               DO ji = 1, jpi
349                  DO jj = 1,jpj
350                     z2d(ji,jj) = z2d(ji,jj) + rau0 * fsaeiw(ji,jj,jk)                 &
351                        &                    * rn2b(ji,jj,jk) * fse3w(ji,jj,jk)        &
352                        &                    * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)    &
353                        &                       + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  ) * wmask(ji,jj,jk)
354                  END DO
355               END DO
356            END DO
357            CALL iom_put( "eketrd_eiv", z2d )
358         ENDIF
359         !
360!!gm  removed from trdken.F90    IF( ln_KE_trd )   CALL trd_dyn(u_eiv, v_eiv, jpdyn_eivke, kt )
361         !
362      ENDIF
363# endif 
364
365# if defined key_diaeiv 
366      CALL wrk_dealloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv, z2d )
367      CALL wrk_dealloc( jpi, jpj, jpk, z3d, z3d_T )
368# else
369      CALL wrk_dealloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv )
370# endif
371      !
372      IF( nn_timing == 1 )  CALL timing_stop( 'tra_adv_eiv')
373      !
374    END SUBROUTINE tra_adv_eiv
375
376#else
377   !!----------------------------------------------------------------------
378   !!   Dummy module :             No rotation of the lateral mixing tensor
379   !!----------------------------------------------------------------------
380CONTAINS
381   SUBROUTINE tra_adv_eiv( kt, kit000, pun, pvn, pwn, cdtype )              ! Empty routine
382      INTEGER  ::   kt   
383      INTEGER  ::   kit000   
384      CHARACTER(len=3) ::   cdtype
385      REAL, DIMENSION(:,:,:) ::   pun, pvn, pwn
386      WRITE(*,*) 'tra_adv_eiv: You should not have seen this print! error?', &
387          &  kt, cdtype, pun(1,1,1), pvn(1,1,1), pwn(1,1,1)
388   END SUBROUTINE tra_adv_eiv
389#endif
390
391   !!==============================================================================
392END MODULE traadv_eiv
Note: See TracBrowser for help on using the repository browser.