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

source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 6 years ago

Remove svn keywords
File size: 10.4 KB
RevLine 
[458]1MODULE traadv_eiv
[503]2   !!======================================================================
[458]3   !!                    ***  MODULE  traadv_eiv  ***
[2528]4   !! Ocean tracers:  advection trend - eddy induced velocity
[503]5   !!======================================================================
[2528]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
[458]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
[1482]22   USE iom
[2528]23   USE trc_oce         ! share passive tracers/Ocean variables
[1756]24# if defined key_diaeiv
25   USE phycst          ! physical constants
26   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
27# endif 
[3294]28   USE wrk_nemo        ! Memory Allocation
29   USE timing          ! Timing
[458]30
31   IMPLICIT NONE
32   PRIVATE
33
[503]34   PUBLIC   tra_adv_eiv   ! routine called by step.F90
[458]35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "ldftra_substitute.h90"
39#  include "ldfeiv_substitute.h90"
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
[2528]42   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
43   !! $Id$
[2715]44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[458]45   !!----------------------------------------------------------------------
46CONTAINS
47
[3294]48   SUBROUTINE tra_adv_eiv( kt, kit000, pun, pvn, pwn, cdtype )
[458]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
[503]65      !!----------------------------------------------------------------------
[2528]66      INTEGER                         , INTENT(in   ) ::   kt       ! ocean time-step index
[3294]67      INTEGER                         , INTENT(in   ) ::   kit000   ! first time step index
[2528]68      CHARACTER(len=3)                , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
69      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun      ! in : 3 ocean velocity components
70      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pvn      ! out: 3 ocean velocity components
71      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pwn      ! increased by the eiv
[458]72      !!
[503]73      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
[2528]74      REAL(wp) ::   zuwk, zuwk1, zuwi, zuwi1   ! local scalars
75      REAL(wp) ::   zvwk, zvwk1, zvwj, zvwj1   !   -      -
76# if defined key_diaeiv 
77      REAL(wp) ::   zztmp                      ! local scalar
[1756]78# endif 
[3294]79      REAL(wp), POINTER, DIMENSION(:,:) :: zu_eiv, zv_eiv, zw_eiv, z2d
[458]80      !!----------------------------------------------------------------------
[3294]81      !
82      IF( nn_timing == 1 )  CALL timing_start( 'tra_adv_eiv')
83      !
[2715]84# if defined key_diaeiv 
[3294]85      CALL wrk_alloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv, z2d )
[2715]86# else
[3294]87      CALL wrk_alloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv )
[2715]88# endif
89
[3294]90      IF( kt == kit000 )  THEN
[458]91         IF(lwp) WRITE(numout,*)
[2528]92         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection on ', cdtype,' :'
[458]93         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component'
[2528]94# if defined key_diaeiv 
95         IF( cdtype == 'TRA') THEN
96            u_eiv(:,:,:) = 0.e0
97            v_eiv(:,:,:) = 0.e0
98            w_eiv(:,:,:) = 0.e0
99         END IF
[458]100# endif
101      ENDIF
[2528]102
103      zu_eiv(:,:) = 0.e0   ;   zv_eiv(:,:) = 0.e0   ;    zw_eiv(:,:) = 0.e0 
104     
105                                                    ! =================
[458]106      DO jk = 1, jpkm1                              !  Horizontal slab
107         !                                          ! =================
108         DO jj = 1, jpjm1
109            DO ji = 1, fs_jpim1   ! vector opt.
110               zuwk = ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk  ) ) * fsaeiu(ji,jj,jk  ) * umask(ji,jj,jk  )
111               zuwk1= ( wslpi(ji,jj,jk+1) + wslpi(ji+1,jj,jk+1) ) * fsaeiu(ji,jj,jk+1) * umask(ji,jj,jk+1)
112               zvwk = ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk  ) ) * fsaeiv(ji,jj,jk  ) * vmask(ji,jj,jk  )
113               zvwk1= ( wslpj(ji,jj,jk+1) + wslpj(ji,jj+1,jk+1) ) * fsaeiv(ji,jj,jk+1) * vmask(ji,jj,jk+1)
114
[2528]115               zu_eiv(ji,jj) = 0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 ) 
116               zv_eiv(ji,jj) = 0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 ) 
[458]117   
[2528]118               pun(ji,jj,jk) = pun(ji,jj,jk) + e2u(ji,jj) * zu_eiv(ji,jj)
119               pvn(ji,jj,jk) = pvn(ji,jj,jk) + e1v(ji,jj) * zv_eiv(ji,jj)
[458]120            END DO
121         END DO
[2528]122# if defined key_diaeiv 
123         IF( cdtype == 'TRA') THEN
124            u_eiv(:,:,jk) = zu_eiv(:,:) / fse3u(:,:,jk)
125            v_eiv(:,:,jk) = zv_eiv(:,:) / fse3v(:,:,jk)
126         END IF
127# endif
[458]128         IF( jk >=2 ) THEN                             ! jk=1 zw_eiv=0, not computed
129            DO jj = 2, jpjm1
130               DO ji = fs_2, fs_jpim1   ! vector opt.
131# if defined key_traldf_c2d || defined key_traldf_c3d
132                  zuwi  = ( wslpi(ji,jj,jk)+wslpi(ji-1,jj,jk) ) * fsaeiu(ji-1,jj,jk) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
133                  zuwi1 = ( wslpi(ji,jj,jk)+wslpi(ji+1,jj,jk) ) * fsaeiu(ji  ,jj,jk) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
134                  zvwj  = ( wslpj(ji,jj,jk)+wslpj(ji,jj-1,jk) ) * fsaeiv(ji,jj-1,jk) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
135                  zvwj1 = ( wslpj(ji,jj,jk)+wslpj(ji,jj+1,jk) ) * fsaeiv(ji,jj  ,jk) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
136 
[2528]137                  zw_eiv(ji,jj) = - 0.5 * tmask(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj ) 
[458]138# else
[503]139                  zuwi  = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) ) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
140                  zuwi1 = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) ) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
141                  zvwj  = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) ) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
142                  zvwj1 = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) ) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
[458]143
[2528]144                  zw_eiv(ji,jj) = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj )
[458]145# endif
[2528]146                  pwn(ji,jj,jk) = pwn(ji,jj,jk) + zw_eiv(ji,jj)
[458]147               END DO
148            END DO
[2528]149# if defined key_diaeiv 
150            IF( cdtype == 'TRA')  w_eiv(:,:,jk) = zw_eiv(:,:) / ( e1t(:,:) * e2t(:,:) )
151# endif
[458]152         ENDIF
153         !                                          ! =================
154      END DO                                        !    End of slab 
155      !                                             ! =================
156
[2528]157# if defined key_diaeiv 
158      IF( cdtype == 'TRA') THEN
159         CALL iom_put( "uoce_eiv", u_eiv )    ! i-eiv current
160         CALL iom_put( "voce_eiv", v_eiv )    ! j-eiv current
161         CALL iom_put( "woce_eiv", w_eiv )    ! vert. eiv current
[4990]162         IF( iom_use('ueiv_heattr') ) THEN
[2528]163            zztmp = 0.5 * rau0 * rcp 
164            z2d(:,:) = 0.e0 
165            DO jk = 1, jpkm1
166               DO jj = 2, jpjm1
167                  DO ji = fs_2, fs_jpim1   ! vector opt.
[4990]168                     z2d(ji,jj) = z2d(ji,jj) + u_eiv(ji,jj,jk) &
[3787]169                       &         * (tsn(ji,jj,jk,jp_tem)+tsn(ji+1,jj,jk,jp_tem)) * e2u(ji,jj) * fse3u(ji,jj,jk) 
[2528]170                  END DO
[1756]171               END DO
172            END DO
[2528]173            CALL lbc_lnk( z2d, 'U', -1. )
[4990]174            CALL iom_put( "ueiv_heattr", zztmp * z2d )                  ! heat transport in i-direction
175         ENDIF
176           
177         IF( iom_use('veiv_heattr') ) THEN
178            zztmp = 0.5 * rau0 * rcp 
[2528]179            z2d(:,:) = 0.e0 
180            DO jk = 1, jpkm1
181               DO jj = 2, jpjm1
182                  DO ji = fs_2, fs_jpim1   ! vector opt.
[4990]183                     z2d(ji,jj) = z2d(ji,jj) + v_eiv(ji,jj,jk) &
[3787]184                     &           * (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * e1v(ji,jj) * fse3v(ji,jj,jk) 
[2528]185                  END DO
[1756]186               END DO
187            END DO
[2528]188            CALL lbc_lnk( z2d, 'V', -1. )
[4990]189            CALL iom_put( "veiv_heattr", zztmp * z2d )                  !  heat transport in i-direction
[2528]190         ENDIF
191    END IF
[1482]192# endif 
193      !
[2715]194# if defined key_diaeiv 
[3294]195      CALL wrk_dealloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv, z2d )
[2715]196# else
[3294]197      CALL wrk_dealloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv )
[2715]198# endif
199      !
[3294]200      IF( nn_timing == 1 )  CALL timing_stop( 'tra_adv_eiv')
201      !
[1482]202    END SUBROUTINE tra_adv_eiv
203
[458]204#else
205   !!----------------------------------------------------------------------
206   !!   Dummy module :             No rotation of the lateral mixing tensor
207   !!----------------------------------------------------------------------
208CONTAINS
[3294]209   SUBROUTINE tra_adv_eiv( kt, kit000, pun, pvn, pwn, cdtype )              ! Empty routine
[2528]210      INTEGER  ::   kt   
[3294]211      INTEGER  ::   kit000   
[2528]212      CHARACTER(len=3) ::   cdtype
[458]213      REAL, DIMENSION(:,:,:) ::   pun, pvn, pwn
[6204]214      WRITE(*,*) 'tra_adv_eiv: You should not have seen this print! error?', &
215          &  kt, cdtype, pun(1,1,1), pvn(1,1,1), pwn(1,1,1)
[458]216   END SUBROUTINE tra_adv_eiv
217#endif
218
219   !!==============================================================================
220END MODULE traadv_eiv
Note: See TracBrowser for help on using the repository browser.