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

source: branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_eiv.F90 @ 791

Last change on this file since 791 was 791, checked in by gm, 16 years ago

dev_001_GM - TRA/traadv : switch from velocity to transport + optimised traadv_eiv2 introduced - compilation OK

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.5 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   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport
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
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   tra_adv_eiv   ! routine called by step.F90
28
29   !! * Substitutions
30#  include "domzgr_substitute.h90"
31#  include "ldftra_substitute.h90"
32#  include "ldfeiv_substitute.h90"
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)
36   !! $Id:$
37   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39
40CONTAINS
41
42   SUBROUTINE tra_adv_eiv( kt, pun, pvn, pwn )
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE tra_adv_eiv  ***
45      !!
46      !! ** Purpose :   Compute the eddy induced transport and add it to the
47      !!              effective transport
48      !!
49      !! ** Method  :   The eddy induced transport is computed from the slope
50      !!      of iso-neutral surfaces computed in routine ldf_slp as follows:
51      !!         zu_eiv =   dk[ aeiu e2u mi(wslpi) ]
52      !!         zv_eiv =   dk[ aeiv e1v mj(wslpj) ]
53      !!         zw_eiv = - { di[ aeiu e2u mi(wslpi) ] + dj[ aeiv e1v mj(wslpj) ] }
54      !!      add the eiv component to the model velocity:
55      !!         p.n = p.n + z._eiv
56      !! CAUTION : the horizontal transports not updated along jpi column and jpj row
57      !!           the vertical   transports not updated along 1 & jpi columns and 1 & jpj rows
58      !!
59      !! ** Action  : - add to p.n the eiv component
60      !!----------------------------------------------------------------------
61      INTEGER , INTENT(in   )                         ::   kt    ! ocean time-step index
62      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun   ! in : 3 ocean transport components
63      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pvn   ! out: 3 ocean transport components
64      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pwn   !      increased by the eiv
65      !!
66      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
67      REAL(wp) ::   zuwk, zuwk1, zuwi, zuwi1   ! temporary scalar
68      REAL(wp) ::   zvwk, zvwk1, zvwj, zvwj1   !    "         "
69      REAL(wp) ::   zu_eiv, zv_eiv, zw_eiv     !    "         "
70      !!----------------------------------------------------------------------
71
72      IF( kt == nit000 ) THEN
73         IF(lwp) WRITE(numout,*)
74         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection :'
75         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component'
76# if defined key_diaeiv
77         u_eiv(:,:,:) = 0.e0
78         v_eiv(:,:,:) = 0.e0
79         w_eiv(:,:,:) = 0.e0
80# endif
81      ENDIF
82      !                                             ! =================
83      DO jk = 1, jpkm1                              !  Horizontal slab
84         !                                          ! =================
85         DO jj = 1, jpjm1
86            DO ji = 1, fs_jpim1   ! vector opt.
87               zuwk = ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk  ) ) * fsaeiu(ji,jj,jk  ) * umask(ji,jj,jk  )
88               zuwk1= ( wslpi(ji,jj,jk+1) + wslpi(ji+1,jj,jk+1) ) * fsaeiu(ji,jj,jk+1) * umask(ji,jj,jk+1)
89               zvwk = ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk  ) ) * fsaeiv(ji,jj,jk  ) * vmask(ji,jj,jk  )
90               zvwk1= ( wslpj(ji,jj,jk+1) + wslpj(ji,jj+1,jk+1) ) * fsaeiv(ji,jj,jk+1) * vmask(ji,jj,jk+1)
91
92               zu_eiv = 0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 )
93               zv_eiv = 0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 )
94   
95               pun(ji,jj,jk) = pun(ji,jj,jk) + e2u(ji,jj) * zu_eiv
96               pvn(ji,jj,jk) = pvn(ji,jj,jk) + e1v(ji,jj) * zv_eiv
97# if defined key_diaeiv
98               u_eiv(ji,jj,jk) = zu_eiv / fse3u(ji,jj,jk)
99               v_eiv(ji,jj,jk) = zv_eiv / fse3v(ji,jj,jk)
100# endif
101            END DO
102         END DO
103         IF( jk >=2 ) THEN                             ! jk=1 zw_eiv=0, not computed
104            DO jj = 2, jpjm1
105               DO ji = fs_2, fs_jpim1   ! vector opt.
106# if defined key_traldf_c2d || defined key_traldf_c3d
107                  zuwi  = ( wslpi(ji,jj,jk)+wslpi(ji-1,jj,jk) ) * fsaeiu(ji-1,jj,jk) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
108                  zuwi1 = ( wslpi(ji,jj,jk)+wslpi(ji+1,jj,jk) ) * fsaeiu(ji  ,jj,jk) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
109                  zvwj  = ( wslpj(ji,jj,jk)+wslpj(ji,jj-1,jk) ) * fsaeiv(ji,jj-1,jk) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
110                  zvwj1 = ( wslpj(ji,jj,jk)+wslpj(ji,jj+1,jk) ) * fsaeiv(ji,jj  ,jk) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
111 
112                  zw_eiv = - 0.5 * tmask(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj ) / ( e1t(ji,jj)*e2t(ji,jj) )
113# else
114                  zuwi  = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) ) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
115                  zuwi1 = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) ) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
116                  zvwj  = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) ) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
117                  zvwj1 = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) ) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
118
119                  zw_eiv = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj )
120# endif
121                  pwn(ji,jj,jk) = pwn(ji,jj,jk) + zw_eiv
122
123# if defined key_diaeiv
124                  w_eiv(ji,jj,jk) = zw_eiv / ( e1t(ji,jj)*e2t(ji,jj) )
125# endif
126               END DO
127            END DO
128         ENDIF
129         !                                          ! =================
130      END DO                                        !    End of slab 
131      !                                             ! =================
132   END SUBROUTINE tra_adv_eiv
133
134
135!!gm   test tra_adv_eiv better (faster) coded?  to be verified
136
137   SUBROUTINE tra_adv_eiv2( kt, pun, pvn, pwn )
138      !!----------------------------------------------------------------------
139      !!                  ***  ROUTINE tra_adv_eiv  ***
140      !!
141      !! ** Purpose :   Compute the eddy induced transport and add it to the
142      !!              effective transport
143      !!
144      !! ** Method  :   The eddy induced transport is computed from the slope
145      !!              of iso-neutral surfaces (see ldfslp.F90) as follows:
146      !!                   zu_eiv =   dk[ aeiu e2u mi(wslpi) ]
147      !!                   zv_eiv =   dk[ aeiv e1v mj(wslpj) ]
148      !!                   zw_eiv = - { di[ aeiu e2u mi(wslpi) ] + dj[ aeiv e1v mj(wslpj) ] }
149      !!              add the eiv component to the model velocity:
150      !!                   p.n = p.n + z._eiv
151      !! CAUTION : the horizontal transports not updated along jpi column and jpj row
152      !!           the vertical   transports not updated along 1 & jpi columns and 1 & jpj rows
153      !!
154      !! ** Action  : - add to p.n the eiv transport component
155      !!----------------------------------------------------------------------
156      INTEGER , INTENT(in   )                         ::   kt    ! ocean time-step index
157      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun   ! in : 3 ocean transport components
158      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pvn   ! out: 3 ocean transport components
159      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pwn   !      increased by the eiv
160      !!
161      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
162      REAL(wp) ::   zuwslpi, zvwslpj           ! temporary scalar
163      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsfu, zsfv   ! eiv stream-function in u and v directions
164      !!----------------------------------------------------------------------
165
166      IF( kt == nit000 ) THEN
167         IF(lwp) WRITE(numout,*)
168         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection :'
169         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component'
170      ENDIF
171
172      ! eiv stream-function in u- and v-directions
173      ! NB: UW-point mask at level k is umask(:,:,k) idem form VW-point mask
174      zsfu(:,:, 1 ) = 0.e0   ;   zsfv(:,:, 1 ) = 0.e0      ! surface value set to zero
175      zsfu(:,:,jpk) = 0.e0   ;   zsfv(:,:,jpk) = 0.e0      ! bottom  value set to zero
176      DO jk = 2, jpkm1
177         DO jj = 1, jpjm1
178            DO ji = 1, fs_jpim1   ! vector opt.
179               zuwslpi = 0.5 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
180               zvwslpj = 0.5 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) ) 
181               zsfu(ji,jj,jk) = zuwslpi * e2u(ji,jj) * fsaeiu(ji,jj,jk) * umask(ji,jj,jk)
182               zsfv(ji,jj,jk) = zvwslpj * e1v(ji,jj) * fsaeiv(ji,jj,jk) * vmask(ji,jj,jk)
183            END DO
184         END DO
185      END DO
186# if defined key_diaeiv
187      ! save eiv stream function in the output
188!!gm  to be done, u_sfeiv and v_sfeiv not defined    ==> new IOM....
189!!gm  and zsfu, zfv notdefined for jpi column and jpj row
190!     u_sfeiv(:,:,:) = zsfu(:,:,:)
191!     v_sfeiv(:,:,:) = zsfu(:,:,:)
192# endif
193
194      ! increase the transport with the eiv transport
195      DO jk = 1, jpkm1
196         DO jj = 1, jpjm1
197            DO ji = 1, fs_jpim1   ! vector opt.
198               pun(ji,jj,jk) = pun(ji,jj,jk) + ( zsfu(ji,jj,jk) - zsfu(ji,jj,jk+1) )
199               pvn(ji,jj,jk) = pvn(ji,jj,jk) + ( zsfv(ji,jj,jk) - zsfv(ji,jj,jk+1) )
200            END DO
201         END DO
202      END DO
203      DO jk = 2, jpkm1
204         DO jj = 2, jpjm1
205            DO ji = fs_2, fs_jpim1   ! vector opt.
206               pwn(ji,jj,jk) = pwn(ji,jj,jk) - (  zsfu(ji,jj,jk) - zsfu(ji-1,jj  ,jk)                        &
207                  &                             + zsfv(ji,jj,jk) - zsfv(ji  ,jj-1,jk)  ) * tmask(ji,jj,jk)
208            END DO
209         END DO
210      END DO
211      !
212   END SUBROUTINE tra_adv_eiv2
213
214#else
215   !!----------------------------------------------------------------------
216   !!   Dummy module :             No rotation of the lateral mixing tensor
217   !!----------------------------------------------------------------------
218CONTAINS
219   SUBROUTINE tra_adv_eiv( kt, pun, pvn, pwn )              ! Empty routine
220      REAL, DIMENSION(:,:,:) ::   pun, pvn, pwn
221      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)
222   END SUBROUTINE tra_adv_eiv
223#endif
224
225   !!==============================================================================
226END MODULE traadv_eiv
Note: See TracBrowser for help on using the repository browser.