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 @ 4401

Last change on this file since 4401 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • 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,jpk)      ! in : 3 ocean velocity components
87      REAL(wp), INTENT(inout) ::   pvn(jpi,jpj,jpk)      ! out: 3 ocean velocity components
88      REAL(wp), INTENT(inout) ::   pwn(jpi,jpj,jpk)      ! 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.