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/trunk/src/OCE/TRA – NEMO

source: NEMO/trunk/src/OCE/TRA/traadv_fct_lf.F90 @ 14072

Last change on this file since 14072 was 14072, checked in by laurent, 3 years ago

Merging branch "2020/dev_r13648_ASINTER-04_laurent_bulk_ice", ticket #2369

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