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_triad.F90 in NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA – NEMO

source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/traldf_triad.F90 @ 12397

Last change on this file since 12397 was 12397, checked in by davestorkey, 4 years ago

2020/KERNEL-03_Storkey_Coward_RK3_stage2 : Consolidation of code to
handle initial Euler timestep in the context of leapfrog
timestepping. This version passes all SETTE tests but fails to bit
compare with the control for several tests (ORCA2_ICE_PISCES, AMM12,
ISOMIP, AGRIF_DEMO, SPITZ12).
  • Property svn:keywords set to Id
File size: 21.1 KB
Line 
1MODULE traldf_triad
2   !!======================================================================
3   !!                   ***  MODULE  traldf_triad  ***
4   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
5   !!======================================================================
6   !! History :  3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)  Griffies operator (original code)
7   !!            3.7  ! 2013-12  (F. Lemarie, G. Madec)  triad operator (Griffies) + Method of Stabilizing Correction
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   tra_ldf_triad : update the tracer trend with the iso-neutral laplacian triad-operator
12   !!----------------------------------------------------------------------
13   USE oce            ! ocean dynamics and active tracers
14   USE dom_oce        ! ocean space and time domain
15   USE phycst         ! physical constants
16   USE trc_oce        ! share passive tracers/Ocean variables
17   USE zdf_oce        ! ocean vertical physics
18   USE ldftra         ! lateral physics: eddy diffusivity
19   USE ldfslp         ! lateral physics: iso-neutral slopes
20   USE traldf_iso     ! lateral diffusion (Madec operator)         (tra_ldf_iso routine)
21   USE diaptr         ! poleward transport diagnostics
22   USE diaar5         ! AR5 diagnostics
23   USE zpshde         ! partial step: hor. derivative     (zps_hde routine)
24   !
25   USE in_out_manager ! I/O manager
26   USE iom            ! I/O library
27   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
28   USE lib_mpp        ! MPP library
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   tra_ldf_triad   ! routine called by traldf.F90
34
35   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   zdkt3d   !: vertical tracer gradient at 2 levels
36
37   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
38   LOGICAL  ::   l_hst   ! flag to compute heat transport
39
40
41   !! * Substitutions
42#  include "do_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
45   !! $Id$
46   !! Software governed by the CeCILL license (see ./LICENSE)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50  SUBROUTINE tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv,               &
51      &                                              pgu , pgv  , pgui, pgvi , &
52      &                                         pt , pt2, pt_rhs, kjpt, kpass )
53      !!----------------------------------------------------------------------
54      !!                  ***  ROUTINE tra_ldf_triad  ***
55      !!
56      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
57      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
58      !!      add it to the general trend of tracer equation.
59      !!
60      !! ** Method  :   The horizontal component of the lateral diffusive trends
61      !!      is provided by a 2nd order operator rotated along neural or geopo-
62      !!      tential surfaces to which an eddy induced advection can be added
63      !!      It is computed using before fields (forward in time) and isopyc-
64      !!      nal or geopotential slopes computed in routine ldfslp.
65      !!
66      !!      see documentation for the desciption
67      !!
68      !! ** Action :   pt_rhs   updated with the before rotated diffusion
69      !!               ah_wslp2 ....
70      !!               akz   stabilizing vertical diffusivity coefficient (used in trazdf_imp)
71      !!----------------------------------------------------------------------
72      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
73      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
74      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
75      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
76      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
77      INTEGER                              , INTENT(in)    ::   Kmm        ! ocean time level indices
78      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
79      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu , pgv  ! tracer gradient at pstep levels
80      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
81      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt         ! tracer (kpass=1) or laplacian of tracer (kpass=2)
82      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt2        ! tracer (only used in kpass=2)
83      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend
84      !
85      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
86      INTEGER  ::  ip,jp,kp         ! dummy loop indices
87      INTEGER  ::  ierr            ! local integer
88      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3    ! local scalars
89      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4    !   -      -
90      REAL(wp) ::  zcoef0, ze3w_2, zsign          !   -      -
91      !
92      REAL(wp) ::   zslope_skew, zslope_iso, zslope2, zbu, zbv
93      REAL(wp) ::   ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt
94      REAL(wp) ::   zah, zah_slp, zaei_slp
95      REAL(wp), DIMENSION(jpi,jpj    ) ::   z2d                                              ! 2D workspace
96      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw   ! 3D     -
97      !!----------------------------------------------------------------------
98      !
99      IF( .NOT.ALLOCATED(zdkt3d) )  THEN
100         ALLOCATE( zdkt3d(jpi,jpj,0:1) , STAT=ierr )
101         CALL mpp_sum ( 'traldf_triad', ierr )
102         IF( ierr > 0 )   CALL ctl_stop('STOP', 'tra_ldf_triad: unable to allocate arrays')
103      ENDIF
104     !
105      IF( kpass == 1 .AND. kt == kit000 )  THEN
106         IF(lwp) WRITE(numout,*)
107         IF(lwp) WRITE(numout,*) 'tra_ldf_triad : rotated laplacian diffusion operator on ', cdtype
108         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
109      ENDIF
110      !   
111      l_hst = .FALSE.
112      l_ptr = .FALSE.
113      IF( cdtype == 'TRA' ) THEN
114         IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf') )      l_ptr = .TRUE. 
115         IF( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR.                   &
116         &   iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  )   l_hst = .TRUE.
117      ENDIF
118      !
119      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
120      ELSE                    ;   zsign = -1._wp
121      ENDIF
122      !   
123      !!----------------------------------------------------------------------
124      !!   0 - calculate  ah_wslp2, akz, and optionally zpsi_uw, zpsi_vw
125      !!----------------------------------------------------------------------
126      !
127      IF( kpass == 1 ) THEN         !==  first pass only  and whatever the tracer is  ==!
128         !
129         akz     (:,:,:) = 0._wp     
130         ah_wslp2(:,:,:) = 0._wp
131         IF( ln_ldfeiv_dia ) THEN
132            zpsi_uw(:,:,:) = 0._wp
133            zpsi_vw(:,:,:) = 0._wp
134         ENDIF
135         !
136         DO ip = 0, 1                            ! i-k triads
137            DO kp = 0, 1
138               DO_3D_10_10( 1, jpkm1 )
139                  ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
140                  zbu   = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
141                  zah   = 0.25_wp * pahu(ji,jj,jk)
142                  zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
143                  ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth)
144                  zslope2 = zslope_skew + ( gdept(ji+1,jj,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
145                  zslope2 = zslope2 *zslope2
146                  ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2
147                  akz     (ji+ip,jj,jk+kp) = akz     (ji+ip,jj,jk+kp) + zah * r1_e1u(ji,jj)       &
148                     &                                                      * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
149                     !
150                 IF( ln_ldfeiv_dia )   zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp)   &
151                     &                                       + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * zslope_skew
152               END_3D
153            END DO
154         END DO
155         !
156         DO jp = 0, 1                            ! j-k triads
157            DO kp = 0, 1
158               DO_3D_10_10( 1, jpkm1 )
159                  ze3wr = 1.0_wp / e3w(ji,jj+jp,jk+kp,Kmm)
160                  zbv   = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
161                  zah   = 0.25_wp * pahv(ji,jj,jk)
162                  zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
163                  ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
164                  !    (do this by *adding* gradient of depth)
165                  zslope2 = zslope_skew + ( gdept(ji,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
166                  zslope2 = zslope2 * zslope2
167                  ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2
168                  akz     (ji,jj+jp,jk+kp) = akz     (ji,jj+jp,jk+kp) + zah * r1_e2v(ji,jj)     &
169                     &                                                      * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
170                  !
171                  IF( ln_ldfeiv_dia )   zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp)   &
172                     &                                       + 0.25 * aeiv(ji,jj,jk) * e1v(ji,jj) * zslope_skew
173               END_3D
174            END DO
175         END DO
176         !
177         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
178            !
179            IF( ln_traldf_blp ) THEN                ! bilaplacian operator
180               DO_3D_10_10( 2, jpkm1 )
181                  akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
182                     &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
183               END_3D
184            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
185               DO_3D_10_10( 2, jpkm1 )
186                  ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
187                  zcoef0 = r2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
188                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_2dt
189               END_3D
190           ENDIF
191           !
192         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
193            akz(:,:,:) = ah_wslp2(:,:,:)     
194         ENDIF
195         !
196         IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
197         !
198      ENDIF                                  !==  end 1st pass only  ==!
199      !
200      !                                                           ! ===========
201      DO jn = 1, kjpt                                             ! tracer loop
202         !                                                        ! ===========
203         ! Zero fluxes for each tracer
204!!gm  this should probably be done outside the jn loop
205         ztfw(:,:,:) = 0._wp
206         zftu(:,:,:) = 0._wp
207         zftv(:,:,:) = 0._wp
208         !
209         DO_3D_10_10( 1, jpkm1 )
210            zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
211            zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
212         END_3D
213         IF( ln_zps .AND. l_grad_zps ) THEN    ! partial steps: correction at top/bottom ocean level
214            DO_2D_10_10
215               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
216               zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
217            END_2D
218            IF( ln_isfcav ) THEN                   ! top level (ocean cavities only)
219               DO_2D_10_10
220                  IF( miku(ji,jj)  > 1 )   zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn) 
221                  IF( mikv(ji,jj)  > 1 )   zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn) 
222               END_2D
223            ENDIF
224         ENDIF
225         !
226         !!----------------------------------------------------------------------
227         !!   II - horizontal trend  (full)
228         !!----------------------------------------------------------------------
229         !
230         DO jk = 1, jpkm1
231            !                    !==  Vertical tracer gradient at level jk and jk+1
232            zdkt3d(:,:,1) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
233            !
234            !                    ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1)
235            IF( jk == 1 ) THEN   ;   zdkt3d(:,:,0) = zdkt3d(:,:,1)
236            ELSE                 ;   zdkt3d(:,:,0) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
237            ENDIF
238            !
239            zaei_slp = 0._wp
240            !
241            IF( ln_botmix_triad ) THEN
242               DO ip = 0, 1              !==  Horizontal & vertical fluxes
243                  DO kp = 0, 1
244                     DO_2D_10_10
245                        ze1ur = r1_e1u(ji,jj)
246                        zdxt  = zdit(ji,jj,jk) * ze1ur
247                        ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
248                        zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
249                        zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
250                        zslope_iso  = triadi  (ji+ip,jj,jk,1-ip,kp)
251                        !
252                        zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
253                        ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked....
254                        zah = pahu(ji,jj,jk)
255                        zah_slp  = zah * zslope_iso
256                        IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew
257                        zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
258                        ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - ( zah_slp + zaei_slp) * zdxt                 * zbu * ze3wr
259                     END_2D
260                  END DO
261               END DO
262               !
263               DO jp = 0, 1
264                  DO kp = 0, 1
265                     DO_2D_10_10
266                        ze2vr = r1_e2v(ji,jj)
267                        zdyt  = zdjt(ji,jj,jk) * ze2vr
268                        ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm)
269                        zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
270                        zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
271                        zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
272                        zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
273                        ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????  ahv is masked...
274                        zah = pahv(ji,jj,jk)
275                        zah_slp = zah * zslope_iso
276                        IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew
277                        zftv(ji,jj   ,jk   ) = zftv(ji,jj   ,jk   ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
278                        ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - ( zah_slp + zaei_slp ) * zdyt                * zbv * ze3wr
279                     END_2D
280                  END DO
281               END DO
282               !
283            ELSE
284               !
285               DO ip = 0, 1               !==  Horizontal & vertical fluxes
286                  DO kp = 0, 1
287                     DO_2D_10_10
288                        ze1ur = r1_e1u(ji,jj)
289                        zdxt  = zdit(ji,jj,jk) * ze1ur
290                        ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
291                        zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
292                        zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
293                        zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp)
294                        !
295                        zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
296                        ! ln_botmix_triad is .F. mask zah for bottom half cells
297                        zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ????
298                        zah_slp  = zah * zslope_iso
299                        IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew        ! aeit(ji+ip,jj,jk)*zslope_skew
300                        zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
301                        ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr
302                     END_2D
303                  END DO
304               END DO
305               !
306               DO jp = 0, 1
307                  DO kp = 0, 1
308                     DO_2D_10_10
309                        ze2vr = r1_e2v(ji,jj)
310                        zdyt  = zdjt(ji,jj,jk) * ze2vr
311                        ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm)
312                        zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
313                        zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
314                        zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
315                        zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
316                        ! ln_botmix_triad is .F. mask zah for bottom half cells
317                        zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! pahv(ji,jj+jp,jk)  ????
318                        zah_slp = zah * zslope_iso
319                        IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew        ! aeit(ji,jj+jp,jk)*zslope_skew
320                        zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
321                        ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr
322                     END_2D
323                  END DO
324               END DO
325            ENDIF
326            !                             !==  horizontal divergence and add to the general trend  ==!
327            DO_2D_00_00
328               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       &
329                  &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   &
330                  &                                        / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  )
331            END_2D
332            !
333         END DO
334         !
335         !                                !==  add the vertical 33 flux  ==!
336         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
337            DO_3D_10_00( 2, jpkm1 )
338               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   &
339                  &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
340                  &                            * (  pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
341            END_3D
342         ELSE                                   ! bilaplacian
343            SELECT CASE( kpass )
344            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
345               DO_3D_10_00( 2, jpkm1 )
346                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             &
347                     &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
348               END_3D
349            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
350               DO_3D_10_00( 2, jpkm1 )
351                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      &
352                     &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
353                     &                               + akz     (ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   )
354               END_3D
355            END SELECT
356         ENDIF
357         !
358         DO_3D_00_00( 1, jpkm1 )
359            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
360               &                                              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
361         END_3D
362         !
363         IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==!
364             ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==!
365            !
366            !                          ! "Poleward" diffusive heat or salt transports (T-S case only)
367            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', zftv(:,:,:)  )
368            !                          ! Diffusive heat transports
369            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', zftu(:,:,:), zftv(:,:,:) )
370            !
371         ENDIF                                                    !== end pass selection  ==!
372         !
373         !                                                        ! ===============
374      END DO                                                      ! end tracer loop
375      !                                                           ! ===============
376   END SUBROUTINE tra_ldf_triad
377
378   !!==============================================================================
379END MODULE traldf_triad
Note: See TracBrowser for help on using the repository browser.