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.
traldf_iso.F90 in branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.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: 23.5 KB
Line 
1MODULE traldf_iso
2   !!======================================================================
3   !!                   ***  MODULE  traldf_iso  ***
4   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
5   !!======================================================================
6   !! History :  OPA  ! 1994-08  (G. Madec, M. Imbard)
7   !!            8.0  ! 1997-05  (G. Madec)  split into traldf and trazdf
8   !!            NEMO ! 2002-08  (G. Madec)  Free form, F90
9   !!            1.0  ! 2005-11  (G. Madec)  merge traldf and trazdf :-)
10   !!            3.3  ! 2010-09  (C. Ethe, G. Madec) Merge TRA-TRC
11   !!            3.7  ! 2014-01  (G. Madec, S. Masson)  restructuration/simplification of aht/aeiv specification
12   !!             -   ! 2014-02  (F. Lemarie, G. Madec)  triad operator (Griffies) + Method of Stabilizing Correction
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   tra_ldf_iso  : update the tracer trend with the horizontal component of a iso-neutral laplacian operator
17   !!                  and with the vertical part of the isopycnal or geopotential s-coord. operator
18   !!----------------------------------------------------------------------
19   USE oce             ! ocean dynamics and active tracers
20   USE dom_oce         ! ocean space and time domain
21   USE trc_oce         ! share passive tracers/Ocean variables
22   USE zdf_oce         ! ocean vertical physics
23   USE ldftra          ! lateral diffusion: tracer eddy coefficients
24   USE ldfslp          ! iso-neutral slopes
25   USE diaptr          ! poleward transport diagnostics
26   !
27   USE in_out_manager  ! I/O manager
28   USE iom             ! I/O library
29   USE phycst          ! physical constants
30   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
31   USE wrk_nemo        ! Memory Allocation
32   USE timing          ! Timing
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   tra_ldf_iso   ! routine called by step.F90
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49  SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
50      &                                                   pgui, pgvi,   &
51      &                                       ptb , ptbb, pta , kjpt, kpass )
52      !!----------------------------------------------------------------------
53      !!                  ***  ROUTINE tra_ldf_iso  ***
54      !!
55      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
56      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
57      !!      add it to the general trend of tracer equation.
58      !!
59      !! ** Method  :   The horizontal component of the lateral diffusive trends
60      !!      is provided by a 2nd order operator rotated along neural or geopo-
61      !!      tential surfaces to which an eddy induced advection can be added
62      !!      It is computed using before fields (forward in time) and isopyc-
63      !!      nal or geopotential slopes computed in routine ldfslp.
64      !!
65      !!      1st part :  masked horizontal derivative of T  ( di[ t ] )
66      !!      ========    with partial cell update if ln_zps=T
67      !!                  with top     cell update if ln_isfcav
68      !!
69      !!      2nd part :  horizontal fluxes of the lateral mixing operator
70      !!      ========   
71      !!         zftu =  pahu e2u*e3u/e1u di[ tb ]
72      !!               - pahu e2u*uslp    dk[ mi(mk(tb)) ]
73      !!         zftv =  pahv e1v*e3v/e2v dj[ tb ]
74      !!               - pahv e2u*vslp    dk[ mj(mk(tb)) ]
75      !!      take the horizontal divergence of the fluxes:
76      !!         difft = 1/(e1e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  }
77      !!      Add this trend to the general trend (ta,sa):
78      !!         ta = ta + difft
79      !!
80      !!      3rd part: vertical trends of the lateral mixing operator
81      !!      ========  (excluding the vertical flux proportional to dk[t] )
82      !!      vertical fluxes associated with the rotated lateral mixing:
83      !!         zftw = - {  mi(mk(pahu)) * e2t*wslpi di[ mi(mk(tb)) ]
84      !!                   + mj(mk(pahv)) * e1t*wslpj dj[ mj(mk(tb)) ]  }
85      !!      take the horizontal divergence of the fluxes:
86      !!         difft = 1/(e1e2t*e3t) dk[ zftw ]
87      !!      Add this trend to the general trend (ta,sa):
88      !!         pta = pta + difft
89      !!
90      !! ** Action :   Update pta arrays with the before rotated diffusion
91      !!----------------------------------------------------------------------
92      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
93      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
94      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
95      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
96      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
97      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
98      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
99      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
100      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2)
101      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2)
102      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
103      !
104      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
105      INTEGER  ::  ierr             ! local integer
106      REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars
107      REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      -
108      REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      -
109#if defined key_diaar5
110      REAL(wp) ::   zztmp   ! local scalar
111#endif
112      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdkt, zdk1t, z2d
113      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdit, zdjt, zftu, zftv, ztfw 
114      !!----------------------------------------------------------------------
115      !
116      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_iso')
117      !
118      CALL wrk_alloc( jpi,jpj,       zdkt, zdk1t, z2d ) 
119      CALL wrk_alloc( jpi,jpj,jpk,   zdit, zdjt , zftu, zftv, ztfw  ) 
120      !
121      IF( kt == kit000 )  THEN
122         IF(lwp) WRITE(numout,*)
123         IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype
124         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
125         !
126         akz     (:,:,:) = 0._wp     
127         ah_wslp2(:,:,:) = 0._wp
128      ENDIF
129      !
130      !                                               ! set time step size (Euler/Leapfrog)
131      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2dt =     rdttra(1)      ! at nit000   (Euler)
132      ELSE                                        ;   z2dt = 2.* rdttra(1)      !             (Leapfrog)
133      ENDIF
134      z1_2dt = 1._wp / z2dt
135      !
136      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
137      ELSE                    ;   zsign = -1._wp
138      ENDIF
139         
140         
141      !!----------------------------------------------------------------------
142      !!   0 - calculate  ah_wslp2 and akz
143      !!----------------------------------------------------------------------
144      !
145      IF( kpass == 1 ) THEN                  !==  first pass only  ==!
146         !
147         DO jk = 2, jpkm1
148            DO jj = 2, jpjm1
149               DO ji = fs_2, fs_jpim1   ! vector opt.
150                  !
151                  zmsku = tmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
152                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
153                  zmskv = tmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
154                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
155                     !
156                  zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    &
157                     &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku
158                  zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    &
159                     &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv
160                     !
161                  ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
162                     &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk)
163               END DO
164            END DO
165         END DO
166         !
167         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
168            DO jk = 2, jpkm1
169               DO jj = 2, jpjm1
170                  DO ji = fs_2, fs_jpim1
171                     akz(ji,jj,jk) = 0.25_wp * (                                                                     &
172                        &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   &
173                        &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   &
174                        &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   &
175                        &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   )
176                  END DO
177               END DO
178            END DO
179            !
180            IF( ln_traldf_blp ) THEN                ! bilaplacian operator
181               DO jk = 2, jpkm1
182                  DO jj = 1, jpjm1
183                     DO ji = 1, fs_jpim1
184                        akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
185                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) )  )
186                     END DO
187                  END DO
188               END DO
189            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
190               DO jk = 2, jpkm1
191                  DO jj = 1, jpjm1
192                     DO ji = 1, fs_jpim1
193                        ze3w_2 = fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
194                        zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
195                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
196                     END DO
197                  END DO
198               END DO
199           ENDIF
200           !
201         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
202            akz(:,:,:) = ah_wslp2(:,:,:)     
203         ENDIF
204      ENDIF
205      !
206      !                                                          ! ===========
207      DO jn = 1, kjpt                                            ! tracer loop
208         !                                                       ! ===========
209         !                                               
210         !!----------------------------------------------------------------------
211         !!   I - masked horizontal derivative
212         !!----------------------------------------------------------------------
213!!gm : bug.... why (x,:,:)?   (1,jpj,:) and (jpi,1,:) should be sufficient....
214         zdit (1,:,:) = 0._wp     ;     zdit (jpi,:,:) = 0._wp
215         zdjt (1,:,:) = 0._wp     ;     zdjt (jpi,:,:) = 0._wp
216         !!end
217
218         ! Horizontal tracer gradient
219         DO jk = 1, jpkm1
220            DO jj = 1, jpjm1
221               DO ji = 1, fs_jpim1   ! vector opt.
222                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
223                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
224               END DO
225            END DO
226         END DO
227         IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient
228            DO jj = 1, jpjm1              ! bottom correction (partial bottom cell)
229               DO ji = 1, fs_jpim1   ! vector opt.
230!!gm  the following anonymous remark is to considered: ! IF useless if zpshde defines pgu everywhere
231                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)         
232                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
233               END DO
234            END DO
235            IF( ln_isfcav ) THEN      ! first wet level beneath a cavity
236               DO jj = 1, jpjm1
237                  DO ji = 1, fs_jpim1   ! vector opt.
238                     IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)         
239                     IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)     
240                  END DO
241               END DO
242            ENDIF
243         ENDIF
244
245         !!----------------------------------------------------------------------
246         !!   II - horizontal trend  (full)
247         !!----------------------------------------------------------------------
248         !
249         DO jk = 1, jpkm1                                 ! Horizontal slab
250            !
251            !                             !== Vertical tracer gradient
252            zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1
253            !
254            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)                          ! surface: zdkt(jk=1)=zdkt(jk=2)
255            ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk)
256            ENDIF
257!!gm I don't understand why we should need this.... since wmask is used instead of tmask
258!         IF ( ln_isfcav ) THEN
259!            DO jj = 1, jpj
260!               DO ji = 1, jpi   ! vector opt.
261!                  ikt = mikt(ji,jj) ! surface level
262!                  zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn  ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1)
263!                  zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt)
264!               END DO
265!            END DO
266!         END IF
267!!gm
268
269            DO jj = 1 , jpjm1            !==  Horizontal fluxes
270               DO ji = 1, fs_jpim1   ! vector opt.
271                  zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk)
272                  zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk)
273                  !
274                  zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   &
275                     &             + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk  ), 1. )
276                  !
277                  zmskv = 1. / MAX(  tmask(ji,jj+1,jk  ) + tmask(ji,jj,jk+1)   &
278                     &             + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1. )
279                  !
280                  zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
281                  zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
282                  !
283                  zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   &
284                     &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      &
285                     &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk)
286                  zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   &
287                     &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      &
288                     &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                 
289               END DO
290            END DO
291            !
292            DO jj = 2 , jpjm1          !== horizontal divergence and add to pta
293               DO ji = fs_2, fs_jpim1   ! vector opt.
294                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      &
295                     &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   &
296                     &                                        / (  e1e2t(ji,jj) * fse3t(ji,jj,jk)  )
297               END DO
298            END DO
299         END DO                                        !   End of slab 
300
301
302         !!----------------------------------------------------------------------
303         !!   III - vertical trend (full)
304         !!----------------------------------------------------------------------
305         
306         ztfw(1,:,:) = 0._wp     ;     ztfw(jpi,:,:) = 0._wp
307         
308         ! Vertical fluxes
309         ! ---------------
310         
311         ! Surface and bottom vertical fluxes set to zero
312         ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp
313         
314         ! interior (2=<jk=<jpk-1)
315         DO jk = 2, jpkm1
316            DO jj = 2, jpjm1
317               DO ji = fs_2, fs_jpim1   ! vector opt.
318                  !
319                  zmsku = tmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
320                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
321                  zmskv = tmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
322                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
323                     !
324                  zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    &
325                     &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku
326                  zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    &
327                     &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv
328                     !
329                  zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked
330                  zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
331                  !
332                  ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      &
333                     &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   &
334                     &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      &
335                     &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  )
336               END DO
337            END DO
338         END DO
339         !
340         !                                !==  add the vertical 33 flux  ==!
341         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
342            DO jk = 2, jpkm1       
343               DO jj = 1, jpjm1
344                  DO ji = fs_2, fs_jpim1
345                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)   &
346                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
347                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
348                  END DO
349               END DO
350            END DO
351            !
352         ELSE                                   ! bilaplacian
353            SELECT CASE( kpass )
354            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
355               DO jk = 2, jpkm1 
356                  DO jj = 1, jpjm1
357                     DO ji = fs_2, fs_jpim1
358                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    &
359                           &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   &
360                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) / fse3w(ji,jj,jk)
361                     END DO
362                  END DO
363               END DO
364            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp.
365               DO jk = 2, jpkm1 
366                  DO jj = 1, jpjm1
367                     DO ji = fs_2, fs_jpim1
368                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)                      &
369                           &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
370                           &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
371                     END DO
372                  END DO
373               END DO
374            END SELECT
375         ENDIF
376         !         
377         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==!
378            DO jj = 2, jpjm1
379               DO ji = fs_2, fs_jpim1   ! vector opt.
380                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   &
381                     &                                        / (  e1e2t(ji,jj) * fse3t_n(ji,jj,jk)  )
382               END DO
383            END DO
384         END DO
385         !
386         IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==!
387             ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==!
388            !
389            !                             ! "Poleward" diffusive heat or salt transports (T-S case only)
390            IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
391               ! note sign is reversed to give down-gradient diffusive transports (#1043)
392               IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( -zftv(:,:,:) )
393               IF( jn == jp_sal)   str_ldf(:) = ptr_sj( -zftv(:,:,:) )
394            ENDIF
395            !
396            IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN
397              !
398              IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN
399                  z2d(:,:) = zftu(ji,jj,1) 
400                  DO jk = 2, jpkm1
401                     DO jj = 2, jpjm1
402                        DO ji = fs_2, fs_jpim1   ! vector opt.
403                           z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
404                        END DO
405                     END DO
406                  END DO
407!!gm CAUTION I think there is an error of sign when using BLP operator....
408!!gm         a multiplication by zsign is required (to be checked twice !)
409                  z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
410                  CALL lbc_lnk( z2d, 'U', -1. )
411                  CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction
412                  !
413                  z2d(:,:) = zftv(ji,jj,1) 
414                  DO jk = 2, jpkm1
415                     DO jj = 2, jpjm1
416                        DO ji = fs_2, fs_jpim1   ! vector opt.
417                           z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
418                        END DO
419                     END DO
420                  END DO
421                  z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
422                  CALL lbc_lnk( z2d, 'V', -1. )
423                  CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction
424               END IF
425               !
426            ENDIF
427            !
428         ENDIF                                                    !== end pass selection  ==!
429         !
430         !                                                        ! ===============
431      END DO                                                      ! end tracer loop
432      !                                                           ! ===============
433      !
434      CALL wrk_dealloc( jpi, jpj,      zdkt, zdk1t, z2d ) 
435      CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt , zftu, zftv, ztfw  ) 
436      !
437      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_iso')
438      !
439   END SUBROUTINE tra_ldf_iso
440
441   !!==============================================================================
442END MODULE traldf_iso
Note: See TracBrowser for help on using the repository browser.