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/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90 @ 7508

Last change on this file since 7508 was 7508, checked in by mocavero, 7 years ago

changes on code duplication and workshare construct
  • Property svn:keywords set to Id
File size: 25.2 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 "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
43   !! $Id$
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48  SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
49      &                                                   pgui, pgvi,   &
50      &                                       ptb , ptbb, pta , kjpt, kpass )
51      !!----------------------------------------------------------------------
52      !!                  ***  ROUTINE tra_ldf_iso  ***
53      !!
54      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
55      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
56      !!      add it to the general trend of tracer equation.
57      !!
58      !! ** Method  :   The horizontal component of the lateral diffusive trends
59      !!      is provided by a 2nd order operator rotated along neural or geopo-
60      !!      tential surfaces to which an eddy induced advection can be added
61      !!      It is computed using before fields (forward in time) and isopyc-
62      !!      nal or geopotential slopes computed in routine ldfslp.
63      !!
64      !!      1st part :  masked horizontal derivative of T  ( di[ t ] )
65      !!      ========    with partial cell update if ln_zps=T
66      !!                  with top     cell update if ln_isfcav
67      !!
68      !!      2nd part :  horizontal fluxes of the lateral mixing operator
69      !!      ========   
70      !!         zftu =  pahu e2u*e3u/e1u di[ tb ]
71      !!               - pahu e2u*uslp    dk[ mi(mk(tb)) ]
72      !!         zftv =  pahv e1v*e3v/e2v dj[ tb ]
73      !!               - pahv e2u*vslp    dk[ mj(mk(tb)) ]
74      !!      take the horizontal divergence of the fluxes:
75      !!         difft = 1/(e1e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  }
76      !!      Add this trend to the general trend (ta,sa):
77      !!         ta = ta + difft
78      !!
79      !!      3rd part: vertical trends of the lateral mixing operator
80      !!      ========  (excluding the vertical flux proportional to dk[t] )
81      !!      vertical fluxes associated with the rotated lateral mixing:
82      !!         zftw = - {  mi(mk(pahu)) * e2t*wslpi di[ mi(mk(tb)) ]
83      !!                   + mj(mk(pahv)) * e1t*wslpj dj[ mj(mk(tb)) ]  }
84      !!      take the horizontal divergence of the fluxes:
85      !!         difft = 1/(e1e2t*e3t) dk[ zftw ]
86      !!      Add this trend to the general trend (ta,sa):
87      !!         pta = pta + difft
88      !!
89      !! ** Action :   Update pta arrays with the before rotated diffusion
90      !!----------------------------------------------------------------------
91      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
92      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
93      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
94      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
95      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
96      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
97      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
98      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
99      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2)
100      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2)
101      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
102      !
103      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
104      INTEGER  ::  ikt
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!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
127         DO jk = 1, jpk
128            DO jj = 1, jpj
129               DO ji = 1, jpi
130                  akz     (ji,jj,jk) = 0._wp     
131                  ah_wslp2(ji,jj,jk) = 0._wp
132               END DO
133            END DO
134         END DO
135      ENDIF
136      !                                               ! set time step size (Euler/Leapfrog)
137      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler)
138      ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog)
139      ENDIF
140      z1_2dt = 1._wp / z2dt
141      !
142      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
143      ELSE                    ;   zsign = -1._wp
144      ENDIF
145         
146      !!----------------------------------------------------------------------
147      !!   0 - calculate  ah_wslp2 and akz
148      !!----------------------------------------------------------------------
149      !
150      IF( kpass == 1 ) THEN                  !==  first pass only  ==!
151         !
152!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zmsku, zmskv, zahu_w, zahv_w)
153         DO jk = 2, jpkm1
154            DO jj = 2, jpjm1
155               DO ji = fs_2, fs_jpim1   ! vector opt.
156                  !
157                  zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
158                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
159                  zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
160                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
161                     !
162                  zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    &
163                     &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku
164                  zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    &
165                     &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv
166                     !
167                  ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
168                     &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk)
169               END DO
170            END DO
171         END DO
172         !
173         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
174!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
175            DO jk = 2, jpkm1
176               DO jj = 2, jpjm1
177                  DO ji = fs_2, fs_jpim1
178                     akz(ji,jj,jk) = 0.25_wp * (                                                                     &
179                        &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   &
180                        &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   &
181                        &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   &
182                        &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   )
183                  END DO
184               END DO
185            END DO
186            !
187            IF( ln_traldf_blp ) THEN                ! bilaplacian operator
188!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
189               DO jk = 2, jpkm1
190                  DO jj = 1, jpjm1
191                     DO ji = 1, fs_jpim1
192                        akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
193                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  )
194                     END DO
195                  END DO
196               END DO
197            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
198!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, ze3w_2, zcoef0)
199               DO jk = 2, jpkm1
200                  DO jj = 1, jpjm1
201                     DO ji = 1, fs_jpim1
202                        ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk)
203                        zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
204                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
205                     END DO
206                  END DO
207               END DO
208           ENDIF
209           !
210         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
211!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
212           DO jk = 1, jpk
213              DO jj = 1, jpj
214                 DO ji = 1, jpi
215                    akz     (ji,jj,jk) = ah_wslp2(ji,jj,jk)
216                 END DO
217              END DO
218           END DO
219         ENDIF
220      ENDIF
221      !
222      !                                                          ! ===========
223      DO jn = 1, kjpt                                            ! tracer loop
224         !                                                       ! ===========
225         !                                               
226         !!----------------------------------------------------------------------
227         !!   I - masked horizontal derivative
228         !!----------------------------------------------------------------------
229!!gm : bug.... why (x,:,:)?   (1,jpj,:) and (jpi,1,:) should be sufficient....
230!$OMP PARALLEL
231!$OMP DO schedule(static) private(jk, jj)
232          DO jk = 1, jpk
233             DO jj = 1, jpj
234                zdit (1,jj,jk) = 0._wp     ;     zdit (jpi,jj,jk) = 0._wp
235                zdjt (1,jj,jk) = 0._wp     ;     zdjt (jpi,jj,jk) = 0._wp
236             END DO
237          END DO
238         !!end
239
240         ! Horizontal tracer gradient
241!$OMP DO schedule(static) private(jk, jj, ji)
242         DO jk = 1, jpkm1
243            DO jj = 1, jpjm1
244               DO ji = 1, fs_jpim1   ! vector opt.
245                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
246                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
247               END DO
248            END DO
249         END DO
250!$OMP END DO NOWAIT
251!$OMP END PARALLEL
252         IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient
253!$OMP PARALLEL DO schedule(static) private(jj, ji)
254            DO jj = 1, jpjm1              ! bottom correction (partial bottom cell)
255               DO ji = 1, fs_jpim1   ! vector opt.
256                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)         
257                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
258               END DO
259            END DO
260            IF( ln_isfcav ) THEN      ! first wet level beneath a cavity
261!$OMP PARALLEL DO schedule(static) private(jj, ji)
262               DO jj = 1, jpjm1
263                  DO ji = 1, fs_jpim1   ! vector opt.
264                     IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)         
265                     IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)     
266                  END DO
267               END DO
268            ENDIF
269         ENDIF
270         !
271         !!----------------------------------------------------------------------
272         !!   II - horizontal trend  (full)
273         !!----------------------------------------------------------------------
274!$OMP PARALLEL
275            DO jk = 1, jpkm1                                 ! Horizontal slab
276!$OMP DO schedule(static) private(jj, ji)
277               DO jj = 1 , jpj
278                  DO ji = 1, jpi
279                     zdk1t(ji,jj) = ( ptb(ji,jj,jk,jn) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1)
280                     IF( jk == 1 ) THEN   ;   zdkt(ji,jj) = zdk1t(ji,jj)               ! surface: zdkt(jk=1)=zdkt(jk=2)
281                     ELSE                 ;   zdkt(ji,jj) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
282                     ENDIF
283                  END DO
284               END DO
285!$OMP DO schedule(static) private(jj, ji, zmsku, zmskv, zabe1, zabe2, zcof1, zcof2)
286               DO jj = 1 , jpjm1            !==  Horizontal fluxes
287                  DO ji = 1, fs_jpim1   ! vector opt.
288                        zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
289                        zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
290                        !
291                        zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   &
292                        &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1.)
293                        !
294                        zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   &
295                        &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1.)
296                        !
297                        zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
298                        zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
299                  !
300                        zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   &
301                     &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj) &
302                     &                          + zdk1t(ji+1,jj) + zdkt (ji,jj))  ) * umask(ji,jj,jk)
303                        zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   &
304                     &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj) &
305                     &                          + zdk1t(ji,jj+1) + zdkt (ji,jj))  ) * vmask(ji,jj,jk)
306                  END DO
307               END DO
308            !
309!$OMP DO schedule(static) private(jj, ji)
310               DO jj = 2 , jpjm1          !== horizontal divergence and add to pta
311                  DO ji = fs_2, fs_jpim1   ! vector opt.
312                        pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      &
313                     &               + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   &
314                     &                          * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
315                  END DO
316               END DO
317!$OMP END DO NOWAIT
318            END DO
319
320
321         !!----------------------------------------------------------------------
322         !!   III - vertical trend (full)
323         !!----------------------------------------------------------------------
324         !
325!$OMP DO schedule(static) private(jk, jj)
326         DO jk = 1, jpk
327             DO jj = 1, jpj
328                ztfw(1,jj,jk) = 0._wp     ;     ztfw(jpi,jj,jk) = 0._wp
329             END DO
330          END DO
331         !
332         ! Vertical fluxes
333         ! ---------------
334         !                          ! Surface and bottom vertical fluxes set to zero
335!$OMP DO schedule(static) private(jj, ji)
336         DO jj = 1, jpj
337             DO ji = 1, jpi
338                ztfw(ji,jj, 1 ) = 0._wp      ;      ztfw(ji,jj,jpk) = 0._wp
339             END DO
340          END DO
341         
342!$OMP DO schedule(static) private(jk, jj, ji, zmsku, zmskv, zahu_w, zahv_w, zcoef3, zcoef4)
343         DO jk = 2, jpkm1           ! interior (2=<jk=<jpk-1)
344            DO jj = 2, jpjm1
345               DO ji = fs_2, fs_jpim1   ! vector opt.
346                  !
347                  zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
348                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
349                  zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
350                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
351                     !
352                  zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    &
353                     &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku
354                  zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    &
355                     &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv
356                     !
357                  zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked
358                  zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
359                  !
360                  ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      &
361                     &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   &
362                     &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      &
363                     &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  )
364               END DO
365            END DO
366         END DO
367!$OMP END DO NOWAIT
368!$OMP END PARALLEL
369         !                                !==  add the vertical 33 flux  ==!
370         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
371!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
372            DO jk = 2, jpkm1       
373               DO jj = 1, jpjm1
374                  DO ji = fs_2, fs_jpim1
375                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   &
376                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
377                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
378                  END DO
379               END DO
380            END DO
381            !
382         ELSE                                   ! bilaplacian
383            SELECT CASE( kpass )
384            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
385!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
386               DO jk = 2, jpkm1 
387                  DO jj = 1, jpjm1
388                     DO ji = fs_2, fs_jpim1
389                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    &
390                           &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   &
391                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
392                     END DO
393                  END DO
394               END DO
395            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp.
396!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
397               DO jk = 2, jpkm1 
398                  DO jj = 1, jpjm1
399                     DO ji = fs_2, fs_jpim1
400                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      &
401                           &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
402                           &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
403                     END DO
404                  END DO
405               END DO
406            END SELECT
407         ENDIF
408         !         
409!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
410         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==!
411            DO jj = 2, jpjm1
412               DO ji = fs_2, fs_jpim1   ! vector opt.
413                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   &
414                     &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
415               END DO
416            END DO
417         END DO
418         !
419         IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==!
420             ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==!
421            !
422            !                             ! "Poleward" diffusive heat or salt transports (T-S case only)
423            IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
424               ! note sign is reversed to give down-gradient diffusive transports (#1043)
425               IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( -zftv(:,:,:) )
426               IF( jn == jp_sal)   str_ldf(:) = ptr_sj( -zftv(:,:,:) )
427            ENDIF
428            !
429            IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN
430              !
431              IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN
432                  z2d(:,:) = zftu(ji,jj,1) 
433!$OMP PARALLEL
434!$OMP DO schedule(static) private(jk, jj, ji)
435                  DO jk = 2, jpkm1
436                     DO jj = 2, jpjm1
437                        DO ji = fs_2, fs_jpim1   ! vector opt.
438                           z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
439                        END DO
440                     END DO
441                  END DO
442!!gm CAUTION I think there is an error of sign when using BLP operator....
443!!gm         a multiplication by zsign is required (to be checked twice !)
444!$OMP DO schedule(static) private(jj, ji)
445         DO jj = 1, jpj
446             DO ji = 1, jpi
447                  z2d(ji,jj) = - rau0_rcp * z2d(ji,jj)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
448             END DO
449         END DO
450!$OMP END PARALLEL
451                  CALL lbc_lnk( z2d, 'U', -1. )
452                  CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction
453                  !
454                  z2d(:,:) = zftv(ji,jj,1) 
455!$OMP PARALLEL
456!$OMP DO schedule(static) private(jk, jj, ji)
457                  DO jk = 2, jpkm1
458                     DO jj = 2, jpjm1
459                        DO ji = fs_2, fs_jpim1   ! vector opt.
460                           z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
461                        END DO
462                     END DO
463                  END DO
464!$OMP DO schedule(static) private(jj, ji)
465         DO jj = 1, jpj
466             DO ji = 1, jpi
467                  z2d(ji,jj) = - rau0_rcp * z2d(ji,jj)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
468             END DO
469         END DO
470!$OMP END PARALLEL
471                  CALL lbc_lnk( z2d, 'V', -1. )
472                  CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction
473               END IF
474               !
475            ENDIF
476            !
477         ENDIF                                                    !== end pass selection  ==!
478         !
479         !                                                        ! ===============
480      END DO                                                      ! end tracer loop
481      !                                                           ! ===============
482      !
483      CALL wrk_dealloc( jpi, jpj,      zdkt, zdk1t, z2d ) 
484      CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt , zftu, zftv, ztfw  ) 
485      !
486      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_iso')
487      !
488   END SUBROUTINE tra_ldf_iso
489
490   !!==============================================================================
491END MODULE traldf_iso
Note: See TracBrowser for help on using the repository browser.