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

Last change on this file was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.
File size: 16.9 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(lwp .AND. lflush) CALL flush(numout)
100# if defined key_diaeiv 
101         IF( cdtype == 'TRA') THEN
102            u_eiv(:,:,:) = 0.e0
103            v_eiv(:,:,:) = 0.e0
104            w_eiv(:,:,:) = 0.e0
105         END IF
106# endif
107      ENDIF
108
109      zu_eiv(:,:) = 0.e0   ;   zv_eiv(:,:) = 0.e0   ;    zw_eiv(:,:) = 0.e0 
110     
111                                                    ! =================
112      DO jk = 1, jpkm1                              !  Horizontal slab
113         !                                          ! =================
114         DO jj = 1, jpjm1
115            DO ji = 1, fs_jpim1   ! vector opt.
116               zuwk = ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk  ) ) * fsaeiu(ji,jj,jk  ) * umask(ji,jj,jk  )
117               zuwk1= ( wslpi(ji,jj,jk+1) + wslpi(ji+1,jj,jk+1) ) * fsaeiu(ji,jj,jk+1) * umask(ji,jj,jk+1)
118               zvwk = ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk  ) ) * fsaeiv(ji,jj,jk  ) * vmask(ji,jj,jk  )
119               zvwk1= ( wslpj(ji,jj,jk+1) + wslpj(ji,jj+1,jk+1) ) * fsaeiv(ji,jj,jk+1) * vmask(ji,jj,jk+1)
120
121               zu_eiv(ji,jj) = 0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 ) 
122               zv_eiv(ji,jj) = 0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 ) 
123   
124               pun(ji,jj,jk) = pun(ji,jj,jk) + e2u(ji,jj) * zu_eiv(ji,jj)
125               pvn(ji,jj,jk) = pvn(ji,jj,jk) + e1v(ji,jj) * zv_eiv(ji,jj)
126            END DO
127         END DO
128# if defined key_diaeiv 
129         IF( cdtype == 'TRA') THEN
130            u_eiv(:,:,jk) = zu_eiv(:,:) / fse3u(:,:,jk)
131            v_eiv(:,:,jk) = zv_eiv(:,:) / fse3v(:,:,jk)
132         END IF
133# endif
134         IF( jk >=2 ) THEN                             ! jk=1 zw_eiv=0, not computed
135            DO jj = 2, jpjm1
136               DO ji = fs_2, fs_jpim1   ! vector opt.
137# if defined key_traldf_c2d || defined key_traldf_c3d
138                  zuwi  = ( wslpi(ji,jj,jk)+wslpi(ji-1,jj,jk) ) * fsaeiu(ji-1,jj,jk) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
139                  zuwi1 = ( wslpi(ji,jj,jk)+wslpi(ji+1,jj,jk) ) * fsaeiu(ji  ,jj,jk) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
140                  zvwj  = ( wslpj(ji,jj,jk)+wslpj(ji,jj-1,jk) ) * fsaeiv(ji,jj-1,jk) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
141                  zvwj1 = ( wslpj(ji,jj,jk)+wslpj(ji,jj+1,jk) ) * fsaeiv(ji,jj  ,jk) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
142 
143                  zw_eiv(ji,jj) = - 0.5 * tmask(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj ) 
144# else
145                  zuwi  = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) ) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
146                  zuwi1 = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) ) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
147                  zvwj  = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) ) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
148                  zvwj1 = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) ) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
149
150                  zw_eiv(ji,jj) = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj )
151# endif
152                  pwn(ji,jj,jk) = pwn(ji,jj,jk) + zw_eiv(ji,jj)
153               END DO
154            END DO
155# if defined key_diaeiv 
156            IF( cdtype == 'TRA')  w_eiv(:,:,jk) = zw_eiv(:,:) / ( e1t(:,:) * e2t(:,:) )
157# endif
158         ENDIF
159         !                                          ! =================
160      END DO                                        !    End of slab 
161      !                                             ! =================
162
163# if defined key_diaeiv 
164      IF( cdtype == 'TRA') THEN
165         CALL iom_put( "uoce_eiv", u_eiv )    ! i-eiv current
166         CALL iom_put( "voce_eiv", v_eiv )    ! j-eiv current
167         CALL iom_put( "woce_eiv", w_eiv )    ! vert. eiv current
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         IF( iom_use("ueiv_masstr") .OR. iom_use("ueiv_heattr") .OR. iom_use('ueiv_heattr3d')        &
176                                    .OR. iom_use("ueiv_salttr") .OR. iom_use('ueiv_salttr3d') ) THEN
177            z3d(:,:,jpk) = 0.e0
178            z2d(:,:) = 0.e0
179            DO jk = 1, jpkm1
180               z3d(:,:,jk) = rau0 * u_eiv(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)
181               z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
182            END DO
183            CALL iom_put( "ueiv_masstr", z3d )                  ! mass transport in i-direction
184         ENDIF
185
186         IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN
187            zztmp = 0.5 * rcp 
188            z2d(:,:) = 0.e0 
189            z3d_T(:,:,:) = 0.e0 
190            DO jk = 1, jpkm1
191               DO jj = 2, jpjm1
192                  DO ji = fs_2, fs_jpim1   ! vector opt.
193                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
194                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk) 
195                  END DO
196               END DO
197            END DO
198            IF (iom_use('ueiv_heattr') ) THEN
199               CALL lbc_lnk( z2d, 'U', -1. )
200               CALL iom_put( "ueiv_heattr", zztmp * z2d )                  ! 2D heat transport in i-direction
201            ENDIF
202            IF (iom_use('ueiv_heattr3d') ) THEN
203               CALL lbc_lnk( z3d_T, 'U', -1. )
204               CALL iom_put( "ueiv_heattr3d", zztmp * z3d_T )              ! 3D heat transport in i-direction
205            ENDIF
206         ENDIF
207
208         IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d') ) THEN
209            zztmp = 0.5 * 0.001
210            z2d(:,:) = 0.e0 
211            z3d_T(:,:,:) = 0.e0 
212            DO jk = 1, jpkm1
213               DO jj = 2, jpjm1
214                  DO ji = fs_2, fs_jpim1   ! vector opt.
215                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
216                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk) 
217                  END DO
218               END DO
219            END DO
220            IF (iom_use('ueiv_salttr') ) THEN
221               CALL lbc_lnk( z2d, 'U', -1. )
222               CALL iom_put( "ueiv_salttr", zztmp * z2d )                  ! 2D salt transport in i-direction
223            ENDIF
224            IF (iom_use('ueiv_salttr3d') ) THEN
225               CALL lbc_lnk( z3d_T, 'U', -1. )
226               CALL iom_put( "ueiv_salttr3d", zztmp * z3d_T )              ! 3D salt transport in i-direction
227            ENDIF
228         ENDIF
229
230         IF( iom_use("veiv_masstr") .OR. iom_use("veiv_heattr") .OR. iom_use('veiv_heattr3d')       &
231                                    .OR. iom_use("veiv_salttr") .OR. iom_use('veiv_salttr3d') ) THEN
232            z3d(:,:,jpk) = 0.e0
233            DO jk = 1, jpkm1
234               z3d(:,:,jk) = rau0 * v_eiv(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
235            END DO
236            CALL iom_put( "veiv_masstr", z3d )                  ! mass transport in j-direction
237         ENDIF
238           
239         IF( iom_use('veiv_heattr') .OR. iom_use('veiv_heattr3d') ) THEN
240            zztmp = 0.5 * rcp 
241            z2d(:,:) = 0.e0 
242            z3d_T(:,:,:) = 0.e0 
243            DO jk = 1, jpkm1
244               DO jj = 2, jpjm1
245                  DO ji = fs_2, fs_jpim1   ! vector opt.
246                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
247                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk) 
248                  END DO
249               END DO
250            END DO
251            IF (iom_use('veiv_heattr') ) THEN
252               CALL lbc_lnk( z2d, 'V', -1. )
253               CALL iom_put( "veiv_heattr", zztmp * z2d )                  ! 2D heat transport in j-direction
254            ENDIF
255            IF (iom_use('veiv_heattr3d') ) THEN
256               CALL lbc_lnk( z3d_T, 'V', -1. )
257               CALL iom_put( "veiv_heattr3d", zztmp * z3d_T )              ! 3D heat transport in j-direction
258            ENDIF
259         ENDIF
260
261         IF( iom_use('veiv_salttr') .OR. iom_use('veiv_salttr3d') ) THEN
262            zztmp = 0.5 * 0.001
263            z2d(:,:) = 0.e0 
264            z3d_T(:,:,:) = 0.e0 
265            DO jk = 1, jpkm1
266               DO jj = 2, jpjm1
267                  DO ji = fs_2, fs_jpim1   ! vector opt.
268                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
269                     z2d(ji,jj) = z2d(ji,jj) + z3d_T(ji,jj,jk)
270                  END DO
271               END DO
272            END DO
273            IF (iom_use('veiv_salttr') ) THEN
274               CALL lbc_lnk( z2d, 'V', -1. )
275               CALL iom_put( "veiv_salttr", zztmp * z2d )                  ! 2D salt transport in i-direction
276            ENDIF
277            IF (iom_use('veiv_salttr3d') ) THEN
278               CALL lbc_lnk( z3d_T, 'V', -1. )
279               CALL iom_put( "veiv_salttr3d", zztmp * z3d_T )              ! 3D salt transport in i-direction
280            ENDIF
281         ENDIF
282
283         IF( iom_use('weiv_masstr') .OR. iom_use('weiv_heattr3d') .OR. iom_use('weiv_salttr3d')) THEN   ! vertical mass transport & its square value
284           z2d(:,:) = rau0 * e12t(:,:)
285           DO jk = 1, jpk
286              z3d(:,:,jk) = w_eiv(:,:,jk) * z2d(:,:)
287           END DO
288           CALL iom_put( "weiv_masstr" , z3d )                  ! mass transport in k-direction
289         ENDIF
290
291         IF( iom_use('weiv_heattr3d') ) THEN
292            zztmp = 0.5 * rcp 
293            DO jk = 1, jpkm1
294               DO jj = 2, jpjm1
295                  DO ji = fs_2, fs_jpim1   ! vector opt.
296                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj,jk+1,jp_tem) )
297                  END DO
298               END DO
299            END DO
300            CALL lbc_lnk( z3d_T, 'T', 1. )
301            CALL iom_put( "weiv_heattr3d", zztmp * z3d_T )                 ! 3D heat transport in k-direction
302         ENDIF
303
304         IF( iom_use('weiv_salttr3d') ) THEN
305            zztmp = 0.5 * 0.001 
306            DO jk = 1, jpkm1
307               DO jj = 2, jpjm1
308                  DO ji = fs_2, fs_jpim1   ! vector opt.
309                     z3d_T(ji,jj,jk) = z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj,jk+1,jp_sal) )
310                  END DO
311               END DO
312            END DO
313            CALL lbc_lnk( z3d_T, 'T', 1. )
314            CALL iom_put( "weiv_salttr3d", zztmp * z3d_T )                 ! 3D salt transport in k-direction
315         ENDIF
316
317    END IF
318!
319    IF( ln_diaptr .AND. cdtype == 'TRA' ) THEN
320       z3d(:,:,:) = 0._wp
321       DO jk = 1, jpkm1
322          DO jj = 2, jpjm1
323             DO ji = fs_2, fs_jpim1   ! vector opt.
324                z3d(ji,jj,jk) = v_eiv(ji,jj,jk) * 0.5 * (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) &
325                &             * e1v(ji,jj) * fse3v(ji,jj,jk)
326             END DO
327          END DO
328       END DO
329       CALL dia_ptr_ohst_components( jp_tem, 'eiv', z3d )
330       z3d(:,:,:) = 0._wp
331       DO jk = 1, jpkm1
332          DO jj = 2, jpjm1
333             DO ji = fs_2, fs_jpim1   ! vector opt.
334                z3d(ji,jj,jk) = v_eiv(ji,jj,jk) * 0.5 * (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) &
335                &             * e1v(ji,jj) * fse3v(ji,jj,jk)
336             END DO
337          END DO
338       END DO
339       CALL dia_ptr_ohst_components( jp_sal, 'eiv', z3d )
340    ENDIF
341
342    IF( ln_KE_trd ) CALL trd_dyn(u_eiv, v_eiv, jpdyn_eivke, kt )
343# endif 
344
345# if defined key_diaeiv 
346      CALL wrk_dealloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv, z2d )
347      CALL wrk_dealloc( jpi, jpj, jpk, z3d, z3d_T )
348# else
349      CALL wrk_dealloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv )
350# endif
351      !
352      IF( nn_timing == 1 )  CALL timing_stop( 'tra_adv_eiv')
353      !
354    END SUBROUTINE tra_adv_eiv
355
356#else
357   !!----------------------------------------------------------------------
358   !!   Dummy module :             No rotation of the lateral mixing tensor
359   !!----------------------------------------------------------------------
360CONTAINS
361   SUBROUTINE tra_adv_eiv( kt, kit000, pun, pvn, pwn, cdtype )              ! Empty routine
362      INTEGER  ::   kt   
363      INTEGER  ::   kit000   
364      CHARACTER(len=3) ::   cdtype
365      REAL, DIMENSION(:,:,:) ::   pun, pvn, pwn
366      WRITE(*,*) 'tra_adv_eiv: You should not have seen this print! error?', &
367          &  kt, cdtype, pun(1,1,1), pvn(1,1,1), pwn(1,1,1)
368   END SUBROUTINE tra_adv_eiv
369#endif
370
371   !!==============================================================================
372END MODULE traadv_eiv
Note: See TracBrowser for help on using the repository browser.