source: NEMO/branches/UKMO/dev_r12745_HPC-02_Daley_Tiling_trial_public/src/OCE/TRA/traldf_iso.F90 @ 12879

Last change on this file since 12879 was 12879, checked in by hadcv, 9 months ago

Changes following feedback

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