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.
ldftra.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90 @ 9169

Last change on this file since 9169 was 9169, checked in by gm, 6 years ago

dev_merge_2017: all SRC: finalize the removal of useless warning when reading namelist_cfg + remove all nn_closea + nn_msh replaced by a logical
  • Property svn:keywords set to Id
File size: 42.4 KB
Line 
1MODULE ldftra
2   !!======================================================================
3   !!                       ***  MODULE  ldftra  ***
4   !! Ocean physics:  lateral diffusivity coefficients
5   !!=====================================================================
6   !! History :       ! 1997-07  (G. Madec)  from inimix.F split in 2 routines
7   !!   NEMO     1.0  ! 2002-09  (G. Madec)  F90: Free form and module
8   !!            2.0  ! 2005-11  (G. Madec) 
9   !!            3.7  ! 2013-12  (F. Lemarie, G. Madec)  restructuration/simplification of aht/aeiv specification,
10   !!                 !                                  add velocity dependent coefficient and optional read in file
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   ldf_tra_init : initialization, namelist read, and parameters control
15   !!   ldf_tra      : update lateral eddy diffusivity coefficients at each time step
16   !!   ldf_eiv_init : initialization of the eiv coeff. from namelist choices
17   !!   ldf_eiv      : time evolution of the eiv coefficients (function of the growth rate of baroclinic instability)
18   !!   ldf_eiv_trp  : add to the input ocean transport the contribution of the EIV parametrization
19   !!   ldf_eiv_dia  : diagnose the eddy induced velocity from the eiv streamfunction
20   !!----------------------------------------------------------------------
21   USE oce             ! ocean dynamics and tracers
22   USE dom_oce         ! ocean space and time domain
23   USE phycst          ! physical constants
24   USE ldfslp          ! lateral diffusion: slope of iso-neutral surfaces
25   USE ldfc1d_c2d      ! lateral diffusion: 1D & 2D cases
26   USE diaptr
27   !
28   USE in_out_manager  ! I/O manager
29   USE iom             ! I/O module for ehanced bottom friction file
30   USE lib_mpp         ! distribued memory computing library
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   ldf_tra_init   ! called by nemogcm.F90
37   PUBLIC   ldf_tra        ! called by step.F90
38   PUBLIC   ldf_eiv_init   ! called by nemogcm.F90
39   PUBLIC   ldf_eiv        ! called by step.F90
40   PUBLIC   ldf_eiv_trp    ! called by traadv.F90
41   PUBLIC   ldf_eiv_dia    ! called by traldf_iso and traldf_iso_triad.F90
42   
43   !                                   !!* Namelist namtra_ldf : lateral mixing on tracers *
44   !                                    != Operator type =!
45   LOGICAL , PUBLIC ::   ln_traldf_NONE      !: no operator: No explicit diffusion
46   LOGICAL , PUBLIC ::   ln_traldf_lap       !: laplacian operator
47   LOGICAL , PUBLIC ::   ln_traldf_blp       !: bilaplacian operator
48   !                                    != Direction of action =!
49   LOGICAL , PUBLIC ::   ln_traldf_lev       !: iso-level direction
50   LOGICAL , PUBLIC ::   ln_traldf_hor       !: horizontal (geopotential) direction
51!  LOGICAL , PUBLIC ::   ln_traldf_iso       !: iso-neutral direction                    (see ldfslp)
52!  LOGICAL , PUBLIC ::   ln_traldf_triad     !: griffies triad scheme                    (see ldfslp)
53   LOGICAL , PUBLIC ::   ln_traldf_msc       !: Method of Stabilizing Correction
54!  LOGICAL , PUBLIC ::   ln_triad_iso        !: pure horizontal mixing in ML             (see ldfslp)
55!  LOGICAL , PUBLIC ::   ln_botmix_triad     !: mixing on bottom                         (see ldfslp)
56!  REAL(wp), PUBLIC ::   rn_sw_triad         !: =1/0 switching triad / all 4 triads used (see ldfslp)
57!  REAL(wp), PUBLIC ::   rn_slpmax           !: slope limit                              (see ldfslp)
58   !                                    !=  Coefficients =!
59   INTEGER , PUBLIC ::   nn_aht_ijk_t        !: choice of time & space variations of the lateral eddy diffusivity coef.
60   REAL(wp), PUBLIC ::   rn_aht_0            !:   laplacian lateral eddy diffusivity [m2/s]
61   REAL(wp), PUBLIC ::   rn_bht_0            !: bilaplacian lateral eddy diffusivity [m4/s]
62
63   !                                   !!* Namelist namtra_ldfeiv : eddy induced velocity param. *
64   !                                    != Use/diagnose eiv =!
65   LOGICAL , PUBLIC ::   ln_ldfeiv           !: eddy induced velocity flag
66   LOGICAL , PUBLIC ::   ln_ldfeiv_dia       !: diagnose & output eiv streamfunction and velocity (IOM)
67   !                                    != Coefficients =!
68   INTEGER , PUBLIC ::   nn_aei_ijk_t        !: choice of time/space variation of the eiv coeff.
69   REAL(wp), PUBLIC ::   rn_aeiv_0           !: eddy induced velocity coefficient [m2/s]
70   
71   LOGICAL , PUBLIC ::   l_ldftra_time = .FALSE.   !: flag for time variation of the lateral eddy diffusivity coef.
72   LOGICAL , PUBLIC ::   l_ldfeiv_time = .FALSE.   ! flag for time variation of the eiv coef.
73
74   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahtu, ahtv   !: eddy diffusivity coef. at U- and V-points   [m2/s]
75   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aeiu, aeiv   !: eddy induced velocity coeff.                [m2/s]
76
77   REAL(wp) ::   r1_4  = 0.25_wp          ! =1/4
78   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! =1/12
79
80   !! * Substitutions
81#  include "vectopt_loop_substitute.h90"
82   !!----------------------------------------------------------------------
83   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
84   !! $Id$
85   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
86   !!----------------------------------------------------------------------
87CONTAINS
88
89   SUBROUTINE ldf_tra_init
90      !!----------------------------------------------------------------------
91      !!                  ***  ROUTINE ldf_tra_init  ***
92      !!
93      !! ** Purpose :   initializations of the tracer lateral mixing coeff.
94      !!
95      !! ** Method  : * the eddy diffusivity coef. specification depends on:
96      !!
97      !!    ln_traldf_lap = T     laplacian operator
98      !!    ln_traldf_blp = T   bilaplacian operator
99      !!
100      !!    nn_aht_ijk_t  =  0 => = constant
101      !!                  !
102      !!                  = 10 => = F(z) : constant with a reduction of 1/4 with depth
103      !!                  !
104      !!                  =-20 => = F(i,j)   = shape read in 'eddy_diffusivity.nc' file
105      !!                  = 20    = F(i,j)   = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
106      !!                  = 21    = F(i,j,t) = F(growth rate of baroclinic instability)
107      !!                  !
108      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_diffusivity.nc' file
109      !!                  = 30    = F(i,j,k)   = 2D (case 20) + decrease with depth (case 10)
110      !!                  = 31    = F(i,j,k,t) = F(local velocity) (  |u|e  /12   laplacian operator
111      !!                                                          or |u|e^3/12 bilaplacian operator )
112      !!              * initialisation of the eddy induced velocity coefficient by a call to ldf_eiv_init
113      !!           
114      !! ** action  : ahtu, ahtv initialized once for all or l_ldftra_time set to true
115      !!              aeiu, aeiv initialized once for all or l_ldfeiv_time set to true
116      !!----------------------------------------------------------------------
117      INTEGER  ::   jk                ! dummy loop indices
118      INTEGER  ::   ierr, inum, ios   ! local integer
119      REAL(wp) ::   zah0              ! local scalar
120      !!
121      NAMELIST/namtra_ldf/ ln_traldf_NONE, ln_traldf_lap  , ln_traldf_blp  ,  &   ! type of operator
122         &                 ln_traldf_lev , ln_traldf_hor  , ln_traldf_triad,  &   ! acting direction of the operator
123         &                 ln_traldf_iso , ln_traldf_msc  ,  rn_slpmax     ,  &   ! option for iso-neutral operator
124         &                 ln_triad_iso  , ln_botmix_triad, rn_sw_triad    ,  &   ! option for triad operator
125         &                 rn_aht_0      , rn_bht_0       , nn_aht_ijk_t          ! lateral eddy coefficient
126      !!----------------------------------------------------------------------
127      !
128      IF(lwp) THEN                      ! control print
129         WRITE(numout,*)
130         WRITE(numout,*) 'ldf_tra_init : lateral tracer physics'
131         WRITE(numout,*) '~~~~~~~~~~~~ '
132      ENDIF
133      !
134      !  Choice of lateral tracer physics
135      ! =================================
136      !
137      REWIND( numnam_ref )              ! Namelist namtra_ldf in reference namelist : Lateral physics on tracers
138      READ  ( numnam_ref, namtra_ldf, IOSTAT = ios, ERR = 901)
139901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_ldf in reference namelist', lwp )
140      REWIND( numnam_cfg )              ! Namelist namtra_ldf in configuration namelist : Lateral physics on tracers
141      READ  ( numnam_cfg, namtra_ldf, IOSTAT = ios, ERR = 902 )
142902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_ldf in configuration namelist', lwp )
143      IF(lwm) WRITE( numond, namtra_ldf )
144      !
145      IF(lwp) THEN                      ! control print
146         WRITE(numout,*) '   Namelist : namtra_ldf --- lateral mixing parameters (type, direction, coefficients)'
147         WRITE(numout,*) '      type :'
148         WRITE(numout,*) '         no explicit diffusion                   ln_traldf_NONE  = ', ln_traldf_NONE
149         WRITE(numout,*) '         laplacian operator                      ln_traldf_lap   = ', ln_traldf_lap
150         WRITE(numout,*) '         bilaplacian operator                    ln_traldf_blp   = ', ln_traldf_blp
151         WRITE(numout,*) '      direction of action :'
152         WRITE(numout,*) '         iso-level                               ln_traldf_lev   = ', ln_traldf_lev
153         WRITE(numout,*) '         horizontal (geopotential)               ln_traldf_hor   = ', ln_traldf_hor
154         WRITE(numout,*) '         iso-neutral Madec operator              ln_traldf_iso   = ', ln_traldf_iso
155         WRITE(numout,*) '         iso-neutral triad operator              ln_traldf_triad = ', ln_traldf_triad
156         WRITE(numout,*) '            iso-neutral (Method of Stab. Corr.)  ln_traldf_msc   = ', ln_traldf_msc
157         WRITE(numout,*) '            maximum isoppycnal slope             rn_slpmax       = ', rn_slpmax
158         WRITE(numout,*) '            pure lateral mixing in ML            ln_triad_iso    = ', ln_triad_iso
159         WRITE(numout,*) '            switching triad or not               rn_sw_triad     = ', rn_sw_triad
160         WRITE(numout,*) '            lateral mixing on bottom             ln_botmix_triad = ', ln_botmix_triad
161         WRITE(numout,*) '      coefficients :'
162         WRITE(numout,*) '         lateral eddy diffusivity   (lap case)   rn_aht_0        = ', rn_aht_0
163         WRITE(numout,*) '         lateral eddy diffusivity (bilap case)   rn_bht_0        = ', rn_bht_0
164         WRITE(numout,*) '         type of time-space variation            nn_aht_ijk_t    = ', nn_aht_ijk_t
165      ENDIF
166      !
167      !                                ! Parameter control
168      !
169      IF( ln_traldf_NONE ) THEN
170         IF(lwp) WRITE(numout,*) '   ==>>>   No diffusive operator selected. ahtu and ahtv are not allocated'
171         l_ldftra_time = .FALSE.
172         RETURN
173      ENDIF
174      !
175      IF( ln_traldf_blp .AND. ( ln_traldf_iso .OR. ln_traldf_triad) ) THEN     ! iso-neutral bilaplacian need MSC
176         IF( .NOT.ln_traldf_msc )   CALL ctl_stop( 'tra_ldf_init: iso-neutral bilaplacian requires ln_traldf_msc=.true.' )
177      ENDIF
178      !
179      !  Space/time variation of eddy coefficients
180      ! ===========================================
181      !                                               ! allocate the aht arrays
182      ALLOCATE( ahtu(jpi,jpj,jpk) , ahtv(jpi,jpj,jpk) , STAT=ierr )
183      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays')
184      !
185      ahtu(:,:,jpk) = 0._wp                           ! last level always 0 
186      ahtv(:,:,jpk) = 0._wp
187      !
188      !                                               ! value of eddy mixing coef.
189      IF    ( ln_traldf_lap ) THEN   ;   zah0 =      rn_aht_0        !   laplacian operator
190      ELSEIF( ln_traldf_blp ) THEN   ;   zah0 = ABS( rn_bht_0 )      ! bilaplacian operator
191      ENDIF
192      !
193      l_ldftra_time = .FALSE.                         ! no time variation except in case defined below
194      !
195      IF( ln_traldf_lap .OR. ln_traldf_blp ) THEN     ! only if a lateral diffusion operator is used
196         !
197         SELECT CASE(  nn_aht_ijk_t  )                   ! Specification of space time variations of ehtu, ahtv
198         !
199         CASE(   0  )      !==  constant  ==!
200            IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = constant = ', rn_aht_0
201            ahtu(:,:,:) = zah0 * umask(:,:,:)
202            ahtv(:,:,:) = zah0 * vmask(:,:,:)
203            !
204         CASE(  10  )      !==  fixed profile  ==!
205            IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( depth )'
206            ahtu(:,:,1) = zah0 * umask(:,:,1)                      ! constant surface value
207            ahtv(:,:,1) = zah0 * vmask(:,:,1)
208            CALL ldf_c1d( 'TRA', r1_4, ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv )
209            !
210         CASE ( -20 )      !== fixed horizontal shape read in file  ==!
211            IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F(i,j) read in eddy_diffusivity.nc file'
212            CALL iom_open( 'eddy_diffusivity_2D.nc', inum )
213            CALL iom_get ( inum, jpdom_data, 'ahtu_2D', ahtu(:,:,1) )
214            CALL iom_get ( inum, jpdom_data, 'ahtv_2D', ahtv(:,:,1) )
215            CALL iom_close( inum )
216            DO jk = 2, jpkm1
217               ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk)
218               ahtv(:,:,jk) = ahtv(:,:,1) * vmask(:,:,jk)
219            END DO
220            !
221         CASE(  20  )      !== fixed horizontal shape  ==!
222            IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or blp case)'
223            IF( ln_traldf_lap )   CALL ldf_c2d( 'TRA', 'LAP', zah0, ahtu, ahtv )    ! surface value proportional to scale factor
224            IF( ln_traldf_blp )   CALL ldf_c2d( 'TRA', 'BLP', zah0, ahtu, ahtv )    ! surface value proportional to scale factor
225            !
226         CASE(  21  )      !==  time varying 2D field  ==!
227            IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( latitude, longitude, time )'
228            IF(lwp) WRITE(numout,*) '                               = F( growth rate of baroclinic instability )'
229            IF(lwp) WRITE(numout,*) '                               min value = 0.1 * rn_aht_0'
230            IF(lwp) WRITE(numout,*) '                               max value = rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21)'
231            IF(lwp) WRITE(numout,*) '                               increased to rn_aht_0 within 20N-20S'
232            !
233            l_ldftra_time = .TRUE.     ! will be calculated by call to ldf_tra routine in step.F90
234            !
235            IF( ln_traldf_blp ) THEN
236               CALL ctl_stop( 'ldf_tra_init: aht=F(growth rate of baroc. insta.) incompatible with bilaplacian operator' )
237            ENDIF
238            !
239         CASE( -30  )      !== fixed 3D shape read in file  ==!
240            IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F(i,j,k) read in eddy_diffusivity.nc file'
241            CALL iom_open( 'eddy_diffusivity_3D.nc', inum )
242            CALL iom_get ( inum, jpdom_data, 'ahtu_3D', ahtu )
243            CALL iom_get ( inum, jpdom_data, 'ahtv_3D', ahtv )
244            CALL iom_close( inum )
245            DO jk = 1, jpkm1
246               ahtu(:,:,jk) = ahtu(:,:,jk) * umask(:,:,jk)
247               ahtv(:,:,jk) = ahtv(:,:,jk) * vmask(:,:,jk)
248            END DO
249            !
250         CASE(  30  )      !==  fixed 3D shape  ==!
251            IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( latitude, longitude, depth )'
252            IF( ln_traldf_lap )   CALL ldf_c2d( 'TRA', 'LAP', zah0, ahtu, ahtv )    ! surface value proportional to scale factor
253            IF( ln_traldf_blp )   CALL ldf_c2d( 'TRA', 'BLP', zah0, ahtu, ahtv )    ! surface value proportional to scale factor
254            !                                                    ! reduction with depth
255            CALL ldf_c1d( 'TRA', r1_4, ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv )
256            !
257         CASE(  31  )      !==  time varying 3D field  ==!
258            IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( latitude, longitude, depth , time )'
259            IF(lwp) WRITE(numout,*) '                                 proportional to the velocity : |u|e/12 or |u|e^3/12'
260            !
261            l_ldftra_time = .TRUE.     ! will be calculated by call to ldf_tra routine in step.F90
262            !
263         CASE DEFAULT
264            CALL ctl_stop('ldf_tra_init: wrong choice for nn_aht_ijk_t, the type of space-time variation of aht')
265         END SELECT
266         !
267         IF( ln_traldf_blp .AND. .NOT. l_ldftra_time ) THEN
268            ahtu(:,:,:) = SQRT( ahtu(:,:,:) )
269            ahtv(:,:,:) = SQRT( ahtv(:,:,:) )
270         ENDIF
271         !
272      ENDIF
273      !
274   END SUBROUTINE ldf_tra_init
275
276
277   SUBROUTINE ldf_tra( kt )
278      !!----------------------------------------------------------------------
279      !!                  ***  ROUTINE ldf_tra  ***
280      !!
281      !! ** Purpose :   update at kt the tracer lateral mixing coeff. (aht and aeiv)
282      !!
283      !! ** Method  :   time varying eddy diffusivity coefficients:
284      !!
285      !!    nn_aei_ijk_t = 21    aeiu, aeiv = F(i,j,  t) = F(growth rate of baroclinic instability)
286      !!                                                   with a reduction to 0 in vicinity of the Equator
287      !!    nn_aht_ijk_t = 21    ahtu, ahtv = F(i,j,  t) = F(growth rate of baroclinic instability)
288      !!
289      !!                 = 31    ahtu, ahtv = F(i,j,k,t) = F(local velocity) (  |u|e  /12   laplacian operator
290      !!                                                                     or |u|e^3/12 bilaplacian operator )
291      !!
292      !! ** action  :   ahtu, ahtv   update at each time step   
293      !!                aeiu, aeiv      -       -     -    -   (if ln_ldfeiv=T)
294      !!----------------------------------------------------------------------
295      INTEGER, INTENT(in) ::   kt   ! time step
296      !
297      INTEGER  ::   ji, jj, jk   ! dummy loop indices
298      REAL(wp) ::   zaht, zahf, zaht_min, z1_f20       ! local scalar
299      !!----------------------------------------------------------------------
300      !
301      IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN       ! eddy induced velocity coefficients
302         !                                ! =F(growth rate of baroclinic instability)
303         !                                ! max value rn_aeiv_0 ; decreased to 0 within 20N-20S
304         CALL ldf_eiv( kt, rn_aeiv_0, aeiu, aeiv )
305      ENDIF
306      !
307      SELECT CASE(  nn_aht_ijk_t  )       ! Eddy diffusivity coefficients
308      !
309      CASE(  21  )       !==  time varying 2D field  ==!   = F( growth rate of baroclinic instability )
310         !                                             !   min value rn_aht_0 / 10
311         !                                             !   max value rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21)
312         !                                             !   increase to rn_aht_0 within 20N-20S
313         IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN   ! use the already computed aei.
314            ahtu(:,:,1) = aeiu(:,:,1)
315            ahtv(:,:,1) = aeiv(:,:,1)
316         ELSE                                            ! compute aht.
317            CALL ldf_eiv( kt, rn_aht_0, ahtu, ahtv )
318         ENDIF
319         !
320         z1_f20   = 1._wp / (  2._wp * omega * SIN( rad * 20._wp )  )      ! 1 / ff(20 degrees)   
321         zaht_min = 0.2_wp * rn_aht_0                                      ! minimum value for aht
322         DO jj = 1, jpj
323            DO ji = 1, jpi
324               !!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg)
325               !!     ==>>>   The Coriolis value is identical for t- & u_points, and for v- and f-points
326               zaht = ( 1._wp -  MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min )
327               zahf = ( 1._wp -  MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min )
328               ahtu(ji,jj,1) = (  MAX( zaht_min, ahtu(ji,jj,1) ) + zaht  ) * umask(ji,jj,1)     ! min value zaht_min
329               ahtv(ji,jj,1) = (  MAX( zaht_min, ahtv(ji,jj,1) ) + zahf  ) * vmask(ji,jj,1)     ! increase within 20S-20N
330            END DO
331         END DO
332         DO jk = 2, jpkm1                             ! deeper value = surface value
333            ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk)
334            ahtv(:,:,jk) = ahtv(:,:,1) * vmask(:,:,jk)
335         END DO
336         !
337      CASE(  31  )       !==  time varying 3D field  ==!   = F( local velocity )
338         IF( ln_traldf_lap     ) THEN          !   laplacian operator |u| e /12
339            DO jk = 1, jpkm1
340               ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12
341               ahtv(:,:,jk) = ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12
342            END DO
343         ELSEIF( ln_traldf_blp ) THEN      ! bilaplacian operator      sqrt( |u| e^3 /12 ) = sqrt( |u| e /12 ) * e
344            DO jk = 1, jpkm1
345               ahtu(:,:,jk) = SQRT(  ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12  ) * e1u(:,:)
346               ahtv(:,:,jk) = SQRT(  ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12  ) * e2v(:,:)
347            END DO
348         ENDIF
349         !
350      END SELECT
351      !
352      CALL iom_put( "ahtu_2d", ahtu(:,:,1) )   ! surface u-eddy diffusivity coeff.
353      CALL iom_put( "ahtv_2d", ahtv(:,:,1) )   ! surface v-eddy diffusivity coeff.
354      CALL iom_put( "ahtu_3d", ahtu(:,:,:) )   ! 3D      u-eddy diffusivity coeff.
355      CALL iom_put( "ahtv_3d", ahtv(:,:,:) )   ! 3D      v-eddy diffusivity coeff.
356      !
357!!gm  : THE IF below is to be checked (comes from Seb)
358      IF( ln_ldfeiv ) THEN
359        CALL iom_put( "aeiu_2d", aeiu(:,:,1) )   ! surface u-EIV coeff.
360        CALL iom_put( "aeiv_2d", aeiv(:,:,1) )   ! surface v-EIV coeff.
361        CALL iom_put( "aeiu_3d", aeiu(:,:,:) )   ! 3D      u-EIV coeff.
362        CALL iom_put( "aeiv_3d", aeiv(:,:,:) )   ! 3D      v-EIV coeff.
363      ENDIF
364      !
365   END SUBROUTINE ldf_tra
366
367
368   SUBROUTINE ldf_eiv_init
369      !!----------------------------------------------------------------------
370      !!                  ***  ROUTINE ldf_eiv_init  ***
371      !!
372      !! ** Purpose :   initialization of the eiv coeff. from namelist choices.
373      !!
374      !! ** Method :
375      !!
376      !! ** Action :   aeiu , aeiv   : EIV coeff. at u- & v-points
377      !!               l_ldfeiv_time : =T if EIV coefficients vary with time
378      !!----------------------------------------------------------------------
379      INTEGER  ::   jk                ! dummy loop indices
380      INTEGER  ::   ierr, inum, ios   ! local integer
381      !
382      NAMELIST/namtra_ldfeiv/ ln_ldfeiv   , ln_ldfeiv_dia,   &    ! eddy induced velocity (eiv)
383         &                    nn_aei_ijk_t, rn_aeiv_0             ! eiv  coefficient
384      !!----------------------------------------------------------------------
385      !
386      IF(lwp) THEN                      ! control print
387         WRITE(numout,*)
388         WRITE(numout,*) 'ldf_eiv_init : eddy induced velocity parametrization'
389         WRITE(numout,*) '~~~~~~~~~~~~ '
390      ENDIF
391      !
392      REWIND( numnam_ref )              ! Namelist namtra_ldfeiv in reference namelist : eddy induced velocity param.
393      READ  ( numnam_ref, namtra_ldfeiv, IOSTAT = ios, ERR = 901)
394901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_ldfeiv in reference namelist', lwp )
395      !
396      REWIND( numnam_cfg )              ! Namelist namtra_ldfeiv in configuration namelist : eddy induced velocity param.
397      READ  ( numnam_cfg, namtra_ldfeiv, IOSTAT = ios, ERR = 902 )
398902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_ldfeiv in configuration namelist', lwp )
399      IF(lwm)  WRITE ( numond, namtra_ldfeiv )
400
401      IF(lwp) THEN                      ! control print
402         WRITE(numout,*) '   Namelist namtra_ldfeiv : '
403         WRITE(numout,*) '      Eddy Induced Velocity (eiv) param.      ln_ldfeiv     = ', ln_ldfeiv
404         WRITE(numout,*) '      eiv streamfunction & velocity diag.     ln_ldfeiv_dia = ', ln_ldfeiv_dia
405         WRITE(numout,*) '      eddy induced velocity coef.             rn_aeiv_0     = ', rn_aeiv_0
406         WRITE(numout,*) '      type of time-space variation            nn_aei_ijk_t  = ', nn_aei_ijk_t
407         WRITE(numout,*)
408      ENDIF
409      !
410      IF( ln_ldfeiv .AND. ln_traldf_blp )   CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' )
411
412      !                                 ! Parameter control
413      l_ldfeiv_time = .FALSE.   
414      !
415      IF( ln_ldfeiv ) THEN                         ! allocate the aei arrays
416         ALLOCATE( aeiu(jpi,jpj,jpk), aeiv(jpi,jpj,jpk), STAT=ierr )
417         IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ldf_eiv: failed to allocate arrays')
418         !
419         SELECT CASE( nn_aei_ijk_t )               ! Specification of space time variations of eaiu, aeiv
420         !
421         CASE(   0  )      !==  constant  ==!
422            IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity coef. = constant = ', rn_aeiv_0
423            aeiu(:,:,:) = rn_aeiv_0
424            aeiv(:,:,:) = rn_aeiv_0
425            !
426         CASE(  10  )      !==  fixed profile  ==!
427            IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity coef. = F( depth )'
428            aeiu(:,:,1) = rn_aeiv_0                                ! constant surface value
429            aeiv(:,:,1) = rn_aeiv_0
430            CALL ldf_c1d( 'TRA', r1_4, aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv )
431            !
432         CASE ( -20 )      !== fixed horizontal shape read in file  ==!
433            IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F(i,j) read in eddy_diffusivity_2D.nc file'
434            CALL iom_open ( 'eddy_induced_velocity_2D.nc', inum )
435            CALL iom_get  ( inum, jpdom_data, 'aeiu', aeiu(:,:,1) )
436            CALL iom_get  ( inum, jpdom_data, 'aeiv', aeiv(:,:,1) )
437            CALL iom_close( inum )
438            DO jk = 2, jpk
439               aeiu(:,:,jk) = aeiu(:,:,1)
440               aeiv(:,:,jk) = aeiv(:,:,1)
441            END DO
442            !
443         CASE(  20  )      !== fixed horizontal shape  ==!
444            IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or bilap case)'
445            CALL ldf_c2d( 'TRA', 'LAP', rn_aeiv_0, aeiu, aeiv )    ! surface value proportional to scale factor
446            !
447         CASE(  21  )       !==  time varying 2D field  ==!
448            IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( latitude, longitude, time )'
449            IF(lwp) WRITE(numout,*) '                               = F( growth rate of baroclinic instability )'
450            !
451            l_ldfeiv_time = .TRUE.     ! will be calculated by call to ldf_tra routine in step.F90
452            !
453         CASE( -30  )      !== fixed 3D shape read in file  ==!
454            IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'
455            CALL iom_open ( 'eddy_induced_velocity_3D.nc', inum )
456            CALL iom_get  ( inum, jpdom_data, 'aeiu', aeiu )
457            CALL iom_get  ( inum, jpdom_data, 'aeiv', aeiv )
458            CALL iom_close( inum )
459            !
460         CASE(  30  )       !==  fixed 3D shape  ==!
461            IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( latitude, longitude, depth )'
462            CALL ldf_c2d( 'TRA', 'LAP', rn_aeiv_0, aeiu, aeiv )    ! surface value proportional to scale factor
463            !                                                 ! reduction with depth
464            CALL ldf_c1d( 'TRA', r1_4, aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv )
465            !
466         CASE DEFAULT
467            CALL ctl_stop('ldf_tra_init: wrong choice for nn_aei_ijk_t, the type of space-time variation of aei')
468         END SELECT
469         !
470      ELSE
471          IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity param is NOT used neither diagnosed'
472          ln_ldfeiv_dia = .FALSE.
473      ENDIF
474      !                   
475   END SUBROUTINE ldf_eiv_init
476
477
478   SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv )
479      !!----------------------------------------------------------------------
480      !!                  ***  ROUTINE ldf_eiv  ***
481      !!
482      !! ** Purpose :   Compute the eddy induced velocity coefficient from the
483      !!              growth rate of baroclinic instability.
484      !!
485      !! ** Method  :   coefficient function of the growth rate of baroclinic instability
486      !!
487      !! Reference : Treguier et al. JPO 1997   ; Held and Larichev JAS 1996
488      !!----------------------------------------------------------------------
489      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index
490      REAL(wp)                        , INTENT(inout) ::   paei0          ! max value            [m2/s]
491      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   paeiu, paeiv   ! eiv coefficient      [m2/s]
492      !
493      INTEGER  ::   ji, jj, jk    ! dummy loop indices
494      REAL(wp) ::   zfw, ze3w, zn2, z1_f20, zaht, zaht_min, zzaei    ! local scalars
495      REAL(wp), DIMENSION(jpi,jpj) ::   zn, zah, zhw, zross, zaeiw   ! 2D workspace
496      !!----------------------------------------------------------------------
497      !
498      zn   (:,:) = 0._wp      ! Local initialization
499      zhw  (:,:) = 5._wp
500      zah  (:,:) = 0._wp
501      zross(:,:) = 0._wp
502      !                       ! Compute lateral diffusive coefficient at T-point
503      IF( ln_traldf_triad ) THEN
504         DO jk = 1, jpk
505            DO jj = 2, jpjm1
506               DO ji = 2, jpim1
507                  ! Take the max of N^2 and zero then take the vertical sum
508                  ! of the square root of the resulting N^2 ( required to compute
509                  ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
510                  zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
511                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk)
512                  ! Compute elements required for the inverse time scale of baroclinic
513                  ! eddies using the isopycnal slopes calculated in ldfslp.F :
514                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
515                  ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
516                  zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w
517                  zhw(ji,jj) = zhw(ji,jj) + ze3w
518               END DO
519            END DO
520         END DO
521      ELSE
522         DO jk = 1, jpk
523            DO jj = 2, jpjm1
524               DO ji = 2, jpim1
525                  ! Take the max of N^2 and zero then take the vertical sum
526                  ! of the square root of the resulting N^2 ( required to compute
527                  ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
528                  zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
529                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk)
530                  ! Compute elements required for the inverse time scale of baroclinic
531                  ! eddies using the isopycnal slopes calculated in ldfslp.F :
532                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
533                  ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
534                  zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
535                     &                            + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w
536                  zhw(ji,jj) = zhw(ji,jj) + ze3w
537               END DO
538            END DO
539         END DO
540      END IF
541
542      DO jj = 2, jpjm1
543         DO ji = fs_2, fs_jpim1   ! vector opt.
544            zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 )
545            ! Rossby radius at w-point taken < 40km and  > 2km
546            zross(ji,jj) = MAX( MIN( .4 * zn(ji,jj) / zfw, 40.e3 ), 2.e3 )
547            ! Compute aeiw by multiplying Ro^2 and T^-1
548            zaeiw(ji,jj) = zross(ji,jj) * zross(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)
549         END DO
550      END DO
551
552      !                                         !==  Bound on eiv coeff.  ==!
553      z1_f20 = 1._wp / (  2._wp * omega * sin( rad * 20._wp )  )
554      DO jj = 2, jpjm1
555         DO ji = fs_2, fs_jpim1   ! vector opt.
556            zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj)       ! tropical decrease
557            zaeiw(ji,jj) = MIN( zzaei , paei0 )                                  ! Max value = paei0
558         END DO
559      END DO
560      CALL lbc_lnk( zaeiw(:,:), 'W', 1. )       ! lateral boundary condition
561      !               
562      DO jj = 2, jpjm1                          !== aei at u- and v-points  ==!
563         DO ji = fs_2, fs_jpim1   ! vector opt.
564            paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj  ) ) * umask(ji,jj,1)
565            paeiv(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji  ,jj+1) ) * vmask(ji,jj,1)
566         END DO
567      END DO
568      CALL lbc_lnk_multi( paeiu(:,:,1), 'U', 1. , paeiv(:,:,1), 'V', 1. )      ! lateral boundary condition
569
570      DO jk = 2, jpkm1                          !==  deeper values equal the surface one  ==!
571         paeiu(:,:,jk) = paeiu(:,:,1) * umask(:,:,jk)
572         paeiv(:,:,jk) = paeiv(:,:,1) * vmask(:,:,jk)
573      END DO
574      ! 
575   END SUBROUTINE ldf_eiv
576
577
578   SUBROUTINE ldf_eiv_trp( kt, kit000, pun, pvn, pwn, cdtype )
579      !!----------------------------------------------------------------------
580      !!                  ***  ROUTINE ldf_eiv_trp  ***
581      !!
582      !! ** Purpose :   add to the input ocean transport the contribution of
583      !!              the eddy induced velocity parametrization.
584      !!
585      !! ** Method  :   The eddy induced transport is computed from a flux stream-
586      !!              function which depends on the slope of iso-neutral surfaces
587      !!              (see ldf_slp). For example, in the i-k plan :
588      !!                   psi_uw = mk(aeiu) e2u mi(wslpi)   [in m3/s]
589      !!                   Utr_eiv = - dk[psi_uw]
590      !!                   Vtr_eiv = + di[psi_uw]
591      !!                ln_ldfeiv_dia = T : output the associated streamfunction,
592      !!                                    velocity and heat transport (call ldf_eiv_dia)
593      !!
594      !! ** Action  : pun, pvn increased by the eiv transport
595      !!----------------------------------------------------------------------
596      INTEGER                         , INTENT(in   ) ::   kt       ! ocean time-step index
597      INTEGER                         , INTENT(in   ) ::   kit000   ! first time step index
598      CHARACTER(len=3)                , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
599      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun      ! in : 3 ocean transport components   [m3/s]
600      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pvn      ! out: 3 ocean transport components   [m3/s]
601      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pwn      ! increased by the eiv                [m3/s]
602      !!
603      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
604      REAL(wp) ::   zuwk, zuwk1, zuwi, zuwi1   ! local scalars
605      REAL(wp) ::   zvwk, zvwk1, zvwj, zvwj1   !   -      -
606      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpsi_uw, zpsi_vw
607      !!----------------------------------------------------------------------
608      !
609      IF( kt == kit000 )  THEN
610         IF(lwp) WRITE(numout,*)
611         IF(lwp) WRITE(numout,*) 'ldf_eiv_trp : eddy induced advection on ', cdtype,' :'
612         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component'
613      ENDIF
614
615     
616      zpsi_uw(:,:, 1 ) = 0._wp   ;   zpsi_vw(:,:, 1 ) = 0._wp
617      zpsi_uw(:,:,jpk) = 0._wp   ;   zpsi_vw(:,:,jpk) = 0._wp
618      !
619      DO jk = 2, jpkm1
620         DO jj = 1, jpjm1
621            DO ji = 1, fs_jpim1   ! vector opt.
622               zpsi_uw(ji,jj,jk) = - 0.25_wp * e2u(ji,jj) * ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk) )   &
623                  &                                       * ( aeiu (ji,jj,jk-1) + aeiu (ji  ,jj,jk) ) * umask(ji,jj,jk)
624               zpsi_vw(ji,jj,jk) = - 0.25_wp * e1v(ji,jj) * ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk) )   &
625                  &                                       * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj  ,jk) ) * vmask(ji,jj,jk)
626            END DO
627         END DO
628      END DO
629      !
630      DO jk = 1, jpkm1
631         DO jj = 1, jpjm1
632            DO ji = 1, fs_jpim1   ! vector opt.               
633               pun(ji,jj,jk) = pun(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
634               pvn(ji,jj,jk) = pvn(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
635            END DO
636         END DO
637      END DO
638      DO jk = 1, jpkm1
639         DO jj = 2, jpjm1
640            DO ji = fs_2, fs_jpim1   ! vector opt.
641               pwn(ji,jj,jk) = pwn(ji,jj,jk) + (  zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj  ,jk)   &
642                  &                             + zpsi_vw(ji,jj,jk) - zpsi_vw(ji  ,jj-1,jk) )
643            END DO
644         END DO
645      END DO
646      !
647      !                              ! diagnose the eddy induced velocity and associated heat transport
648      IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw )
649      !
650    END SUBROUTINE ldf_eiv_trp
651
652
653   SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw )
654      !!----------------------------------------------------------------------
655      !!                  ***  ROUTINE ldf_eiv_dia  ***
656      !!
657      !! ** Purpose :   diagnose the eddy induced velocity and its associated
658      !!              vertically integrated heat transport.
659      !!
660      !! ** Method :
661      !!
662      !!----------------------------------------------------------------------
663      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   psi_uw, psi_vw   ! streamfunction   [m3/s]
664      !
665      INTEGER  ::   ji, jj, jk    ! dummy loop indices
666      REAL(wp) ::   zztmp   ! local scalar
667      REAL(wp), DIMENSION(jpi,jpj)     ::   zw2d   ! 2D workspace
668      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d   ! 3D workspace
669      !!----------------------------------------------------------------------
670      !
671!!gm I don't like this routine....   Crazy  way of doing things, not optimal at all...
672!!gm     to be redesigned....   
673      !                                                  !==  eiv stream function: output  ==!
674      CALL lbc_lnk_multi( psi_uw, 'U', -1. , psi_vw, 'V', -1. )
675      !
676!!gm      CALL iom_put( "psi_eiv_uw", psi_uw )                 ! output
677!!gm      CALL iom_put( "psi_eiv_vw", psi_vw )
678      !
679      !                                                  !==  eiv velocities: calculate and output  ==!
680      !
681      zw3d(:,:,jpk) = 0._wp                                    ! bottom value always 0
682      !
683      DO jk = 1, jpkm1                                         ! e2u e3u u_eiv = -dk[psi_uw]
684         zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u_n(:,:,jk) )
685      END DO
686      CALL iom_put( "uoce_eiv", zw3d )
687      !
688      DO jk = 1, jpkm1                                         ! e1v e3v v_eiv = -dk[psi_vw]
689         zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v_n(:,:,jk) )
690      END DO
691      CALL iom_put( "voce_eiv", zw3d )
692      !
693      DO jk = 1, jpkm1                                         ! e1 e2 w_eiv = dk[psix] + dk[psix]
694         DO jj = 2, jpjm1
695            DO ji = fs_2, fs_jpim1  ! vector opt.
696               zw3d(ji,jj,jk) = (  psi_vw(ji,jj,jk) - psi_vw(ji  ,jj-1,jk)    &
697                  &              + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj  ,jk)  ) / e1e2t(ji,jj)
698            END DO
699         END DO
700      END DO
701      CALL lbc_lnk( zw3d, 'T', 1. )      ! lateral boundary condition
702      CALL iom_put( "woce_eiv", zw3d )
703      !
704      !
705      zztmp = 0.5_wp * rau0 * rcp 
706      IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN
707        zw2d(:,:)   = 0._wp 
708        zw3d(:,:,:) = 0._wp 
709        DO jk = 1, jpkm1
710           DO jj = 2, jpjm1
711              DO ji = fs_2, fs_jpim1   ! vector opt.
712                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
713                    &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji+1,jj,jk,jp_tem) ) 
714                 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
715              END DO
716           END DO
717        END DO
718        CALL lbc_lnk( zw2d, 'U', -1. )
719        CALL lbc_lnk( zw3d, 'U', -1. )
720        CALL iom_put( "ueiv_heattr"  , zztmp * zw2d )                  ! heat transport in i-direction
721        CALL iom_put( "ueiv_heattr3d", zztmp * zw3d )                  ! heat transport in i-direction
722      ENDIF
723      zw2d(:,:)   = 0._wp 
724      zw3d(:,:,:) = 0._wp 
725      DO jk = 1, jpkm1
726         DO jj = 2, jpjm1
727            DO ji = fs_2, fs_jpim1   ! vector opt.
728               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
729                  &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji,jj+1,jk,jp_tem) ) 
730               zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
731            END DO
732         END DO
733      END DO
734      CALL lbc_lnk( zw2d, 'V', -1. )
735      CALL iom_put( "veiv_heattr", zztmp * zw2d )                  !  heat transport in j-direction
736      CALL iom_put( "veiv_heattr", zztmp * zw3d )                  !  heat transport in j-direction
737      !
738      IF( ln_diaptr )  CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d )
739      !
740      zztmp = 0.5_wp * 0.5
741      IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN
742        zw2d(:,:) = 0._wp 
743        zw3d(:,:,:) = 0._wp 
744        DO jk = 1, jpkm1
745           DO jj = 2, jpjm1
746              DO ji = fs_2, fs_jpim1   ! vector opt.
747                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
748                    &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji+1,jj,jk,jp_sal) ) 
749                 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
750              END DO
751           END DO
752        END DO
753        CALL lbc_lnk( zw2d, 'U', -1. )
754        CALL lbc_lnk( zw3d, 'U', -1. )
755        CALL iom_put( "ueiv_salttr", zztmp * zw2d )                  ! salt transport in i-direction
756        CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                  ! salt transport in i-direction
757      ENDIF
758      zw2d(:,:) = 0._wp 
759      zw3d(:,:,:) = 0._wp 
760      DO jk = 1, jpkm1
761         DO jj = 2, jpjm1
762            DO ji = fs_2, fs_jpim1   ! vector opt.
763               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
764                  &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji,jj+1,jk,jp_sal) ) 
765               zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
766            END DO
767         END DO
768      END DO
769      CALL lbc_lnk( zw2d, 'V', -1. )
770      CALL iom_put( "veiv_salttr", zztmp * zw2d )                  !  salt transport in j-direction
771      CALL iom_put( "veiv_salttr", zztmp * zw3d )                  !  salt transport in j-direction
772      !
773      IF( ln_diaptr ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d )
774      !
775      !
776   END SUBROUTINE ldf_eiv_dia
777
778   !!======================================================================
779END MODULE ldftra
Note: See TracBrowser for help on using the repository browser.