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

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90 @ 2524

Last change on this file since 2524 was 1756, checked in by smasson, 14 years ago

implement AR5 diagnostics, see ticket:610
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 9.0 KB
Line 
1MODULE traadv_eiv
2   !!======================================================================
3   !!                    ***  MODULE  traadv_eiv  ***
4   !! Ocean active tracers:  advection trend - eddy induced velocity
5   !!======================================================================
6   !! History :  9.0  !  05-11  (G. Madec)  Original code, from traldf and zdf _iso
7   !!----------------------------------------------------------------------
8#if defined key_traldf_eiv   ||   defined key_esopa
9   !!----------------------------------------------------------------------
10   !!   'key_traldf_eiv'                  rotation of the lateral mixing tensor
11   !!----------------------------------------------------------------------
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# if defined key_diaeiv
24   USE phycst          ! physical constants
25   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
26   USE diaar5, ONLY :   lk_diaar5
27# endif 
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   tra_adv_eiv   ! routine called by step.F90
33
34   !! * Substitutions
35#  include "domzgr_substitute.h90"
36#  include "ldftra_substitute.h90"
37#  include "ldfeiv_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !!  OPA 9.0 , LOCEAN-IPSL (2006)
41   !! $Id$
42   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44
45CONTAINS
46
47   SUBROUTINE tra_adv_eiv( kt, pun, pvn, pwn )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE tra_adv_eiv  ***
50      !!
51      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
52      !!      trend and add it to the general trend of tracer equation.
53      !!
54      !! ** Method  :   The eddy induced advection is computed from the slope
55      !!      of iso-neutral surfaces computed in routine ldf_slp as follows:
56      !!         zu_eiv =  1/(e2u e3u)   dk[ aeiu e2u mi(wslpi) ]
57      !!         zv_eiv =  1/(e1v e3v)   dk[ aeiv e1v mj(wslpj)
58      !!         zw_eiv = -1/(e1t e2t) { di[ aeiu e2u mi(wslpi) ]
59      !!                               + dj[ aeiv e1v mj(wslpj) ] }
60      !!      add the eiv component to the model velocity:
61      !!         p.n = p.n + z._eiv
62      !!
63      !! ** Action  : - add to p.n the eiv component
64      !!----------------------------------------------------------------------
65      INTEGER , INTENT(in   )                         ::   kt    ! ocean time-step index
66      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun   ! in : 3 ocean velocity components
67      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pvn   ! out: 3 ocean velocity components
68      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pwn   !      increased by the eiv
69      !!
70      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
71      REAL(wp) ::   zuwk, zuwk1, zuwi, zuwi1   ! temporary scalar
72      REAL(wp) ::   zvwk, zvwk1, zvwj, zvwj1   !    "         "
73      REAL(wp) ::   zu_eiv, zv_eiv, zw_eiv     !    "         "
74# if defined key_diaeiv
75      REAL(wp) ::   zztmp                      !    "         "
76      REAL(wp), DIMENSION(jpi,jpj) ::   z2d    !    "         "
77# endif 
78      !!----------------------------------------------------------------------
79
80      IF( kt == nit000 ) THEN
81         IF(lwp) WRITE(numout,*)
82         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection :'
83         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component'
84# if defined key_diaeiv
85         u_eiv(:,:,:) = 0.e0
86         v_eiv(:,:,:) = 0.e0
87         w_eiv(:,:,:) = 0.e0
88# endif
89      ENDIF
90      !                                             ! =================
91      DO jk = 1, jpkm1                              !  Horizontal slab
92         !                                          ! =================
93         DO jj = 1, jpjm1
94            DO ji = 1, fs_jpim1   ! vector opt.
95               zuwk = ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk  ) ) * fsaeiu(ji,jj,jk  ) * umask(ji,jj,jk  )
96               zuwk1= ( wslpi(ji,jj,jk+1) + wslpi(ji+1,jj,jk+1) ) * fsaeiu(ji,jj,jk+1) * umask(ji,jj,jk+1)
97               zvwk = ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk  ) ) * fsaeiv(ji,jj,jk  ) * vmask(ji,jj,jk  )
98               zvwk1= ( wslpj(ji,jj,jk+1) + wslpj(ji,jj+1,jk+1) ) * fsaeiv(ji,jj,jk+1) * vmask(ji,jj,jk+1)
99
100               zu_eiv = 0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 ) / fse3u(ji,jj,jk)
101               zv_eiv = 0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 ) / fse3v(ji,jj,jk)
102   
103               pun(ji,jj,jk) = pun(ji,jj,jk) + zu_eiv
104               pvn(ji,jj,jk) = pvn(ji,jj,jk) + zv_eiv
105# if defined key_diaeiv
106               u_eiv(ji,jj,jk) = zu_eiv
107               v_eiv(ji,jj,jk) = zv_eiv
108# endif
109            END DO
110         END DO
111         IF( jk >=2 ) THEN                             ! jk=1 zw_eiv=0, not computed
112            DO jj = 2, jpjm1
113               DO ji = fs_2, fs_jpim1   ! vector opt.
114# if defined key_traldf_c2d || defined key_traldf_c3d
115                  zuwi  = ( wslpi(ji,jj,jk)+wslpi(ji-1,jj,jk) ) * fsaeiu(ji-1,jj,jk) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
116                  zuwi1 = ( wslpi(ji,jj,jk)+wslpi(ji+1,jj,jk) ) * fsaeiu(ji  ,jj,jk) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
117                  zvwj  = ( wslpj(ji,jj,jk)+wslpj(ji,jj-1,jk) ) * fsaeiv(ji,jj-1,jk) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
118                  zvwj1 = ( wslpj(ji,jj,jk)+wslpj(ji,jj+1,jk) ) * fsaeiv(ji,jj  ,jk) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
119 
120                  zw_eiv = - 0.5 * tmask(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj ) / ( e1t(ji,jj)*e2t(ji,jj) )
121# else
122                  zuwi  = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) ) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
123                  zuwi1 = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) ) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
124                  zvwj  = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) ) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
125                  zvwj1 = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) ) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
126
127                  zw_eiv = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj )  &
128                     &                                                / ( e1t(ji,jj)*e2t(ji,jj) )
129# endif
130                  pwn(ji,jj,jk) = pwn(ji,jj,jk) + zw_eiv
131
132# if defined key_diaeiv
133                  w_eiv(ji,jj,jk) = zw_eiv
134# endif
135               END DO
136            END DO
137         ENDIF
138         !                                          ! =================
139      END DO                                        !    End of slab 
140      !                                             ! =================
141
142# if defined key_diaeiv
143      CALL iom_put( "uoce_eiv", u_eiv )    ! i-eiv current
144      CALL iom_put( "voce_eiv", v_eiv )    ! j-eiv current
145      CALL iom_put( "woce_eiv", w_eiv )    ! vert. eiv current
146      IF( lk_diaar5 ) THEN
147         zztmp = 0.5 * rau0 * rcp 
148         z2d(:,:) = 0.e0 
149         DO jk = 1, jpkm1
150            DO jj = 2, jpjm1
151               DO ji = fs_2, fs_jpim1   ! vector opt.
152                  z2d(ji,jj) = z2d(ji,jj) + zztmp * u_eiv(ji,jj,jk) * (tn(ji,jj,jk)+tn(ji+1,jj,jk)) * e1u(ji,jj) * fse3u(ji,jj,jk) 
153               END DO
154            END DO
155         END DO
156         CALL lbc_lnk( z2d, 'U', -1. )
157         CALL iom_put( "ueiv_heattr", z2d )                  ! heat transport in i-direction
158         z2d(:,:) = 0.e0 
159         DO jk = 1, jpkm1
160            DO jj = 2, jpjm1
161               DO ji = fs_2, fs_jpim1   ! vector opt.
162                  z2d(ji,jj) = z2d(ji,jj) + zztmp * v_eiv(ji,jj,jk) * (tn(ji,jj,jk)+tn(ji,jj+1,jk)) * e2v(ji,jj) * fse3v(ji,jj,jk) 
163               END DO
164            END DO
165         END DO
166         CALL lbc_lnk( z2d, 'V', -1. )
167         CALL iom_put( "veiv_heattr", z2d )                  !  heat transport in i-direction
168      ENDIF
169# endif 
170      !
171    END SUBROUTINE tra_adv_eiv
172
173#else
174   !!----------------------------------------------------------------------
175   !!   Dummy module :             No rotation of the lateral mixing tensor
176   !!----------------------------------------------------------------------
177CONTAINS
178   SUBROUTINE tra_adv_eiv( kt, pun, pvn, pwn )              ! Empty routine
179      REAL, DIMENSION(:,:,:) ::   pun, pvn, pwn
180      WRITE(*,*) 'tra_adv_eiv: You should not have seen this print! error?', kt, pun(1,1,1), pvn(1,1,1), pwn(1,1,1)
181   END SUBROUTINE tra_adv_eiv
182#endif
183
184   !!==============================================================================
185END MODULE traadv_eiv
Note: See TracBrowser for help on using the repository browser.