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.
traadv_fct_lf.F90 in NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r13508_HPC-09_loop_fusion/src/OCE/TRA/traadv_fct_lf.F90 @ 13881

Last change on this file since 13881 was 13881, checked in by francesca, 3 years ago

loop fusion v1 - mus and fct advection schemes - ticket #2367
File size: 29.3 KB
Line 
1MODULE traadv_fct_lf
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_fct  ***
4   !! Ocean  tracers:  horizontal & vertical advective trend (2nd/4th order Flux Corrected Transport method)
5   !!==============================================================================
6   !! History :  3.7  !  2015-09  (L. Debreu, G. Madec)  original code (inspired from traadv_tvd.F90)
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!  tra_adv_fct-lf : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme
11   !!                   with sub-time-stepping in the vertical direction - loop fusion version
12   !!  nonosc_lf      : compute monotonic tracer fluxes by a non-oscillatory algorithm - loop fusion version
13   !!----------------------------------------------------------------------
14   USE oce            ! ocean dynamics and active tracers
15   USE dom_oce        ! ocean space and time domain
16   USE trc_oce        ! share passive tracers/Ocean variables
17   USE trd_oce        ! trends: ocean variables
18   USE trdtra         ! tracers trends
19   USE diaptr         ! poleward transport diagnostics
20   USE diaar5         ! AR5 diagnostics
21   USE phycst  , ONLY : rho0_rcp
22   USE zdf_oce , ONLY : ln_zad_Aimp
23   !
24   USE in_out_manager ! I/O manager
25   USE iom            !
26   USE lib_mpp        ! MPP library
27   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
28   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
29   USE traadv_fct
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   tra_adv_fct_lf     ! called by traadv.F90
35
36   LOGICAL  ::   l_trd   ! flag to compute trends
37   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
38   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport
39   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
40
41   !                                        ! tridiag solver associated indices:
42   INTEGER, PARAMETER ::   np_NH   = 0   ! Neumann homogeneous boundary condition
43   INTEGER, PARAMETER ::   np_CEN2 = 1   ! 2nd order centered  boundary condition
44
45   !! * Substitutions
46#  include "do_loop_substitute.h90"
47#  include "domzgr_substitute.h90"
48
49#define search_in_neighbour(out,OP,vec,ji,jj,jk) \
50out = OP(vec(ji,jj,jk),vec(ji-1,jj,jk),vec(ji+1,jj,jk),vec(ji,jj-1,jk),vec(ji,jj+1,jk),vec(ji,jj,MAX(jk-1,1)),vec(ji,jj,jk+1))
51
52#define pos_part_of_flux(ji,jj,jk,out) \
53out = MAX(0.,paa_in(ji-1,jj,jk)) - MIN(0.,paa_in(ji,jj,jk)) \
54+ MAX(0.,pbb_in(ji,jj-1,jk)) - MIN(0.,pbb_in(ji,jj,jk)) \
55+ MAX(0.,pcc_in(ji,jj,jk+1)) - MIN(0.,pcc_in(ji,jj,jk)) 
56
57#define neg_part_of_flux(ji,jj,jk,out) \
58out = MAX( 0.,paa_in(ji,jj,jk) ) - MIN( 0., paa_in(ji-1,jj,jk)) \ 
59+ MAX( 0.,pbb_in(ji,jj,jk) ) - MIN( 0., pbb_in(ji,jj-1,jk)) \
60+ MAX( 0.,pcc_in(ji,jj,jk) ) - MIN( 0., pcc_in(ji,jj,jk+1))
61
62#define beta_terms(bt,betup,betdo,up,pos,do,neg,ji,jj,jk) \
63bt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt ; \
64betup = ( up            - paft(ji,jj,jk) ) / ( pos + zrtrn ) * bt ; \
65betdo = ( paft(ji,jj,jk) - do            ) / ( neg + zrtrn ) * bt
66
67#define monotonic_flux(a,b,c,betup_p1,betdo_p1,vec,vec_in,jk) \
68a = MIN( 1._wp, zbetdo, betup_p1 ) ; \
69b = MIN( 1._wp, zbetup, betdo_p1 ) ; \
70c = ( 0.5_wp  + SIGN( 0.5_wp , vec_in(ji,jj,jk) ) ) ; \
71vec(ji,jj,jk) = vec_in(ji,jj,jk) * ( c * a + ( 1._wp - c) * b )
72
73#define monotonic_flux_k(a,b,c,betup_p1,betdo_p1,vec,vec_in,jk) \
74a = MIN( 1._wp, betdo_p1, zbetup ) ; \
75b = MIN( 1._wp, betup_p1, zbetdo ) ; \
76c = ( 0.5 + SIGN( 0.5_wp , vec_in(ji,jj,jk+1) ) ) ; \
77vec(ji,jj,jk+1) = vec_in(ji,jj,jk+1) * ( c * a + ( 1._wp - c) * b )
78
79   !!----------------------------------------------------------------------
80   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
81   !! $Id: traadv_fct.F90 13660 2020-10-22 10:47:32Z francesca $
82   !! Software governed by the CeCILL license (see ./LICENSE)
83   !!----------------------------------------------------------------------
84CONTAINS
85
86   SUBROUTINE tra_adv_fct_lf( kt, kit000, cdtype, p2dt, pU, pV, pW,       &
87      &                    Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v )
88      !!----------------------------------------------------------------------
89      !!                  ***  ROUTINE tra_adv_fct  ***
90      !!
91      !! **  Purpose :   Compute the now trend due to total advection of tracers
92      !!               and add it to the general trend of tracer equations
93      !!
94      !! **  Method  : - 2nd or 4th FCT scheme on the horizontal direction
95      !!               (choice through the value of kn_fct)
96      !!               - on the vertical the 4th order is a compact scheme
97      !!               - corrected flux (monotonic correction)
98      !!
99      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
100      !!             - send trends to trdtra module for further diagnostics (l_trdtra=T)
101      !!             - poleward advective heat and salt transport (ln_diaptr=T)
102      !!----------------------------------------------------------------------
103      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
104      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
105      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
106      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
107      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
108      INTEGER                                  , INTENT(in   ) ::   kn_fct_h        ! order of the FCT scheme (=2 or 4)
109      INTEGER                                  , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4)
110      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
111      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(inout) ::   pU, pV, pW      ! 3 ocean volume flux components
112      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
113      !
114      INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices 
115      REAL(wp) ::   ztra                                     ! local scalar
116      REAL(wp) ::   zwx, zwx_im1, zfp_ui, zfp_ui_m1, zfp_vj, zfp_vj_m1, zfp_wk, zC2t_u, zC4t_u   !   -      -
117      REAL(wp) ::   zwy, zwy_jm1, zfm_ui, zfm_ui_m1, zfm_vj, zfm_vj_m1, zfm_wk, zC2t_v, zC4t_v   !   -      -
118      REAL(wp) ::   ztu_im1, ztu_ip1, ztv_jm1, ztv_jp1, zltu_ip1, zltv_jp1, zltu, zltv
119      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwx_3d, zwy_3d, zwz, ztw, zltu_3d, zltv_3d, ztu, ztv
120      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdx, ztrdy, ztrdz, zptry
121      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zwinf, zwdia, zwsup
122      LOGICAL  ::   ll_zAimp                                 ! flag to apply adaptive implicit vertical advection
123      !!----------------------------------------------------------------------
124      !
125      IF( kt == kit000 )  THEN
126         IF(lwp) WRITE(numout,*)
127         IF(lwp) WRITE(numout,*) 'tra_adv_fct_lf : FCT advection scheme on ', cdtype
128         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
129      ENDIF
130      !! -- init to 0
131      !! -- init to 0
132      zwi(:,:,:) = 0._wp
133      zwx_3d(:,:,:) = 0._wp
134      zwy_3d(:,:,:) = 0._wp
135      zwz(:,:,:) = 0._wp
136      zwi(:,:,:) = 0._wp
137      !
138      l_trd = .FALSE.            ! set local switches
139      l_hst = .FALSE.
140      l_ptr = .FALSE.
141      ll_zAimp = .FALSE.
142      IF( ( cdtype == 'TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
143      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) )    l_ptr = .TRUE. 
144      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR.  &
145         &                         iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
146      !
147      IF( l_trd .OR. l_hst )  THEN
148         ALLOCATE( ztrdx(jpi,jpj,jpk), ztrdy(jpi,jpj,jpk), ztrdz(jpi,jpj,jpk) )
149         ztrdx(:,:,:) = 0._wp   ;    ztrdy(:,:,:) = 0._wp   ;   ztrdz(:,:,:) = 0._wp
150      ENDIF
151      !
152      IF( l_ptr ) THEN 
153         ALLOCATE( zptry(jpi,jpj,jpk) )
154         zptry(:,:,:) = 0._wp
155      ENDIF
156      !
157      ! If adaptive vertical advection, check if it is needed on this PE at this time
158      IF( ln_zad_Aimp ) THEN
159         IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE.
160      END IF
161      ! If active adaptive vertical advection, build tridiagonal matrix
162      IF( ll_zAimp ) THEN
163         ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk))
164         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
165            zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )   &
166            &                               / e3t(ji,jj,jk,Krhs)
167            zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t(ji,jj,jk,Krhs)
168            zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs)
169         END_3D
170      END IF
171      !
172      DO jn = 1, kjpt            !==  loop over the tracers  ==!
173         !
174         !        !==  upstream advection with initial mass fluxes & intermediate update  ==!
175         !                    !* upstream tracer flux in the i and j direction
176         !                               !* upstream tracer flux in the k direction *!
177         DO_3D( 1, 1, 1, 1, 2, jpkm1 )      ! Interior value ( multiplied by wmask)
178            zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
179            zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
180            zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb) ) * wmask(ji,jj,jk)
181         END_3D
182         IF( ln_linssh ) THEN               ! top ocean value (only in linear free surface as zwz has been w-masked)
183            IF( ln_isfcav ) THEN                        ! top of the ice-shelf cavities and at the ocean surface
184               DO_2D( 1, 1, 1, 1 )
185                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface
186               END_2D
187            ELSE                                        ! no cavities: only at the ocean surface
188               DO_2D( 1, 1, 1, 1 )
189                  zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb)
190               END_2D
191            ENDIF
192         ENDIF
193         !               
194         DO jk = 1, jpkm1
195            DO jj = 1, jpj-1
196               zfp_ui = pU(1,jj,jk) + ABS( pU(1,jj,jk) )
197               zfm_ui = pU(1,jj,jk) - ABS( pU(1,jj,jk) )
198               zwx = 0.5 * ( zfp_ui * pt(1,jj,jk,jn,Kbb) + zfm_ui * pt(2,jj  ,jk,jn,Kbb) )
199               zwx_3d(1,jj,jk) = zwx
200
201               zfp_vj = pV(1,jj,jk) + ABS( pV(1,jj,jk) )
202               zfm_vj = pV(1,jj,jk) - ABS( pV(1,jj,jk) )
203               zwy = 0.5 * ( zfp_vj * pt(1,jj,jk,jn,Kbb) + zfm_vj * pt(1  ,jj+1,jk,jn,Kbb) )
204               zwy_3d(1,jj,jk) = zwy
205            END DO
206            DO ji = 1, jpi-1
207               zfp_ui = pU(ji,1,jk) + ABS( pU(ji,1,jk) )
208               zfm_ui = pU(ji,1,jk) - ABS( pU(ji,1,jk) )
209               zwx = 0.5 * ( zfp_ui * pt(ji,1,jk,jn,Kbb) + zfm_ui * pt(ji+1,1  ,jk,jn,Kbb) )
210               zwx_3d(ji,1,jk) = zwx
211
212               zfp_vj = pV(ji,1,jk) + ABS( pV(ji,1,jk) )
213               zfm_vj = pV(ji,1,jk) - ABS( pV(ji,1,jk) )
214               zwy = 0.5 * ( zfp_vj * pt(ji,1,jk,jn,Kbb) + zfm_vj * pt(ji  ,2,jk,jn,Kbb) )
215               zwy_3d(ji,1,jk) = zwy
216            END DO
217            DO_2D( 1, 1, 1, 1 )
218               zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )
219               zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) )
220               zwx = 0.5 * ( zfp_ui * pt(ji,jj,jk,jn,Kbb) + zfm_ui * pt(ji+1,jj  ,jk,jn,Kbb) )
221               zwx_3d(ji,jj,jk) = zwx
222
223               zfp_ui_m1 = pU(ji-1,jj,jk) + ABS( pU(ji-1,jj,jk) )
224               zfm_ui_m1 = pU(ji-1,jj,jk) - ABS( pU(ji-1,jj,jk) )
225               zwx_im1 = 0.5 * ( zfp_ui_m1 * pt(ji-1,jj,jk,jn,Kbb) + zfm_ui_m1 * pt(ji,jj  ,jk,jn,Kbb) )
226
227               zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) )
228               zfm_vj = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) )
229               zwy = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji  ,jj+1,jk,jn,Kbb) )
230               zwy_3d(ji,jj,jk) = zwy
231
232               zfp_vj_m1 = pV(ji,jj-1,jk) + ABS( pV(ji,jj-1,jk) )
233               zfm_vj_m1 = pV(ji,jj-1,jk) - ABS( pV(ji,jj-1,jk) )
234               zwy_jm1 = 0.5 * ( zfp_vj_m1 * pt(ji,jj-1,jk,jn,Kbb) + zfm_vj_m1 * pt(ji,jj,jk,jn,Kbb) )
235               !                               ! total intermediate advective trends
236               ztra = - (  zwx - zwx_im1 + zwy - zwy_jm1 + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
237               !                               ! update and guess with monotonic sheme
238               pt(ji,jj,jk,jn,Krhs) =                   pt(ji,jj,jk,jn,Krhs) +       ztra   &
239                  &                                  / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk)
240               zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) &
241                  &                                  / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
242            END_2D
243         END DO
244         
245         IF ( ll_zAimp ) THEN
246            CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
247            !
248            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
249            DO_3D( 1, 1, 1, 1, 2, jpkm1 )       ! Interior value ( multiplied by wmask)
250               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
251               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
252               ztw(ji,jj,jk) =  0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk)
253               zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes
254            END_3D
255            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
256               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
257                  &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
258            END_3D
259            !
260         END IF
261         !               
262         IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes)
263            ztrdx(:,:,:) = zwx_3d(:,:,:)   ;   ztrdy(:,:,:) = zwy_3d(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:)
264         END IF
265         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes)
266         IF( l_ptr )   zptry(:,:,:) = zwy_3d(:,:,:) 
267         !
268         !        !==  anti-diffusive flux : high order minus low order  ==!
269         !
270         SELECT CASE( kn_fct_h )    !* horizontal anti-diffusive fluxes
271         !
272         CASE(  2  )                   !- 2nd order centered
273            DO_3D( 2, 1, 2, 1, 1, jpkm1 )
274               zwx_3d(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx_3d(ji,jj,jk)
275               zwy_3d(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy_3d(ji,jj,jk)
276            END_3D
277            !
278         CASE(  4  )                   !- 4th order centered
279            zltu_3d(:,:,jpk) = 0._wp            ! Bottom value : flux set to zero
280            zltv_3d(:,:,jpk) = 0._wp
281            DO jk = 1, jpkm1                 ! Laplacian
282               DO_2D( 1, 0, 1, 0 )                 ! 1st derivative (gradient)
283                  ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
284                  ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
285               END_2D
286               DO_2D( 0, 0, 0, 0 )                 ! 2nd derivative * 1/ 6
287                  zltu_3d(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6
288                  zltv_3d(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6
289               END_2D
290            END DO
291            CALL lbc_lnk_multi( 'traadv_fct', zltu_3d, 'T', 1.0_wp , zltv_3d, 'T', 1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
292!            !
293            DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )
294               zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points
295               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
296               !                                                        ! C4 minus upstream advective fluxes
297               zwx_3d(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu_3d(ji,jj,jk) - zltu_3d(ji+1,jj,jk) ) - zwx_3d(ji,jj,jk)
298               zwy_3d(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv_3d(ji,jj,jk) - zltv_3d(ji,jj+1,jk) ) - zwy_3d(ji,jj,jk)
299            END_3D
300            !
301            CALL lbc_lnk_multi( 'traadv_fct', zwx_3d, 'U', -1.0_wp , zwy_3d, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
302         CASE(  41 )                   !- 4th order centered       ==>>   !!gm coding attempt   need to be tested
303            DO_3D( 0, 0, 0, 0, 1, jpkm1 )    ! Horizontal advective fluxes
304               ztu_im1 = ( pt(ji  ,jj  ,jk,jn,Kmm) - pt(ji-1,jj,jk,jn,Kmm) ) * umask(ji-1,jj,jk)
305               ztu_ip1 = ( pt(ji+2,jj  ,jk,jn,Kmm) - pt(ji+1,jj,jk,jn,Kmm) ) * umask(ji+1,jj,jk)
306
307               ztv_jm1 = ( pt(ji,jj  ,jk,jn,Kmm) - pt(ji,jj-1,jk,jn,Kmm) ) * vmask(ji,jj-1,jk)
308               ztv_jp1 = ( pt(ji,jj+2,jk,jn,Kmm) - pt(ji,jj+1,jk,jn,Kmm) ) * vmask(ji,jj+1,jk)
309
310               zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points (x2)
311               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
312               !                                                  ! C4 interpolation of T at u- & v-points (x2)
313               zC4t_u =  zC2t_u + r1_6 * ( ztu_im1 - ztu_ip1 )
314               zC4t_v =  zC2t_v + r1_6 * ( ztv_jm1 - ztv_jp1 )
315               !                                                  ! C4 minus upstream advective fluxes
316               zwx_3d(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx_3d(ji,jj,jk)
317               zwy_3d(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy_3d(ji,jj,jk)
318            END_3D
319            CALL lbc_lnk_multi( 'traadv_fct', zwx_3d, 'U', -1.0_wp , zwy_3d, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
320            !
321         END SELECT
322         !                     
323         SELECT CASE( kn_fct_v )    !* vertical anti-diffusive fluxes (w-masked interior values)
324         !
325         CASE(  2  )                   !- 2nd order centered
326            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
327               zwz(ji,jj,jk) =  (  pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
328                  &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk)
329            END_3D
330            !
331         CASE(  4  )                   !- 4th order COMPACT
332            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )   ! zwt = COMPACT interpolation of T at w-point
333            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
334               zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
335            END_3D
336            !
337         END SELECT
338         IF( ln_linssh ) THEN    ! top ocean value: high order = upstream  ==>>  zwz=0
339            zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked
340         ENDIF
341         !         
342         CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp)
343         !
344         IF ( ll_zAimp ) THEN
345            DO_3D( 1, 1, 1, 1, 1, jpkm1 )    !* trend and after field with monotonic scheme
346               !                                                ! total intermediate advective trends
347               ztra = - (  zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj  ,jk  )   &
348                  &      + zwy_3d(ji,jj,jk) - zwy_3d(ji  ,jj-1,jk  )   &
349                  &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
350               ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
351            END_3D
352            !
353            CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
354            !
355            DO_3D( 1, 1, 1, 1, 2, jpkm1 )       ! Interior value ( multiplied by wmask)
356               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
357               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
358               zwz(ji,jj,jk) = zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztw(ji,jj,jk) + zfm_wk * ztw(ji,jj,jk-1) ) * wmask(ji,jj,jk)
359            END_3D
360         END IF
361         !
362         !        !==  monotonicity algorithm  ==!
363         !
364#if defined key_agrif
365         CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx_3d, zwy_3d, zwz, zwi, p2dt )
366#else
367         CALL nonosc_lf( Kmm, pt(:,:,:,jn,Kbb), zwx_3d, zwy_3d, zwz, zwi, p2dt )
368#endif
369         !
370         !        !==  final trend with corrected fluxes  ==!
371         !
372         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
373            ztra = - (  zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj  ,jk  )   &
374               &      + zwy_3d(ji,jj,jk) - zwy_3d(ji  ,jj-1,jk  )   &
375               &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
376            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm)
377            zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
378         END_3D
379         !
380         IF ( ll_zAimp ) THEN
381            !
382            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp
383            DO_3D( 0, 0, 0, 0, 2, jpkm1 )      ! Interior value ( multiplied by wmask)
384               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
385               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
386               ztw(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk)
387               zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic
388            END_3D
389            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
390               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
391                  &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
392            END_3D
393         END IF         
394         ! NOT TESTED - NEED l_trd OR l_hst TRUE
395         IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport
396            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx_3d(:,:,:)  ! <<< add anti-diffusive fluxes
397            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy_3d(:,:,:)  !     to upstream fluxes
398            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  !
399            !
400            IF( l_trd ) THEN              ! trend diagnostics
401               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx, pU, pt(:,:,:,jn,Kmm) )
402               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy, pV, pt(:,:,:,jn,Kmm) )
403               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz, pW, pt(:,:,:,jn,Kmm) )
404            ENDIF
405            !                             ! heat/salt transport
406            IF( l_hst )   CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
407            !
408         ENDIF
409         ! NOT TESTED - NEED l_ptr TRUE
410         IF( l_ptr ) THEN              ! "Poleward" transports
411            zptry(:,:,:) = zptry(:,:,:) + zwy_3d(:,:,:)  ! <<< add anti-diffusive fluxes
412            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
413         ENDIF
414         !
415      END DO                     ! end of tracer loop
416      !
417      IF ( ll_zAimp ) THEN
418         DEALLOCATE( zwdia, zwinf, zwsup )
419      ENDIF
420      IF( l_trd .OR. l_hst ) THEN
421         DEALLOCATE( ztrdx, ztrdy, ztrdz )
422      ENDIF
423      IF( l_ptr ) THEN
424         DEALLOCATE( zptry )
425      ENDIF
426      !
427   END SUBROUTINE tra_adv_fct_lf
428
429   SUBROUTINE nonosc_lf( Kmm, pbef, paa, pbb, pcc, paft, p2dt )
430      !!---------------------------------------------------------------------
431      !!                    ***  ROUTINE nonosc  ***
432      !!     
433      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
434      !!       scheme and the before field by a nonoscillatory algorithm
435      !!
436      !! **  Method  :   ... ???
437      !!       warning : pbef and paft must be masked, but the boundaries
438      !!       conditions on the fluxes are not necessary zalezak (1979)
439      !!       drange (1995) multi-dimensional forward-in-time and upstream-
440      !!       in-space based differencing for fluid
441      !!----------------------------------------------------------------------
442      INTEGER                          , INTENT(in   ) ::   Kmm             ! time level index
443      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! tracer time-step
444      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
445      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb, pcc   ! monotonic fluxes in the 3 directions
446      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   paa_in, pbb_in, pcc_in   ! monotonic fluxes in the 3 directions
447      !
448      REAL(dp), DIMENSION (jpi,jpj,jpk) :: zbup, zbdo
449      !
450      INTEGER  ::   ji, jj, jk   ! dummy loop indices
451      INTEGER  ::   ikm1         ! local integer
452      REAL(dp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars
453      REAL(dp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
454      REAL(dp) ::   zbetup, zbetdo
455      REAL(dp) ::   zbt_ip1, zpos_ip1, zneg_ip1, zup_ip1, zdo_ip1, zbetup_ip1, zbetdo_ip1 
456      REAL(dp) ::   zbt_jp1, zpos_jp1, zneg_jp1, zup_jp1, zdo_jp1, zbetup_jp1, zbetdo_jp1
457      REAL(dp) ::   zbt_kp1, zpos_kp1, zneg_kp1, zup_kp1, zdo_kp1, zbetup_kp1, zbetdo_kp1
458      !!----------------------------------------------------------------------
459      !
460      zbig  = 1.e+40_dp
461      zrtrn = 1.e-15_dp
462
463      paa_in(:,:,:) = paa(:,:,:)
464      pbb_in(:,:,:) = pbb(:,:,:)
465      pcc_in(:,:,:) = pcc(:,:,:)
466      ! Search local extrema
467      ! --------------------
468      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
469      zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ),   &
470         &        paft * tmask - zbig * ( 1._wp - tmask )  )
471      zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ),   &
472         &        paft * tmask + zbig * ( 1._wp - tmask )  )
473
474      DO_3D( 1, 0, 1, 0, 1, jpk-2 )
475
476        ! search maximum in neighbourhood
477        search_in_neighbour(zup,MAX,zbup,ji,jj,jk)
478        search_in_neighbour(zup_ip1,MAX,zbup,ji+1,jj,jk)
479        search_in_neighbour(zup_jp1,MAX,zbup,ji,jj+1,jk)
480        search_in_neighbour(zup_kp1,MAX,zbup,ji,jj,jk+1)
481
482        ! search minimum in neighbourhood
483        search_in_neighbour(zdo,MIN,zbdo,ji,jj,jk)
484        search_in_neighbour(zdo_ip1,MIN,zbdo,ji+1,jj,jk)
485        search_in_neighbour(zdo_jp1,MIN,zbdo,ji,jj+1,jk)
486        search_in_neighbour(zdo_kp1,MIN,zbdo,ji,jj,jk+1)
487
488        ! positive part of the flux
489        pos_part_of_flux(ji,jj,jk,zpos)
490        pos_part_of_flux(ji+1,jj,jk,zpos_ip1)
491        pos_part_of_flux(ji,jj+1,jk,zpos_jp1)
492        pos_part_of_flux(ji,jj,jk+1,zpos_kp1)
493
494        ! negative part of the flux
495        neg_part_of_flux(ji,jj,jk,zneg)
496        neg_part_of_flux(ji+1,jj,jk,zneg_ip1)
497        neg_part_of_flux(ji,jj+1,jk,zneg_jp1)
498        neg_part_of_flux(ji,jj,jk+1,zneg_kp1)
499
500        ! up & down beta terms
501        beta_terms(zbt,zbetup,zbetdo,zup,zpos,zdo,zneg,ji,jj,jk)
502        beta_terms(zbt_ip1,zbetup_ip1,zbetdo_ip1,zup_ip1,zpos_ip1,zdo_ip1,zneg_ip1,ji+1,jj,jk)
503        beta_terms(zbt_jp1,zbetup_jp1,zbetdo_jp1,zup_jp1,zpos_jp1,zdo_jp1,zneg_jp1,ji,jj+1,jk)
504        beta_terms(zbt_kp1,zbetup_kp1,zbetdo_kp1,zup_kp1,zpos_kp1,zdo_kp1,zneg_kp1,ji,jj,jk+1)
505
506        ! 3. monotonic flux in the i & j (paa & pbb)
507        ! ----------------------------------------
508        monotonic_flux(zau,zbu,zcu,zbetup_ip1,zbetdo_ip1,paa,paa_in,jk)
509        monotonic_flux(zav,zbv,zcv,zbetup_jp1,zbetdo_jp1,pbb,pbb_in,jk)
510
511        ! monotonic flux in the k direction, i.e. pcc
512        ! -------------------------------------------
513        monotonic_flux_k(za,zb,zc,zbetup_kp1,zbetdo_kp1,pcc,pcc_in,jk)
514      END_3D
515      !
516      DO_2D( 1, 0, 1, 0 )
517
518        ! search maximum in neighbourhood
519        search_in_neighbour(zup,MAX,zbup,ji,jj,jpk-1)
520        search_in_neighbour(zup_ip1,MAX,zbup,ji+1,jj,jpk-1)
521        search_in_neighbour(zup_jp1,MAX,zbup,ji,jj+1,jpk-1)
522
523        ! search minimum in neighbourhood
524        search_in_neighbour(zdo,MIN,zbdo,ji,jj,jk)
525        search_in_neighbour(zdo_ip1,MIN,zbdo,ji+1,jj,jpk-1)
526        search_in_neighbour(zdo_jp1,MIN,zbdo,ji,jj+1,jpk-1)
527
528        ! positive part of the flux
529        pos_part_of_flux(ji,jj,jpk-1,zpos)
530        pos_part_of_flux(ji+1,jj,jpk-1,zpos_ip1)
531        pos_part_of_flux(ji,jj+1,jpk-1,zpos_jp1)
532
533        ! negative part of the flux
534        neg_part_of_flux(ji,jj,jpk-1,zneg)
535        neg_part_of_flux(ji+1,jj,jpk-1,zneg_ip1)
536        neg_part_of_flux(ji,jj+1,jpk-1,zneg_jp1)
537
538        ! up & down beta terms
539        beta_terms(zbt,zbetup,zbetdo,zup,zpos,zdo,zneg,ji,jj,jpk-1)
540        beta_terms(zbt_ip1,zbetup_ip1,zbetdo_ip1,zup_ip1,zpos_ip1,zdo_ip1,zneg_ip1,ji+1,jj,jpk-1)
541        beta_terms(zbt_jp1,zbetup_jp1,zbetdo_jp1,zup_jp1,zpos_jp1,zdo_jp1,zneg_jp1,ji,jj+1,jpk-1)
542
543        zbetup_kp1 = 0._dp 
544        zbetdo_kp1 = 0._dp
545
546        ! 3. monotonic flux in the i & j direction (paa & pbb)
547        ! ----------------------------------------
548        monotonic_flux(zau,zbu,zcu,zbetup_ip1,zbetdo_ip1,paa,paa_in,jpk-1)
549        monotonic_flux(zav,zbv,zcv,zbetup_jp1,zbetdo_jp1,pbb,pbb_in,jpk-1)
550
551        ! monotonic flux in the k direction, i.e. pcc
552        ! -------------------------------------------
553        monotonic_flux_k(za,zb,zc,zbetup_kp1,zbetdo_kp1,pcc,pcc_in,jpk-1)
554      END_2D
555   END SUBROUTINE nonosc_lf
556
557END MODULE traadv_fct_lf
Note: See TracBrowser for help on using the repository browser.