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

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_eiv.F90 @ 2024

Last change on this file since 2024 was 2024, checked in by cetlod, 14 years ago

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