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.
traadv_eiv.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90 @ 7061

Last change on this file since 7061 was 7061, checked in by davestorkey, 8 years ago

Enable 3D tracer trends diagnostics ,ie. diagnostics obtained for ln_tra_trd=T, in UKMO/dev_r5518_GO6_package branch.
Remove error check for VVL and implement thickness weighting.
Also implement new diagnostics:
ttrd_totad : the total trend due to advection (correct alternative to ttrd_ad!)
ttrd_iso_x/y/z1 : individual components of trend due to isopycnal diffusion. Note the "z1" diagnostic is only part of the vertical trend. The other part is ttrd_zdf minus ttrd_zdfp.
ttrd_evd : the trend due to the EVD convection term (which is contained in ttrd_zdf and ttrd_zdfp).
ttrd_tot : the total model trend.
u/v/weiv_masstr3d : the components of the mass transport due to the Gent-McWilliams? scheme.
u/v/weiv_heattr3d : the components of the heat transport due to the Gent-McWilliams? scheme.
u/v/weiv_salttr3d : the components of the salt transport due to the Gent-McWilliams? scheme.
File size: 15.5 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
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   tra_adv_eiv   ! routine called by step.F90
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "ldftra_substitute.h90"
40#  include "ldfeiv_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE tra_adv_eiv( kt, kit000, pun, pvn, pwn, cdtype )
50      !!----------------------------------------------------------------------
51      !!                  ***  ROUTINE tra_adv_eiv  ***
52      !!
53      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
54      !!      trend and add it to the general trend of tracer equation.
55      !!
56      !! ** Method  :   The eddy induced advection is computed from the slope
57      !!      of iso-neutral surfaces computed in routine ldf_slp as follows:
58      !!         zu_eiv =  1/(e2u e3u)   dk[ aeiu e2u mi(wslpi) ]
59      !!         zv_eiv =  1/(e1v e3v)   dk[ aeiv e1v mj(wslpj)
60      !!         zw_eiv = -1/(e1t e2t) { di[ aeiu e2u mi(wslpi) ]
61      !!                               + dj[ aeiv e1v mj(wslpj) ] }
62      !!      add the eiv component to the model velocity:
63      !!         p.n = p.n + z._eiv
64      !!
65      !! ** Action  : - add to p.n the eiv component
66      !!----------------------------------------------------------------------
67      INTEGER                         , INTENT(in   ) ::   kt       ! ocean time-step index
68      INTEGER                         , INTENT(in   ) ::   kit000   ! first time step index
69      CHARACTER(len=3)                , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
70      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun      ! in : 3 ocean velocity components
71      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pvn      ! out: 3 ocean velocity components
72      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pwn      ! increased by the eiv
73      !!
74      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
75      REAL(wp) ::   zuwk, zuwk1, zuwi, zuwi1   ! local scalars
76      REAL(wp) ::   zvwk, zvwk1, zvwj, zvwj1   !   -      -
77# if defined key_diaeiv 
78      REAL(wp) ::   zztmp                      ! local scalar
79# endif 
80      REAL(wp), POINTER, DIMENSION(:,:) :: zu_eiv, zv_eiv, zw_eiv, z2d
81      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, z3d_T
82      !!----------------------------------------------------------------------
83      !
84      IF( nn_timing == 1 )  CALL timing_start( 'tra_adv_eiv')
85      !
86# if defined key_diaeiv 
87      CALL wrk_alloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv, z2d )
88      CALL wrk_alloc( jpi, jpj, jpk, z3d, z3d_T )
89# else
90      CALL wrk_alloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv )
91# endif
92
93      IF( kt == kit000 )  THEN
94         IF(lwp) WRITE(numout,*)
95         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection on ', cdtype,' :'
96         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component'
97# if defined key_diaeiv 
98         IF( cdtype == 'TRA') THEN
99            u_eiv(:,:,:) = 0.e0
100            v_eiv(:,:,:) = 0.e0
101            w_eiv(:,:,:) = 0.e0
102         END IF
103# endif
104      ENDIF
105
106      zu_eiv(:,:) = 0.e0   ;   zv_eiv(:,:) = 0.e0   ;    zw_eiv(:,:) = 0.e0 
107     
108                                                    ! =================
109      DO jk = 1, jpkm1                              !  Horizontal slab
110         !                                          ! =================
111         DO jj = 1, jpjm1
112            DO ji = 1, fs_jpim1   ! vector opt.
113               zuwk = ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk  ) ) * fsaeiu(ji,jj,jk  ) * umask(ji,jj,jk  )
114               zuwk1= ( wslpi(ji,jj,jk+1) + wslpi(ji+1,jj,jk+1) ) * fsaeiu(ji,jj,jk+1) * umask(ji,jj,jk+1)
115               zvwk = ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk  ) ) * fsaeiv(ji,jj,jk  ) * vmask(ji,jj,jk  )
116               zvwk1= ( wslpj(ji,jj,jk+1) + wslpj(ji,jj+1,jk+1) ) * fsaeiv(ji,jj,jk+1) * vmask(ji,jj,jk+1)
117
118               zu_eiv(ji,jj) = 0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 ) 
119               zv_eiv(ji,jj) = 0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 ) 
120   
121               pun(ji,jj,jk) = pun(ji,jj,jk) + e2u(ji,jj) * zu_eiv(ji,jj)
122               pvn(ji,jj,jk) = pvn(ji,jj,jk) + e1v(ji,jj) * zv_eiv(ji,jj)
123            END DO
124         END DO
125# if defined key_diaeiv 
126         IF( cdtype == 'TRA') THEN
127            u_eiv(:,:,jk) = zu_eiv(:,:) / fse3u(:,:,jk)
128            v_eiv(:,:,jk) = zv_eiv(:,:) / fse3v(:,:,jk)
129         END IF
130# endif
131         IF( jk >=2 ) THEN                             ! jk=1 zw_eiv=0, not computed
132            DO jj = 2, jpjm1
133               DO ji = fs_2, fs_jpim1   ! vector opt.
134# if defined key_traldf_c2d || defined key_traldf_c3d
135                  zuwi  = ( wslpi(ji,jj,jk)+wslpi(ji-1,jj,jk) ) * fsaeiu(ji-1,jj,jk) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
136                  zuwi1 = ( wslpi(ji,jj,jk)+wslpi(ji+1,jj,jk) ) * fsaeiu(ji  ,jj,jk) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
137                  zvwj  = ( wslpj(ji,jj,jk)+wslpj(ji,jj-1,jk) ) * fsaeiv(ji,jj-1,jk) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
138                  zvwj1 = ( wslpj(ji,jj,jk)+wslpj(ji,jj+1,jk) ) * fsaeiv(ji,jj  ,jk) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
139 
140                  zw_eiv(ji,jj) = - 0.5 * tmask(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj ) 
141# else
142                  zuwi  = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) ) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
143                  zuwi1 = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) ) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
144                  zvwj  = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) ) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
145                  zvwj1 = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) ) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
146
147                  zw_eiv(ji,jj) = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj )
148# endif
149                  pwn(ji,jj,jk) = pwn(ji,jj,jk) + zw_eiv(ji,jj)
150               END DO
151            END DO
152# if defined key_diaeiv 
153            IF( cdtype == 'TRA')  w_eiv(:,:,jk) = zw_eiv(:,:) / ( e1t(:,:) * e2t(:,:) )
154# endif
155         ENDIF
156         !                                          ! =================
157      END DO                                        !    End of slab 
158      !                                             ! =================
159
160# if defined key_diaeiv 
161      IF( cdtype == 'TRA') THEN
162         CALL iom_put( "uoce_eiv", u_eiv )    ! i-eiv current
163         CALL iom_put( "voce_eiv", v_eiv )    ! j-eiv current
164         CALL iom_put( "woce_eiv", w_eiv )    ! vert. eiv current
165
166         IF( iom_use("ueiv_masstr") .OR. iom_use("ueiv_heattr") .OR. iom_use('ueiv_heattr3d')        &
167                                    .OR. iom_use("ueiv_salttr") .OR. iom_use('ueiv_salttr3d') ) THEN
168            z3d(:,:,jpk) = 0.e0
169            DO jk = 1, jpkm1
170               z3d(:,:,jk) = rau0 * u_eiv(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)
171            END DO
172            CALL iom_put( "ueiv_masstr", z3d )                  ! mass transport in i-direction
173         ENDIF
174
175         IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN
176            zztmp = 0.5 * rcp 
177            z2d(:,:) = 0.e0 
178            z3d_T(:,:,:) = 0.e0 
179            DO jk = 1, jpkm1
180               DO jj = 2, jpjm1
181                  DO ji = fs_2, fs_jpim1   ! vector opt.
182                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
183                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk) 
184                  END DO
185               END DO
186            END DO
187            IF (iom_use('ueiv_heattr') ) THEN
188               CALL lbc_lnk( z2d, 'U', -1. )
189               CALL iom_put( "ueiv_heattr", zztmp * z2d )                  ! 2D heat transport in i-direction
190            ENDIF
191            IF (iom_use('ueiv_heattr3d') ) THEN
192               CALL lbc_lnk( z3d_T, 'U', -1. )
193               CALL iom_put( "ueiv_heattr3d", zztmp * z3d_T )              ! 3D heat transport in i-direction
194            ENDIF
195         ENDIF
196
197         IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d') ) THEN
198            zztmp = 0.5 * 0.001
199            z2d(:,:) = 0.e0 
200            z3d_T(:,:,:) = 0.e0 
201            DO jk = 1, jpkm1
202               DO jj = 2, jpjm1
203                  DO ji = fs_2, fs_jpim1   ! vector opt.
204                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
205                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk) 
206                  END DO
207               END DO
208            END DO
209            IF (iom_use('ueiv_salttr') ) THEN
210               CALL lbc_lnk( z2d, 'U', -1. )
211               CALL iom_put( "ueiv_salttr", zztmp * z2d )                  ! 2D salt transport in i-direction
212            ENDIF
213            IF (iom_use('ueiv_salttr3d') ) THEN
214               CALL lbc_lnk( z3d_T, 'U', -1. )
215               CALL iom_put( "ueiv_salttr3d", zztmp * z3d_T )              ! 3D salt transport in i-direction
216            ENDIF
217         ENDIF
218
219         IF( iom_use("veiv_masstr") .OR. iom_use("veiv_heattr") .OR. iom_use('veiv_heattr3d')       &
220                                    .OR. iom_use("veiv_salttr") .OR. iom_use('veiv_salttr3d') ) THEN
221            z3d(:,:,jpk) = 0.e0
222            DO jk = 1, jpkm1
223               z3d(:,:,jk) = rau0 * v_eiv(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
224            END DO
225            CALL iom_put( "veiv_masstr", z3d )                  ! mass transport in j-direction
226         ENDIF
227           
228         IF( iom_use('veiv_heattr') .OR. iom_use('veiv_heattr3d') ) THEN
229            zztmp = 0.5 * rcp 
230            z2d(:,:) = 0.e0 
231            z3d_T(:,:,:) = 0.e0 
232            DO jk = 1, jpkm1
233               DO jj = 2, jpjm1
234                  DO ji = fs_2, fs_jpim1   ! vector opt.
235                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
236                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk) 
237                  END DO
238               END DO
239            END DO
240            IF (iom_use('veiv_heattr') ) THEN
241               CALL lbc_lnk( z2d, 'V', -1. )
242               CALL iom_put( "veiv_heattr", zztmp * z2d )                  ! 2D heat transport in j-direction
243            ENDIF
244            IF (iom_use('veiv_heattr3d') ) THEN
245               CALL lbc_lnk( z3d_T, 'V', -1. )
246               CALL iom_put( "veiv_heattr3d", zztmp * z3d_T )              ! 3D heat transport in j-direction
247            ENDIF
248         ENDIF
249
250         IF( iom_use('veiv_salttr') .OR. iom_use('veiv_salttr3d') ) THEN
251            zztmp = 0.5 * 0.001
252            z2d(:,:) = 0.e0 
253            z3d_T(:,:,:) = 0.e0 
254            DO jk = 1, jpkm1
255               DO jj = 2, jpjm1
256                  DO ji = fs_2, fs_jpim1   ! vector opt.
257                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
258                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk)
259                  END DO
260               END DO
261            END DO
262            IF (iom_use('veiv_salttr') ) THEN
263               CALL lbc_lnk( z2d, 'V', -1. )
264               CALL iom_put( "veiv_salttr", zztmp * z2d )                  ! 2D salt transport in i-direction
265            ENDIF
266            IF (iom_use('veiv_salttr3d') ) THEN
267               CALL lbc_lnk( z3d_T, 'V', -1. )
268               CALL iom_put( "veiv_salttr3d", zztmp * z3d_T )              ! 3D salt transport in i-direction
269            ENDIF
270         ENDIF
271
272         IF( iom_use('weiv_masstr') .OR. iom_use('weiv_heattr3d') .OR. iom_use('weiv_salttr3d')) THEN   ! vertical mass transport & its square value
273           z2d(:,:) = rau0 * e12t(:,:)
274           DO jk = 1, jpk
275              z3d(:,:,jk) = w_eiv(:,:,jk) * z2d(:,:)
276           END DO
277           CALL iom_put( "weiv_masstr" , z3d )                  ! mass transport in k-direction
278         ENDIF
279
280         IF( iom_use('weiv_heattr3d') ) THEN
281            zztmp = 0.5 * rcp 
282            DO jk = 1, jpkm1
283               DO jj = 2, jpjm1
284                  DO ji = fs_2, fs_jpim1   ! vector opt.
285                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj,jk+1,jp_tem) )
286                  END DO
287               END DO
288            END DO
289            CALL lbc_lnk( z3d_T, 'T', 1. )
290            CALL iom_put( "weiv_heattr3d", zztmp * z3d_T )                 ! 3D heat transport in k-direction
291         ENDIF
292
293         IF( iom_use('weiv_salttr3d') ) THEN
294            zztmp = 0.5 * 0.001 
295            DO jk = 1, jpkm1
296               DO jj = 2, jpjm1
297                  DO ji = fs_2, fs_jpim1   ! vector opt.
298                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj,jk+1,jp_sal) )
299                  END DO
300               END DO
301            END DO
302            CALL lbc_lnk( z3d_T, 'T', 1. )
303            CALL iom_put( "weiv_salttr3d", zztmp * z3d_T )                 ! 3D salt transport in k-direction
304         ENDIF
305
306    END IF
307# endif 
308
309# if defined key_diaeiv 
310      CALL wrk_dealloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv, z2d )
311      CALL wrk_dealloc( jpi, jpj, jpk, z3d, z3d_T )
312# else
313      CALL wrk_dealloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv )
314# endif
315      !
316      IF( nn_timing == 1 )  CALL timing_stop( 'tra_adv_eiv')
317      !
318    END SUBROUTINE tra_adv_eiv
319
320#else
321   !!----------------------------------------------------------------------
322   !!   Dummy module :             No rotation of the lateral mixing tensor
323   !!----------------------------------------------------------------------
324CONTAINS
325   SUBROUTINE tra_adv_eiv( kt, kit000, pun, pvn, pwn, cdtype )              ! Empty routine
326      INTEGER  ::   kt   
327      INTEGER  ::   kit000   
328      CHARACTER(len=3) ::   cdtype
329      REAL, DIMENSION(:,:,:) ::   pun, pvn, pwn
330      WRITE(*,*) 'tra_adv_eiv: You should not have seen this print! error?', &
331          &  kt, cdtype, pun(1,1,1), pvn(1,1,1), pwn(1,1,1)
332   END SUBROUTINE tra_adv_eiv
333#endif
334
335   !!==============================================================================
336END MODULE traadv_eiv
Note: See TracBrowser for help on using the repository browser.