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/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

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