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/r12377_ticket2386/src/OCE/TRA – NEMO

source: NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traldf_triad.F90 @ 12451

Last change on this file since 12451 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
  • Property svn:keywords set to Id
File size: 21.5 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, z2dt, z1_2dt  !   -      -
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' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )      l_ptr = .TRUE. 
114      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
115         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
116      !
117      !                                                        ! set time step size (Euler/Leapfrog)
118      IF( neuler == 0 .AND. kt == kit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler)
119      ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog)
120      ENDIF
121      z1_2dt = 1._wp / z2dt
122      !
123      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
124      ELSE                    ;   zsign = -1._wp
125      ENDIF
126      !   
127      !!----------------------------------------------------------------------
128      !!   0 - calculate  ah_wslp2, akz, and optionally zpsi_uw, zpsi_vw
129      !!----------------------------------------------------------------------
130      !
131      IF( kpass == 1 ) THEN         !==  first pass only  and whatever the tracer is  ==!
132         !
133         akz     (:,:,:) = 0._wp     
134         ah_wslp2(:,:,:) = 0._wp
135         IF( ln_ldfeiv_dia ) THEN
136            zpsi_uw(:,:,:) = 0._wp
137            zpsi_vw(:,:,:) = 0._wp
138         ENDIF
139         !
140         DO ip = 0, 1                            ! i-k triads
141            DO kp = 0, 1
142               DO_3D_10_10( 1, jpkm1 )
143                  ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
144                  zbu   = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
145                  zah   = 0.25_wp * pahu(ji,jj,jk)
146                  zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
147                  ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth)
148                  zslope2 = zslope_skew + ( gdept(ji+1,jj,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
149                  zslope2 = zslope2 *zslope2
150                  ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2
151                  akz     (ji+ip,jj,jk+kp) = akz     (ji+ip,jj,jk+kp) + zah * r1_e1u(ji,jj)       &
152                     &                                                      * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
153                     !
154                 IF( ln_ldfeiv_dia )   zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp)   &
155                     &                                       + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * zslope_skew
156               END_3D
157            END DO
158         END DO
159         !
160         DO jp = 0, 1                            ! j-k triads
161            DO kp = 0, 1
162               DO_3D_10_10( 1, jpkm1 )
163                  ze3wr = 1.0_wp / e3w(ji,jj+jp,jk+kp,Kmm)
164                  zbv   = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
165                  zah   = 0.25_wp * pahv(ji,jj,jk)
166                  zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
167                  ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
168                  !    (do this by *adding* gradient of depth)
169                  zslope2 = zslope_skew + ( gdept(ji,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
170                  zslope2 = zslope2 * zslope2
171                  ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2
172                  akz     (ji,jj+jp,jk+kp) = akz     (ji,jj+jp,jk+kp) + zah * r1_e2v(ji,jj)     &
173                     &                                                      * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
174                  !
175                  IF( ln_ldfeiv_dia )   zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp)   &
176                     &                                       + 0.25 * aeiv(ji,jj,jk) * e1v(ji,jj) * zslope_skew
177               END_3D
178            END DO
179         END DO
180         !
181         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
182            !
183            IF( ln_traldf_blp ) THEN                ! bilaplacian operator
184               DO_3D_10_10( 2, jpkm1 )
185                  akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
186                     &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
187               END_3D
188            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
189               DO_3D_10_10( 2, jpkm1 )
190                  ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
191                  zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
192                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
193               END_3D
194           ENDIF
195           !
196         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
197            akz(:,:,:) = ah_wslp2(:,:,:)     
198         ENDIF
199         !
200         IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
201         !
202      ENDIF                                  !==  end 1st pass only  ==!
203      !
204      !                                                           ! ===========
205      DO jn = 1, kjpt                                             ! tracer loop
206         !                                                        ! ===========
207         ! Zero fluxes for each tracer
208!!gm  this should probably be done outside the jn loop
209         ztfw(:,:,:) = 0._wp
210         zftu(:,:,:) = 0._wp
211         zftv(:,:,:) = 0._wp
212         !
213         DO_3D_10_10( 1, jpkm1 )
214            zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
215            zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
216         END_3D
217         IF( ln_zps .AND. l_grad_zps ) THEN    ! partial steps: correction at top/bottom ocean level
218            DO_2D_10_10
219               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
220               zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
221            END_2D
222            IF( ln_isfcav ) THEN                   ! top level (ocean cavities only)
223               DO_2D_10_10
224                  IF( miku(ji,jj)  > 1 )   zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn) 
225                  IF( mikv(ji,jj)  > 1 )   zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn) 
226               END_2D
227            ENDIF
228         ENDIF
229         !
230         !!----------------------------------------------------------------------
231         !!   II - horizontal trend  (full)
232         !!----------------------------------------------------------------------
233         !
234         DO jk = 1, jpkm1
235            !                    !==  Vertical tracer gradient at level jk and jk+1
236            zdkt3d(:,:,1) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
237            !
238            !                    ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1)
239            IF( jk == 1 ) THEN   ;   zdkt3d(:,:,0) = zdkt3d(:,:,1)
240            ELSE                 ;   zdkt3d(:,:,0) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
241            ENDIF
242            !
243            zaei_slp = 0._wp
244            !
245            IF( ln_botmix_triad ) THEN
246               DO ip = 0, 1              !==  Horizontal & vertical fluxes
247                  DO kp = 0, 1
248                     DO_2D_10_10
249                        ze1ur = r1_e1u(ji,jj)
250                        zdxt  = zdit(ji,jj,jk) * ze1ur
251                        ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
252                        zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
253                        zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
254                        zslope_iso  = triadi  (ji+ip,jj,jk,1-ip,kp)
255                        !
256                        zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
257                        ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked....
258                        zah = pahu(ji,jj,jk)
259                        zah_slp  = zah * zslope_iso
260                        IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew
261                        zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
262                        ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - ( zah_slp + zaei_slp) * zdxt                 * zbu * ze3wr
263                     END_2D
264                  END DO
265               END DO
266               !
267               DO jp = 0, 1
268                  DO kp = 0, 1
269                     DO_2D_10_10
270                        ze2vr = r1_e2v(ji,jj)
271                        zdyt  = zdjt(ji,jj,jk) * ze2vr
272                        ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm)
273                        zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
274                        zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
275                        zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
276                        zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
277                        ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????  ahv is masked...
278                        zah = pahv(ji,jj,jk)
279                        zah_slp = zah * zslope_iso
280                        IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew
281                        zftv(ji,jj   ,jk   ) = zftv(ji,jj   ,jk   ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
282                        ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - ( zah_slp + zaei_slp ) * zdyt                * zbv * ze3wr
283                     END_2D
284                  END DO
285               END DO
286               !
287            ELSE
288               !
289               DO ip = 0, 1               !==  Horizontal & vertical fluxes
290                  DO kp = 0, 1
291                     DO_2D_10_10
292                        ze1ur = r1_e1u(ji,jj)
293                        zdxt  = zdit(ji,jj,jk) * ze1ur
294                        ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
295                        zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
296                        zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
297                        zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp)
298                        !
299                        zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
300                        ! ln_botmix_triad is .F. mask zah for bottom half cells
301                        zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ????
302                        zah_slp  = zah * zslope_iso
303                        IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew        ! aeit(ji+ip,jj,jk)*zslope_skew
304                        zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
305                        ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr
306                     END_2D
307                  END DO
308               END DO
309               !
310               DO jp = 0, 1
311                  DO kp = 0, 1
312                     DO_2D_10_10
313                        ze2vr = r1_e2v(ji,jj)
314                        zdyt  = zdjt(ji,jj,jk) * ze2vr
315                        ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm)
316                        zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
317                        zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
318                        zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
319                        zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
320                        ! ln_botmix_triad is .F. mask zah for bottom half cells
321                        zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! pahv(ji,jj+jp,jk)  ????
322                        zah_slp = zah * zslope_iso
323                        IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew        ! aeit(ji,jj+jp,jk)*zslope_skew
324                        zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
325                        ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr
326                     END_2D
327                  END DO
328               END DO
329            ENDIF
330            !                             !==  horizontal divergence and add to the general trend  ==!
331            DO_2D_00_00
332               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       &
333                  &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   &
334                  &                                        / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  )
335            END_2D
336            !
337         END DO
338         !
339         !                                !==  add the vertical 33 flux  ==!
340         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
341            DO_3D_10_00( 2, jpkm1 )
342               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   &
343                  &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
344                  &                            * (  pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
345            END_3D
346         ELSE                                   ! bilaplacian
347            SELECT CASE( kpass )
348            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
349               DO_3D_10_00( 2, jpkm1 )
350                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             &
351                     &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
352               END_3D
353            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
354               DO_3D_10_00( 2, jpkm1 )
355                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      &
356                     &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
357                     &                               + akz     (ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   )
358               END_3D
359            END SELECT
360         ENDIF
361         !
362         DO_3D_00_00( 1, jpkm1 )
363            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
364               &                                              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
365         END_3D
366         !
367         IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==!
368             ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==!
369            !
370            !                          ! "Poleward" diffusive heat or salt transports (T-S case only)
371            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', zftv(:,:,:)  )
372            !                          ! Diffusive heat transports
373            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', zftu(:,:,:), zftv(:,:,:) )
374            !
375         ENDIF                                                    !== end pass selection  ==!
376         !
377         !                                                        ! ===============
378      END DO                                                      ! end tracer loop
379      !                                                           ! ===============
380   END SUBROUTINE tra_ldf_triad
381
382   !!==============================================================================
383END MODULE traldf_triad
Note: See TracBrowser for help on using the repository browser.