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_tile.F90 in NEMO/branches/UKMO/BPC_miniapp/Master – NEMO

source: NEMO/branches/UKMO/BPC_miniapp/Master/traldf_iso_tile.F90 @ 10831

Last change on this file since 10831 was 10831, checked in by wayne_gaudin, 5 years ago

Ticket #2197 - added the Master directory with code and makefile

File size: 23.7 KB
Line 
1MODULE traldf_iso_tile
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_tile   : 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 par_kind
20   USE len_oce            ! ocean dynamics and active tracers
21   USE tiletype
22   !
23   USE in_out_manager ! I/O manager
24
25   IMPLICIT NONE
26
27   PUBLIC   tra_ldf_iso_tile   ! routine called by step.F90
28
29   !! * Substitutions
30#  include "vectopt_loop_substitute.h90"
31   !!----------------------------------------------------------------------
32   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
33   !! $Id$
34   !! Software governed by the CeCILL licence     (./LICENSE)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38  SUBROUTINE tra_ldf_iso_tile( tile, ln_traldf_msc, ln_traldf_blp, ln_traldf_lap, ln_zps, ln_isfcav, kjpt, kpass, p2dt, &
39      &                   e1t, e1u, e1v, e2t, e2u, e2v, e1e2t, e1_e2v, e2_e1u, r1_e1e2t, mbku, mbkv, miku, mikv, &
40      &                   umask, vmask, wmask, e3t_n, e3u_n, e3v_n, e3w_n, wslpi, wslpj, uslp, vslp, &
41      &                   pahu, pahv, pgu, pgv, pgui, pgvi, ptb , ptbb, pta, &
42      &                   l_hst, l_ptr, pakz, pah_wslp2, pftu, pftv )
43
44      !!----------------------------------------------------------------------
45      !!                  ***  ROUTINE tra_ldf_iso_tile  ***
46      !!
47      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
48      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
49      !!      add it to the general trend of tracer equation.
50      !!
51      !! ** Method  :   The horizontal component of the lateral diffusive trends
52      !!      is provided by a 2nd order operator rotated along neural or geopo-
53      !!      tential surfaces to which an eddy induced advection can be added
54      !!      It is computed using before fields (forward in time) and isopyc-
55      !!      nal or geopotential slopes computed in routine ldfslp.
56      !!
57      !!      1st part :  masked horizontal derivative of T  ( di[ t ] )
58      !!      ========    with partial cell update if ln_zps=T
59      !!                  with top     cell update if ln_isfcav
60      !!
61      !!      2nd part :  horizontal fluxes of the lateral mixing operator
62      !!      ========   
63      !!         zftu =  pahu e2u*e3u/e1u di[ tb ]
64      !!               - pahu e2u*uslp    dk[ mi(mk(tb)) ]
65      !!         zftv =  pahv e1v*e3v/e2v dj[ tb ]
66      !!               - pahv e2u*vslp    dk[ mj(mk(tb)) ]
67      !!      take the horizontal divergence of the fluxes:
68      !!         difft = 1/(e1e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  }
69      !!      Add this trend to the general trend (ta,sa):
70      !!         ta = ta + difft
71      !!
72      !!      3rd part: vertical trends of the lateral mixing operator
73      !!      ========  (excluding the vertical flux proportional to dk[t] )
74      !!      vertical fluxes associated with the rotated lateral mixing:
75      !!         zftw = - {  mi(mk(pahu)) * e2t*wslpi di[ mi(mk(tb)) ]
76      !!                   + mj(mk(pahv)) * e1t*wslpj dj[ mj(mk(tb)) ]  }
77      !!      take the horizontal divergence of the fluxes:
78      !!         difft = 1/(e1e2t*e3t) dk[ zftw ]
79      !!      Add this trend to the general trend (ta,sa):
80      !!         pta = pta + difft
81      !!
82      !! ** Action :   Update pta arrays with the before rotated diffusion
83      !!----------------------------------------------------------------------
84      TYPE(tile_type)                      , INTENT(in   ) ::   tile
85      LOGICAL                              , INTENT(in   ) ::   ln_traldf_msc, ln_traldf_blp, ln_traldf_lap, ln_zps, ln_isfcav
86      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers passed in
87      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
88      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! 2 . dt
89      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(in   ) ::   e1t, e1u, e1v, e2t, e2u, e2v, e1e2t, e1_e2v, e2_e1u, r1_e1e2t     
90      INTEGER, DIMENSION(jpi,jpj)          , INTENT(in   ) ::   mbku, mbkv, miku, mikv
91      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   umask, vmask, wmask
92      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   e3t_n, e3u_n, e3v_n, e3w_n, wslpi, wslpj, uslp, vslp
93
94      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
95      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
96      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
97      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2)
98      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2)
99      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
100
101      LOGICAL                              , INTENT(in   ) ::   l_ptr      ! flag to compute poleward transport
102      LOGICAL                              , INTENT(in   ) ::   l_hst      ! flag to compute heat transport
103      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(inout) ::   pakz, pah_wslp2 !       
104      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pftu, pftv ! diagnostic outputs
105
106      !
107      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
108      INTEGER  ::  ikt
109      INTEGER  ::  ierr             ! local integer
110      REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars
111      REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      -
112      REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      -
113
114      REAL(wp) ,ALLOCATABLE, DIMENSION(:,:)     ::   zdkt, zdk1t, z2d      ! Make these use global scratch at some point
115      REAL(wp) ,ALLOCATABLE, DIMENSION(:,:,:)   ::   zdit, zdjt, zftu, zftv, ztfw,  zakz, zah_wslp2 
116
117      ALLOCATE(zdkt (tile%jsih:tile%jeih, tile%jsjh:tile%jejh))
118      ALLOCATE(zdk1t(tile%jsih:tile%jeih, tile%jsjh:tile%jejh))
119      ALLOCATE(z2d  (tile%jsih:tile%jeih, tile%jsjh:tile%jejh))
120
121      ALLOCATE(zdit(     tile%jsih:tile%jeih, tile%jsjh:tile%jejh,tile%jskh:tile%jekh))
122      ALLOCATE(zdjt(     tile%jsih:tile%jeih, tile%jsjh:tile%jejh,tile%jskh:tile%jekh))
123      ALLOCATE(zftu(     tile%jsih:tile%jeih, tile%jsjh:tile%jejh,tile%jskh:tile%jekh))
124      ALLOCATE(zftv(     tile%jsih:tile%jeih, tile%jsjh:tile%jejh,tile%jskh:tile%jekh))
125      ALLOCATE(ztfw(     tile%jsih:tile%jeih, tile%jsjh:tile%jejh,tile%jskh:tile%jekh))
126      ALLOCATE(zakz(     tile%jsih:tile%jeih, tile%jsjh:tile%jejh,tile%jskh:tile%jekh))
127      ALLOCATE(zah_wslp2(tile%jsih:tile%jeih, tile%jsjh:tile%jejh,tile%jskh:tile%jekh))
128
129write(0,*)"tile%jsih",tile%jsih
130write(0,*)"tile%jeih",tile%jeih
131write(0,*)"tile%jsjh",tile%jsjh
132write(0,*)"tile%jejh",tile%jejh
133write(0,*)"tile%jskh",tile%jskh
134write(0,*)"tile%jekh",tile%jekh
135
136!$ACC KERNELS
137      !!----------------------------------------------------------------------
138      !
139      !   
140      !
141      !                                            ! set time step size (Euler/Leapfrog)
142      z1_2dt = 1._wp / p2dt   !  note this was dependent on kt and neuler before
143      !
144      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
145      ELSE                    ;   zsign = -1._wp
146      ENDIF
147         
148      !!----------------------------------------------------------------------
149      !!   0 - calculate  ah_wslp2 and akz
150      !!----------------------------------------------------------------------
151      !
152      IF( kpass == 1 ) THEN                  !==  first pass only  ==!
153         !
154         DO jk = tile % jsk_2, tile % jekp1 
155            DO jj = tile % jsj, tile % jej
156               DO ji = tile % jsi, tile % jei   ! vector opt.
157
158                  !
159                  zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
160                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
161                  zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
162                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
163                     !
164                  zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    &
165                     &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku
166                  zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    &
167                     &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv
168                     !
169                  zah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
170                     &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk)
171               END DO
172            END DO
173         END DO
174         !
175         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
176            DO jk = tile % jsk_2, tile % jekp1 
177               DO jj = tile % jsj, tile % jej
178                  DO ji = tile % jsi, tile % jei
179                     zakz(ji,jj,jk) = 0.25_wp * (                                                                     &
180                        &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   &
181                        &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   &
182                        &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   &
183                        &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   )
184                  END DO
185               END DO
186            END DO
187            !
188            IF( ln_traldf_blp ) THEN                ! bilaplacian operator
189               DO jk = tile % jsk_2, tile % jekp1 
190                  DO jj = tile % jsj, tile % jej
191                     DO ji = tile % jsi, tile % jei
192                        zakz(ji,jj,jk) = 16._wp * zah_wslp2(ji,jj,jk)   &
193                           &          * (  zakz(ji,jj,jk) + zah_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               DO jk = tile % jsk_2, tile % jekp1
199                  DO jj = tile % jsj, tile % jej
200                     DO ji = tile % jsi, tile % jei
201                        ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk)
202                        zcoef0 = p2dt * (  zakz(ji,jj,jk) + zah_wslp2(ji,jj,jk) / ze3w_2  )
203                        zakz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
204                     END DO
205                  END DO
206               END DO
207           ENDIF
208           !
209         ELSE                                    ! 33 flux set to zero with zakz=zah_wslp2 ==>> computed in full implicit
210           DO jk = tile % jsk_2, tile % jekp1
211              DO jj = tile % jsj, tile % jej
212                 DO ji = tile % jsi, tile % jei
213                    zakz(ji,jj,jk) = zah_wslp2(ji,jj,jk)
214                 END DO
215               END DO
216            END DO
217         ENDIF
218      ENDIF
219
220! Save zakz and zah_wslp2 (private arrays) in non-overlapping elements of akz and ah_wslp2 arrays (for re-use in trazdf)
221      DO jk = tile % jsk, tile % jek
222         DO jj = tile % jsj, tile % jej
223            DO ji = tile % jsi, tile % jei
224               pakz(ji,jj,jk) = zakz(ji,jj,jk) ; pah_wslp2(ji,jj,jk) = zah_wslp2(ji,jj,jk)
225            END DO
226         END DO
227      END DO
228
229      !
230      !                                                          ! ===========
231      DO jn = 1, kjpt                                            ! tracer loop
232         !                                                       ! ===========
233         !                                               
234         !!----------------------------------------------------------------------
235         !!   I - masked horizontal derivative
236         !!----------------------------------------------------------------------
237
238         ! Horizontal tracer gradient
239         DO jk = tile % jskm1, tile % jekp1
240            DO jj = tile % jsjm1, tile % jej
241               DO ji = tile % jsim1, tile % jei   ! vector opt.
242                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
243                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
244               END DO
245            END DO
246         END DO
247         IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient
248            DO jj = tile % jsjm1, tile % jej              ! bottom correction (partial bottom cell)
249               DO ji = tile % jsim1, tile % jei   ! vector opt.
250                  IF( tile % jskm1 <= mbku(ji,jj) .AND. mbku(ji,jj) <= tile % jekp1 ) zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)         
251                  IF( tile % jskm1 <= mbkv(ji,jj) .AND. mbkv(ji,jj) <= tile % jekp1 ) zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
252               END DO
253            END DO
254            IF( ln_isfcav ) THEN      ! first wet level beneath a cavity
255               DO jj = tile % jsjm1, tile % jej
256                  DO ji = tile % jsim1, tile % jei   ! vector opt.
257                     IF( tile % jskm1 <= miku(ji,jj) .AND. miku(ji,jj) <= tile % jekp1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)         
258                     IF( tile % jskm1 <= mikv(ji,jj) .AND. mikv(ji,jj) <= tile % jekp1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)
259                  END DO
260               END DO
261            ENDIF
262         ENDIF
263         !
264         !!----------------------------------------------------------------------
265         !!   II - horizontal trend  (full)
266         !!----------------------------------------------------------------------
267         !
268         DO jk = tile % jsk, tile % jek                                 ! Horizontal slab
269            !
270            !                             !== Vertical tracer gradient (loops expanded for low-level tiling)
271            DO jj = tile % jsjm1, tile % jejp1
272              DO ji = tile % jsim1, tile % jeip1
273                zdk1t(ji,jj) = ( ptb(ji,jj,jk,jn) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1)     ! level jk+1
274              END DO
275            END DO
276            !
277            IF( jk == 1 ) THEN   
278              DO jj = tile % jsjm1, tile % jejp1
279                DO ji = tile % jsim1, tile % jeip1
280                  zdkt(ji,jj) = zdk1t(ji,jj)    ! surface: zdkt(jk=1)=zdkt(jk=2)
281                END DO
282              END DO
283            ELSE                 
284              DO jj = tile % jsjm1, tile % jejp1
285                DO ji = tile % jsim1, tile % jeip1
286                  zdkt(ji,jj) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
287                END DO
288              END DO
289            ENDIF
290            DO jj = tile % jsjm1, tile % jej            !==  Horizontal fluxes
291               DO ji = tile % jsim1, tile % jei   ! vector opt.
292                  zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
293                  zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
294                  !
295                  zmsku = 1._wp / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   &
296                     &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1._wp )
297                  !
298                  zmskv = 1._wp / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   &
299                     &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1._wp )
300                  !
301                  zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
302                  zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
303                  !
304                  zftu(ji,jj,jk) = (  zabe1 * zdit(ji,jj,jk)   &
305                     &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      &
306                     &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk)
307                  zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   &
308                     &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      &
309                     &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                 
310               END DO
311            END DO
312            !
313            DO jj = tile % jsj , tile % jej          !== horizontal divergence and add to pta
314               DO ji = tile % jsi, tile % jei    ! vector opt.
315                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      &
316                     &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   &
317                     &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
318               END DO
319            END DO
320         END DO                                        !   End of slab 
321
322         !!----------------------------------------------------------------------
323         !!   III - vertical trend (full)
324         !!----------------------------------------------------------------------
325         !
326         ! Vertical fluxes
327         ! ---------------
328         !                          ! Surface and bottom vertical fluxes set to zero
329         IF ( tile % jsk == 1  ) THEN
330            DO jj = tile % jsj, tile % jej
331               DO ji = tile % jsi, tile % jei   ! vector opt.     
332                 ztfw(ji,jj,1) = 0._wp   
333               END DO
334             END DO
335         END IF   
336         IF ( tile % jek == jpkm1 ) THEN
337            DO jj = tile % jsj, tile % jej
338               DO ji = tile % jsi, tile % jei   ! vector opt.     
339                 ztfw(ji,jj,jpk) = 0._wp   
340               END DO
341             END DO
342         END IF   
343
344         DO jk = tile % jsk_2, tile % jekp1           ! interior (2=<jk=<jpk-1)
345            DO jj = tile % jsj, tile % jej
346               DO ji = tile % jsi, tile % jei   ! vector opt.
347                  !
348                  zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
349                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
350                  zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
351                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
352                     !
353                  zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    &
354                     &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku
355                  zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    &
356                     &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv
357                     !
358                  zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked
359                  zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
360                  !
361                  ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      &
362                     &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   &
363                     &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      &
364                     &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  )
365               END DO
366            END DO
367         END DO
368
369         !                                !==  add the vertical 33 flux  ==!
370         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = zah_wslp2 - zakz
371            DO jk = tile % jsk_2, tile % jekp1       
372               DO jj = tile % jsj, tile % jej
373                  DO ji = tile % jsi, tile % jei
374                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   &
375                        &                            * ( zah_wslp2(ji,jj,jk) - zakz(ji,jj,jk) )           &
376                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
377                  END DO
378               END DO
379            END DO
380            !
381         ELSE                                   ! bilaplacian
382            SELECT CASE( kpass )
383            CASE(  1  )                            ! 1st pass : eddy coef = zah_wslp2
384               DO jk = tile % jsk_2, tile % jekp1 
385                  DO jj = tile % jsj, tile % jej
386                     DO ji = tile % jsi, tile % jei
387                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk)                       &
388                           &           + zah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   &
389                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
390                     END DO
391                  END DO
392               END DO
393            CASE(  2  )                         ! 2nd pass : eddy flux = zah_wslp2 and akz applied on ptb  and ptbb gradients, resp.
394               DO jk = tile % jsk_2, tile % jekp1
395                  DO jj = tile % jsj, tile % jej
396                     DO ji = tile % jsi, tile % jei
397                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                       &
398                           &                            * (  zah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
399                           &                               +      zakz(ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
400                     END DO
401                  END DO
402               END DO
403            END SELECT
404         ENDIF
405         !         
406         DO jk = tile % jsk, tile % jek                 !==  Divergence of vertical fluxes added to pta  ==!
407            DO jj = tile % jsj, tile % jej
408               DO ji = tile % jsi, tile % jei   ! vector opt.
409                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   &
410                     &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
411               END DO
412            END DO
413         END DO
414         !
415         IF ( l_ptr .OR. l_hst ) THEN
416            DO jk = tile % jsk, tile % jek                 
417               DO jj = tile % jsj, tile % jej
418                  DO ji = tile % jsi, tile % jei   ! vector opt.
419                     pftu(ji,jj,jk,jn) = zftu(ji,jj,jk) ; pftv(ji,jj,jk,jn) = zftv(ji,jj,jk)
420                  END DO
421               END DO
422            END DO
423         END IF 
424 
425         !                                                        ! ===============
426      END DO                                                      ! end tracer loop
427      !
428!$ACC END KERNELS
429   END SUBROUTINE tra_ldf_iso_tile
430
431   !!==============================================================================
432END MODULE traldf_iso_tile
Note: See TracBrowser for help on using the repository browser.