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_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90 @ 12555

Last change on this file since 12555 was 12555, checked in by charris, 4 years ago

Changes from GO6 package branch (GMED ticket 450):

svn merge -r 11035:11101 ​svn+ssh://charris@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/dev_r5518_GO6_package
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.