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 NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

source: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_iso.F90 @ 14680

Last change on this file since 14680 was 14680, checked in by hadcv, 3 years ago

#2600: Merge in dev_r14393_HPC-03_Mele_Comm_Cleanup [14667]
  • Property svn:keywords set to Id
File size: 26.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 domutl, ONLY : is_tile
22   USE trc_oce        ! share passive tracers/Ocean variables
23   USE zdf_oce        ! ocean vertical physics
24   USE ldftra         ! lateral diffusion: tracer eddy coefficients
25   USE ldfslp         ! iso-neutral slopes
26   USE diaptr         ! poleward transport diagnostics
27   USE diaar5         ! AR5 diagnostics
28   !
29   USE in_out_manager ! I/O manager
30   USE iom            ! I/O library
31   USE phycst         ! physical constants
32   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   tra_ldf_iso   ! routine called by step.F90
38
39   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
40   LOGICAL  ::   l_hst   ! flag to compute heat transport
41
42   !! * Substitutions
43#  include "do_loop_substitute.h90"
44#  include "domzgr_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
47   !! $Id$
48   !! Software governed by the CeCILL license (see ./LICENSE)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE tra_ldf_iso( kt, Kmm, kit000, cdtype, pahu, pahv,             &
53      &                                             pgu , pgv , pgui, pgvi, &
54      &                                             pt, pt2, pt_rhs, kjpt, kpass )
55      !!
56      INTEGER                     , INTENT(in   ) ::   kt         ! ocean time-step index
57      INTEGER                     , INTENT(in   ) ::   kit000     ! first time step index
58      CHARACTER(len=3)            , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
59      INTEGER                     , INTENT(in   ) ::   kjpt       ! number of tracers
60      INTEGER                     , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
61      INTEGER                     , INTENT(in   ) ::   Kmm        ! ocean time level index
62      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
63      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
64      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
65      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pt         ! tracer (kpass=1) or laplacian of tracer (kpass=2)
66      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pt2        ! tracer (only used in kpass=2)
67      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pt_rhs     ! tracer trend
68      !!
69      CALL tra_ldf_iso_t( kt, Kmm, kit000, cdtype, pahu, pahv, is_tile(pahu),                             &
70         &                                         pgu , pgv , is_tile(pgu) , pgui, pgvi, is_tile(pgui),  &
71         &                                         pt, is_tile(pt), pt2, is_tile(pt2), pt_rhs, is_tile(pt_rhs), kjpt, kpass )
72   END SUBROUTINE tra_ldf_iso
73
74
75  SUBROUTINE tra_ldf_iso_t( kt, Kmm, kit000, cdtype, pahu, pahv, ktah,                    &
76      &                                              pgu , pgv , ktg , pgui, pgvi, ktgi,  &
77      &                                              pt, ktt, pt2, ktt2, pt_rhs, ktt_rhs, kjpt, kpass )
78      !!----------------------------------------------------------------------
79      !!                  ***  ROUTINE tra_ldf_iso  ***
80      !!
81      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
82      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
83      !!      add it to the general trend of tracer equation.
84      !!
85      !! ** Method  :   The horizontal component of the lateral diffusive trends
86      !!      is provided by a 2nd order operator rotated along neural or geopo-
87      !!      tential surfaces to which an eddy induced advection can be added
88      !!      It is computed using before fields (forward in time) and isopyc-
89      !!      nal or geopotential slopes computed in routine ldfslp.
90      !!
91      !!      1st part :  masked horizontal derivative of T  ( di[ t ] )
92      !!      ========    with partial cell update if ln_zps=T
93      !!                  with top     cell update if ln_isfcav
94      !!
95      !!      2nd part :  horizontal fluxes of the lateral mixing operator
96      !!      ========
97      !!         zftu =  pahu e2u*e3u/e1u di[ tb ]
98      !!               - pahu e2u*uslp    dk[ mi(mk(tb)) ]
99      !!         zftv =  pahv e1v*e3v/e2v dj[ tb ]
100      !!               - pahv e2u*vslp    dk[ mj(mk(tb)) ]
101      !!      take the horizontal divergence of the fluxes:
102      !!         difft = 1/(e1e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  }
103      !!      Add this trend to the general trend (ta,sa):
104      !!         ta = ta + difft
105      !!
106      !!      3rd part: vertical trends of the lateral mixing operator
107      !!      ========  (excluding the vertical flux proportional to dk[t] )
108      !!      vertical fluxes associated with the rotated lateral mixing:
109      !!         zftw = - {  mi(mk(pahu)) * e2t*wslpi di[ mi(mk(tb)) ]
110      !!                   + mj(mk(pahv)) * e1t*wslpj dj[ mj(mk(tb)) ]  }
111      !!      take the horizontal divergence of the fluxes:
112      !!         difft = 1/(e1e2t*e3t) dk[ zftw ]
113      !!      Add this trend to the general trend (ta,sa):
114      !!         pt_rhs = pt_rhs + difft
115      !!
116      !! ** Action :   Update pt_rhs arrays with the before rotated diffusion
117      !!----------------------------------------------------------------------
118      INTEGER                                   , INTENT(in   ) ::   kt         ! ocean time-step index
119      INTEGER                                   , INTENT(in   ) ::   kit000     ! first time step index
120      CHARACTER(len=3)                          , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
121      INTEGER                                   , INTENT(in   ) ::   kjpt       ! number of tracers
122      INTEGER                                   , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
123      INTEGER                                   , INTENT(in   ) ::   Kmm        ! ocean time level index
124      INTEGER                                   , INTENT(in   ) ::   ktah, ktg, ktgi, ktt, ktt2, ktt_rhs
125      REAL(wp), DIMENSION(A2D_T(ktah)   ,JPK)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
126      REAL(wp), DIMENSION(A2D_T(ktg)        ,KJPT), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
127      REAL(wp), DIMENSION(A2D_T(ktgi)       ,KJPT), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
128      REAL(wp), DIMENSION(A2D_T(ktt)    ,JPK,KJPT), INTENT(in   ) ::   pt         ! tracer (kpass=1) or laplacian of tracer (kpass=2)
129      REAL(wp), DIMENSION(A2D_T(ktt2)   ,JPK,KJPT), INTENT(in   ) ::   pt2        ! tracer (only used in kpass=2)
130      REAL(wp), DIMENSION(A2D_T(ktt_rhs),JPK,KJPT), INTENT(inout) ::   pt_rhs     ! tracer trend
131      !
132      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
133      INTEGER  ::  ikt
134      INTEGER  ::  ierr, iij        ! local integer
135      REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars
136      REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      -
137      REAL(wp) ::  zcoef0, ze3w_2, zsign                 !   -      -
138      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zdkt, zdk1t, z2d
139      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zdit, zdjt, zftu, zftv, ztfw
140      !!----------------------------------------------------------------------
141      !
142      IF( kpass == 1 .AND. kt == kit000 )  THEN
143         IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
144            IF(lwp) WRITE(numout,*)
145            IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype
146            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
147         ENDIF
148         !
149         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpk )
150         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
151            akz     (ji,jj,jk) = 0._wp
152            ah_wslp2(ji,jj,jk) = 0._wp
153         END_3D
154      ENDIF
155      !
156      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                           ! Do only on the first tile
157         l_hst = .FALSE.
158         l_ptr = .FALSE.
159         IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )     l_ptr = .TRUE.
160         IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
161            &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
162      ENDIF
163      !
164      ! Define pt_rhs halo points for multi-point haloes in bilaplacian case
165      IF( nldf_tra == np_blp_i .AND. kpass == 1 ) THEN ; iij = nn_hls
166      ELSE                                             ; iij = 1
167      ENDIF
168
169      !
170      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
171      ELSE                    ;   zsign = -1._wp
172      ENDIF
173
174      !!----------------------------------------------------------------------
175      !!   0 - calculate  ah_wslp2 and akz
176      !!----------------------------------------------------------------------
177      !
178      IF( kpass == 1 ) THEN                  !==  first pass only  ==!
179         !
180         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
181         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
182            !
183            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
184               &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
185            zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
186               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
187            ! round brackets added to fix the order of floating point operations
188            ! needed to ensure halo 1 - halo 2 compatibility
189            zahu_w = ( (  pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)                    &
190               &       )                                                           & ! bracket for halo 1 - halo 2 compatibility
191               &       + ( pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)                   &
192               &         ) ) * zmsku                                                 ! bracket for halo 1 - halo 2 compatibility
193            zahv_w = ( (  pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)                    &
194               &       )                                                           & ! bracket for halo 1 - halo 2 compatibility
195               &       + ( pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)                   &
196               &         ) ) * zmskv                                                 ! bracket for halo 1 - halo 2 compatibility
197               !
198            ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
199               &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk)
200         END_3D
201         !
202         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
203            ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
204            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
205               ! round brackets added to fix the order of floating point operations
206               ! needed to ensure halo 1 - halo 2 compatibility
207               akz(ji,jj,jk) = 0.25_wp * (                                                                     &
208                  &            ( ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   &
209                  &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   &
210                  &            )                                                                               & ! bracket for halo 1 - halo 2 compatibility
211                  &            + ( ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) ) &
212                  &              + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) &
213                  &              ) )                                                                             ! bracket for halo 1 - halo 2 compatibility
214            END_3D
215            !
216            IF( ln_traldf_blp ) THEN                ! bilaplacian operator
217               ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
218               DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
219                  akz(ji,jj,jk) = 16._wp   &
220                     &   * ah_wslp2   (ji,jj,jk)   &
221                     &   * (  akz     (ji,jj,jk)   &
222                     &      + ah_wslp2(ji,jj,jk)   &
223                     &        / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
224               END_3D
225            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
226               ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
227               DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
228                  ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
229                  zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
230                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_Dt
231               END_3D
232           ENDIF
233           !
234         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
235            ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpk )
236            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
237               akz(ji,jj,jk) = ah_wslp2(ji,jj,jk)
238            END_3D
239         ENDIF
240      ENDIF
241      !
242      !                                                          ! ===========
243      DO jn = 1, kjpt                                            ! tracer loop
244         !                                                       ! ===========
245         !
246         !!----------------------------------------------------------------------
247         !!   I - masked horizontal derivative
248         !!----------------------------------------------------------------------
249         zdit(:,:,:) = 0._wp
250         zdjt(:,:,:) = 0._wp
251
252         ! Horizontal tracer gradient
253         ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 1, jpkm1 )
254         DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 )
255            zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
256            zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
257         END_3D
258         IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient
259            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )           ! bottom correction (partial bottom cell)
260            DO_2D( iij, iij-1, iij, iij-1 )            ! bottom correction (partial bottom cell)
261               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
262               zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
263            END_2D
264            IF( ln_isfcav ) THEN      ! first wet level beneath a cavity
265               ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )
266               DO_2D( iij, iij-1, iij, iij-1 )
267                  IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)
268                  IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)
269               END_2D
270            ENDIF
271         ENDIF
272         !
273         !!----------------------------------------------------------------------
274         !!   II - horizontal trend  (full)
275         !!----------------------------------------------------------------------
276         !
277         DO jk = 1, jpkm1                                 ! Horizontal slab
278            !
279            ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
280            DO_2D( iij, iij, iij, iij )
281               !                             !== Vertical tracer gradient
282               zdk1t(ji,jj) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1)     ! level jk+1
283               !
284               IF( jk == 1 ) THEN   ;   zdkt(ji,jj) = zdk1t(ji,jj)                            ! surface: zdkt(jk=1)=zdkt(jk=2)
285               ELSE                 ;   zdkt(ji,jj) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
286               ENDIF
287            END_2D
288            !
289            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )           !==  Horizontal fluxes
290            DO_2D( iij, iij-1, iij, iij-1 )           !==  Horizontal fluxes
291               zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)
292               zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
293               !
294               zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   &
295                  &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. )
296               !
297               zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   &
298                  &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. )
299               !
300               zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
301               zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
302               !
303               ! round brackets added to fix the order of floating point operations
304               ! needed to ensure halo 1 - halo 2 compatibility
305               zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)                       &
306                  &               + zcof1 * ( ( zdkt (ji+1,jj) + zdk1t(ji,jj)    &
307                  &                           )                                  & ! bracket for halo 1 - halo 2 compatibility
308                  &                         + ( zdk1t(ji+1,jj) + zdkt (ji,jj)    &
309                  &                           )                                  & ! bracket for halo 1 - halo 2 compatibility
310                  &                         ) ) * umask(ji,jj,jk)
311               zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)                        &
312                  &              + zcof2 * ( ( zdkt (ji,jj+1) + zdk1t(ji,jj)     &
313                  &                           )                                  & ! bracket for halo 1 - halo 2 compatibility
314                  &                         + ( zdk1t(ji,jj+1) + zdkt (ji,jj)    &
315                  &                           )                                  & ! bracket for halo 1 - halo 2 compatibility
316                  &                         ) ) * vmask(ji,jj,jk)
317            END_2D
318            !
319            ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )           !== horizontal divergence and add to pta
320            DO_2D( iij-1, iij-1, iij-1, iij-1 )           !== horizontal divergence and add to pta
321               ! round brackets added to fix the order of floating point operations
322               ! needed to ensure halo 1 - halo 2 compatibility
323               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)                         &
324                  &       + zsign * ( ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk)        &
325                  &                   )                                          & ! bracket for halo 1 - halo 2 compatibility
326                  &                 + ( zftv(ji,jj,jk) - zftv(ji,jj-1,jk)        &
327                  &                   )                                          & ! bracket for halo 1 - halo 2 compatibility
328                  &                 ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
329            END_2D
330         END DO                                        !   End of slab
331
332         !!----------------------------------------------------------------------
333         !!   III - vertical trend (full)
334         !!----------------------------------------------------------------------
335         !
336         ! Vertical fluxes
337         ! ---------------
338         !                          ! Surface and bottom vertical fluxes set to zero
339         ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp
340
341         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )    ! interior (2=<jk=<jpk-1)
342         DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 )    ! interior (2=<jk=<jpk-1)
343            !
344            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
345               &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
346            zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
347               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
348               !
349            zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    &
350               &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku
351            zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    &
352               &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv
353               !
354            zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked
355            zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
356            !
357            ! round brackets added to fix the order of floating point operations
358            ! needed to ensure halo 1 - halo 2 compatibility
359            ztfw(ji,jj,jk) = zcoef3 * ( ( zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)    &
360                  &                     )                                              & ! bracket for halo 1 - halo 2 compatibility
361                  &                   + ( zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)    &
362                  &                     )                                              & ! bracket for halo 1 - halo 2 compatibility
363                  &                   )                                                &
364                  &        + zcoef4 * ( ( zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)    &
365                  &                     )                                              & ! bracket for halo 1 - halo 2 compatibility
366                  &                   + ( zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)    &
367                  &                     )                                              & ! bracket for halo 1 - halo 2 compatibility
368                  &                   )
369         END_3D
370         !                                !==  add the vertical 33 flux  ==!
371         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
372            ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
373            DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 )
374               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)   &
375                  &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )               &
376                  &                            * (  pt(ji,jj,jk-1,jn) -  pt(ji,jj,jk,jn) )
377            END_3D
378            !
379         ELSE                                   ! bilaplacian
380            SELECT CASE( kpass )
381            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
382               ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )
383               DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 )
384                  ztfw(ji,jj,jk) =   &
385                     &  ztfw(ji,jj,jk) + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   &
386                     &           * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
387               END_3D
388            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
389               DO_3D( 0, 0, 0, 0, 2, jpkm1 )
390                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)                  &
391                     &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
392                     &                            +         akz(ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   )
393               END_3D
394            END SELECT
395         ENDIF
396         !
397         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !==  Divergence of vertical fluxes added to pta  ==!
398         DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 )    !==  Divergence of vertical fluxes added to pta  ==!
399            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * r1_e1e2t(ji,jj)   &
400               &                                             / e3t(ji,jj,jk,Kmm)
401         END_3D
402         !
403         IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==!
404             ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==!
405            !
406            !                             ! "Poleward" diffusive heat or salt transports (T-S case only)
407               ! note sign is reversed to give down-gradient diffusive transports )
408            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -zftv(:,:,:)  )
409            !                          ! Diffusive heat transports
410            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', -zftu(:,:,:), -zftv(:,:,:) )
411            !
412         ENDIF                                                    !== end pass selection  ==!
413         !
414         !                                                        ! ===============
415      END DO                                                      ! end tracer loop
416      !
417   END SUBROUTINE tra_ldf_iso_t
418
419   !!==============================================================================
420END MODULE traldf_iso
Note: See TracBrowser for help on using the repository browser.