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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90 @ 4431

Last change on this file since 4431 was 4409, checked in by trackstand2, 10 years ago

Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.
  • Property svn:keywords set to Id
File size: 11.2 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   USE diaar5, ONLY:   lk_diaar5
28# endif 
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   tra_adv_eiv   ! routine called by step.F90
34
35   !! * Control permutation of array indices
36#  include "oce_ftrans.h90"
37#  include "dom_oce_ftrans.h90"
38#  include "trc_oce_ftrans.h90"
39#  include "ldftra_oce_ftrans.h90"
40#  include "ldfslp_ftrans.h90"
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "ldftra_substitute.h90"
45#  include "ldfeiv_substitute.h90"
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
49   !! $Id$
50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE tra_adv_eiv( kt, pun, pvn, pwn, cdtype )
55      !!----------------------------------------------------------------------
56      !!                  ***  ROUTINE tra_adv_eiv  ***
57      !!
58      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
59      !!      trend and add it to the general trend of tracer equation.
60      !!
61      !! ** Method  :   The eddy induced advection is computed from the slope
62      !!      of iso-neutral surfaces computed in routine ldf_slp as follows:
63      !!         zu_eiv =  1/(e2u e3u)   dk[ aeiu e2u mi(wslpi) ]
64      !!         zv_eiv =  1/(e1v e3v)   dk[ aeiv e1v mj(wslpj)
65      !!         zw_eiv = -1/(e1t e2t) { di[ aeiu e2u mi(wslpi) ]
66      !!                               + dj[ aeiv e1v mj(wslpj) ] }
67      !!      add the eiv component to the model velocity:
68      !!         p.n = p.n + z._eiv
69      !!
70      !! ** Action  : - add to p.n the eiv component
71      !!----------------------------------------------------------------------
72      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
73      USE wrk_nemo, ONLY:   zu_eiv => wrk_2d_1 , zv_eiv => wrk_2d_2 , zw_eiv => wrk_2d_3   ! 2D workspace
74# if defined key_diaeiv 
75      USE wrk_nemo, ONLY:   z2d => wrk_2d_4   ! 2D workspace
76#endif
77      INTEGER                         , INTENT(in   ) ::   kt       ! ocean time-step index
78      CHARACTER(len=3)                , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
79
80      !! DCSE_NEMO: This style defeats ftrans
81!     REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun      ! in : 3 ocean velocity components
82!     REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pvn      ! out: 3 ocean velocity components
83!     REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pwn      ! increased by the eiv
84
85!FTRANS pun pvn pwn :I :I :z
86      REAL(wp), INTENT(inout) ::   pun(jpi,jpj,jpkorig)      ! in : 3 ocean velocity components
87      REAL(wp), INTENT(inout) ::   pvn(jpi,jpj,jpkorig)      ! out: 3 ocean velocity components
88      REAL(wp), INTENT(inout) ::   pwn(jpi,jpj,jpkorig)      ! increased by the eiv
89      !!
90      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
91      REAL(wp) ::   zuwk, zuwk1, zuwi, zuwi1   ! local scalars
92      REAL(wp) ::   zvwk, zvwk1, zvwj, zvwj1   !   -      -
93# if defined key_diaeiv 
94      REAL(wp) ::   zztmp                      ! local scalar
95# endif 
96      !!----------------------------------------------------------------------
97
98# if defined key_diaeiv 
99      IF( wrk_in_use(2, 1,2,3,4) ) THEN
100# else
101      IF( wrk_in_use(2, 1,2,3)   ) THEN
102# endif
103         CALL ctl_stop('tra_adv_eiv: requested workspace arrays are unavailable')   ;   RETURN
104      ENDIF
105
106      IF( kt == nit000 )  THEN
107         IF(lwp) WRITE(numout,*)
108         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection on ', cdtype,' :'
109         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component'
110# if defined key_diaeiv 
111         IF( cdtype == 'TRA') THEN
112            u_eiv(:,:,:) = 0.e0
113            v_eiv(:,:,:) = 0.e0
114            w_eiv(:,:,:) = 0.e0
115         END IF
116# endif
117      ENDIF
118
119      zu_eiv(:,:) = 0.e0   ;   zv_eiv(:,:) = 0.e0   ;    zw_eiv(:,:) = 0.e0 
120     
121!!DCSE_NEMO: TODO - restucture loop(s) so that loop over levels is innermost
122                                                    ! =================
123      DO jk = 1, jpkm1                              !  Horizontal slab
124         !                                          ! =================
125         DO jj = 1, jpjm1
126            DO ji = 1, fs_jpim1   ! vector opt.
127               zuwk = ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk  ) ) * fsaeiu(ji,jj,jk  ) * umask(ji,jj,jk  )
128               zuwk1= ( wslpi(ji,jj,jk+1) + wslpi(ji+1,jj,jk+1) ) * fsaeiu(ji,jj,jk+1) * umask(ji,jj,jk+1)
129               zvwk = ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk  ) ) * fsaeiv(ji,jj,jk  ) * vmask(ji,jj,jk  )
130               zvwk1= ( wslpj(ji,jj,jk+1) + wslpj(ji,jj+1,jk+1) ) * fsaeiv(ji,jj,jk+1) * vmask(ji,jj,jk+1)
131
132               zu_eiv(ji,jj) = 0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 ) 
133               zv_eiv(ji,jj) = 0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 ) 
134   
135               pun(ji,jj,jk) = pun(ji,jj,jk) + e2u(ji,jj) * zu_eiv(ji,jj)
136               pvn(ji,jj,jk) = pvn(ji,jj,jk) + e1v(ji,jj) * zv_eiv(ji,jj)
137            END DO
138         END DO
139# if defined key_diaeiv 
140         IF( cdtype == 'TRA') THEN
141            u_eiv(:,:,jk) = zu_eiv(:,:) / fse3u(:,:,jk)
142            v_eiv(:,:,jk) = zv_eiv(:,:) / fse3v(:,:,jk)
143         END IF
144# endif
145         IF( jk >=2 ) THEN                             ! jk=1 zw_eiv=0, not computed
146            DO jj = 2, jpjm1
147               DO ji = fs_2, fs_jpim1   ! vector opt.
148# if defined key_traldf_c2d || defined key_traldf_c3d
149                  zuwi  = ( wslpi(ji,jj,jk)+wslpi(ji-1,jj,jk) ) * fsaeiu(ji-1,jj,jk) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
150                  zuwi1 = ( wslpi(ji,jj,jk)+wslpi(ji+1,jj,jk) ) * fsaeiu(ji  ,jj,jk) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
151                  zvwj  = ( wslpj(ji,jj,jk)+wslpj(ji,jj-1,jk) ) * fsaeiv(ji,jj-1,jk) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
152                  zvwj1 = ( wslpj(ji,jj,jk)+wslpj(ji,jj+1,jk) ) * fsaeiv(ji,jj  ,jk) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
153 
154                  zw_eiv(ji,jj) = - 0.5 * tmask(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj ) 
155# else
156                  zuwi  = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) ) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
157                  zuwi1 = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) ) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
158                  zvwj  = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) ) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
159                  zvwj1 = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) ) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
160
161                  zw_eiv(ji,jj) = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj )
162# endif
163                  pwn(ji,jj,jk) = pwn(ji,jj,jk) + zw_eiv(ji,jj)
164               END DO
165            END DO
166# if defined key_diaeiv 
167            IF( cdtype == 'TRA')  w_eiv(:,:,jk) = zw_eiv(:,:) / ( e1t(:,:) * e2t(:,:) )
168# endif
169         ENDIF
170         !                                          ! =================
171      END DO                                        !    End of slab 
172      !                                             ! =================
173
174# if defined key_diaeiv 
175      IF( cdtype == 'TRA') THEN
176         CALL iom_put( "uoce_eiv", u_eiv )    ! i-eiv current
177         CALL iom_put( "voce_eiv", v_eiv )    ! j-eiv current
178         CALL iom_put( "woce_eiv", w_eiv )    ! vert. eiv current
179         IF( lk_diaar5 ) THEN
180            zztmp = 0.5 * rau0 * rcp 
181            z2d(:,:) = 0.e0 
182#if defined key_z_first
183            DO jj = 2, jpjm1
184               DO ji = fs_2, fs_jpim1   ! vector opt.
185                  DO jk = 1, jpkm1
186#else
187            DO jk = 1, jpkm1
188               DO jj = 2, jpjm1
189                  DO ji = fs_2, fs_jpim1   ! vector opt.
190#endif
191                     z2d(ji,jj) = z2d(ji,jj) + zztmp * u_eiv(ji,jj,jk) &
192                       &         * (tsn(ji,jj,jk,jp_tem)+tsn(ji+1,jj,jk,jp_tem)) * e1u(ji,jj) * fse3u(ji,jj,jk) 
193                  END DO
194               END DO
195            END DO
196            CALL lbc_lnk( z2d, 'U', -1. )
197            CALL iom_put( "ueiv_heattr", z2d )                  ! heat transport in i-direction
198            z2d(:,:) = 0.e0 
199#if defined key_z_first
200            DO jj = 2, jpjm1
201               DO ji = fs_2, fs_jpim1   ! vector opt.
202                  DO jk = 1, jpkm1
203#else
204            DO jk = 1, jpkm1
205               DO jj = 2, jpjm1
206                  DO ji = fs_2, fs_jpim1   ! vector opt.
207#endif
208                     z2d(ji,jj) = z2d(ji,jj) + zztmp * v_eiv(ji,jj,jk) &
209                     &           * (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * e2v(ji,jj) * fse3v(ji,jj,jk) 
210                  END DO
211               END DO
212            END DO
213            CALL lbc_lnk( z2d, 'V', -1. )
214            CALL iom_put( "veiv_heattr", z2d )                  !  heat transport in i-direction
215         ENDIF
216    END IF
217# endif 
218      !
219# if defined key_diaeiv 
220      IF( wrk_not_released(2, 1,2,3,4) )   CALL ctl_stop('tra_adv_eiv: failed to release workspace arrays')
221# else
222      IF( wrk_not_released(2, 1,2,3)   )   CALL ctl_stop('tra_adv_eiv: failed to release workspace arrays')
223# endif
224      !
225    END SUBROUTINE tra_adv_eiv
226
227#else
228   !!----------------------------------------------------------------------
229   !!   Dummy module :             No rotation of the lateral mixing tensor
230   !!----------------------------------------------------------------------
231CONTAINS
232   SUBROUTINE tra_adv_eiv( kt, pun, pvn, pwn, cdtype )              ! Empty routine
233      INTEGER  ::   kt   
234      CHARACTER(len=3) ::   cdtype
235      REAL, DIMENSION(:,:,:) ::   pun, pvn, pwn
236      WRITE(*,*) 'tra_adv_eiv: You should not have seen this print! error?', kt, cdtype, pun(1,1,1), pvn(1,1,1), pwn(1,1,1)
237   END SUBROUTINE tra_adv_eiv
238#endif
239
240   !!==============================================================================
241END MODULE traadv_eiv
Note: See TracBrowser for help on using the repository browser.