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.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traadv_fct.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.
  • Property svn:keywords set to Id
File size: 32.5 KB
Line 
1MODULE traadv_fct
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    : 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
12   !!  nonosc         : compute monotonic tracer fluxes by a non-oscillatory algorithm
13   !!  interp_4th_cpt : 4th order compact scheme for the vertical component of the advection
14   !!----------------------------------------------------------------------
15   USE oce            ! ocean dynamics and active tracers
16   USE dom_oce        ! ocean space and time domain
17   USE trc_oce        ! share passive tracers/Ocean variables
18   USE trd_oce        ! trends: ocean variables
19   USE trdtra         ! tracers trends
20   USE diaptr         ! poleward transport diagnostics
21   USE diaar5         ! AR5 diagnostics
22   USE phycst  , ONLY : rau0_rcp
23   USE zdf_oce , ONLY : ln_zad_Aimp
24   !
25   USE in_out_manager ! I/O manager
26   USE iom            !
27   USE lib_mpp        ! MPP library
28   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
29   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   tra_adv_fct        ! called by traadv.F90
35   PUBLIC   interp_4th_cpt     ! called by traadv_cen.F90
36
37   LOGICAL  ::   l_trd   ! flag to compute trends
38   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
39   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport
40   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
41
42   !                                        ! tridiag solver associated indices:
43   INTEGER, PARAMETER ::   np_NH   = 0   ! Neumann homogeneous boundary condition
44   INTEGER, PARAMETER ::   np_CEN2 = 1   ! 2nd order centered  boundary condition
45
46   !! * Substitutions
47#  include "vectopt_loop_substitute.h90"
48#  include "do_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
51   !! $Id$
52   !! Software governed by the CeCILL license (see ./LICENSE)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pU, pV, pW,       &
57      &                    Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v )
58      !!----------------------------------------------------------------------
59      !!                  ***  ROUTINE tra_adv_fct  ***
60      !!
61      !! **  Purpose :   Compute the now trend due to total advection of tracers
62      !!               and add it to the general trend of tracer equations
63      !!
64      !! **  Method  : - 2nd or 4th FCT scheme on the horizontal direction
65      !!               (choice through the value of kn_fct)
66      !!               - on the vertical the 4th order is a compact scheme
67      !!               - corrected flux (monotonic correction)
68      !!
69      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
70      !!             - send trends to trdtra module for further diagnostics (l_trdtra=T)
71      !!             - poleward advective heat and salt transport (ln_diaptr=T)
72      !!----------------------------------------------------------------------
73      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
74      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
75      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
76      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
77      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
78      INTEGER                                  , INTENT(in   ) ::   kn_fct_h        ! order of the FCT scheme (=2 or 4)
79      INTEGER                                  , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4)
80      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
81      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components
82      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
83      !
84      INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices 
85      REAL(wp) ::   ztra                                     ! local scalar
86      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u   !   -      -
87      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v   !   -      -
88      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw
89      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdx, ztrdy, ztrdz, zptry
90      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zwinf, zwdia, zwsup
91      LOGICAL  ::   ll_zAimp                                 ! flag to apply adaptive implicit vertical advection
92      !!----------------------------------------------------------------------
93      !
94      IF( kt == kit000 )  THEN
95         IF(lwp) WRITE(numout,*)
96         IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype
97         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
98      ENDIF
99      !
100      l_trd = .FALSE.            ! set local switches
101      l_hst = .FALSE.
102      l_ptr = .FALSE.
103      ll_zAimp = .FALSE.
104      IF( ( cdtype == 'TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
105      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) )    l_ptr = .TRUE. 
106      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR.  &
107         &                         iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
108      !
109      IF( l_trd .OR. l_hst )  THEN
110         ALLOCATE( ztrdx(jpi,jpj,jpk), ztrdy(jpi,jpj,jpk), ztrdz(jpi,jpj,jpk) )
111         ztrdx(:,:,:) = 0._wp   ;    ztrdy(:,:,:) = 0._wp   ;   ztrdz(:,:,:) = 0._wp
112      ENDIF
113      !
114      IF( l_ptr ) THEN 
115         ALLOCATE( zptry(jpi,jpj,jpk) )
116         zptry(:,:,:) = 0._wp
117      ENDIF
118      !                          ! surface & bottom value : flux set to zero one for all
119      zwz(:,:, 1 ) = 0._wp           
120      zwx(:,:,jpk) = 0._wp   ;   zwy(:,:,jpk) = 0._wp    ;    zwz(:,:,jpk) = 0._wp
121      !
122      zwi(:,:,:) = 0._wp       
123      !
124      ! If adaptive vertical advection, check if it is needed on this PE at this time
125      IF( ln_zad_Aimp ) THEN
126         IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE.
127      END IF
128      ! If active adaptive vertical advection, build tridiagonal matrix
129      IF( ll_zAimp ) THEN
130         ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk))
131         DO_3D_00_00( 1, jpkm1 )
132            zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t(ji,jj,jk,Krhs)
133            zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t(ji,jj,jk,Krhs)
134            zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs)
135         END_3D
136      END IF
137      !
138      DO jn = 1, kjpt            !==  loop over the tracers  ==!
139         !
140         !        !==  upstream advection with initial mass fluxes & intermediate update  ==!
141         !                    !* upstream tracer flux in the i and j direction
142         DO_3D_10_10( 1, jpkm1 )
143            ! upstream scheme
144            zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )
145            zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) )
146            zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) )
147            zfm_vj = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) )
148            zwx(ji,jj,jk) = 0.5 * ( zfp_ui * pt(ji,jj,jk,jn,Kbb) + zfm_ui * pt(ji+1,jj  ,jk,jn,Kbb) )
149            zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji  ,jj+1,jk,jn,Kbb) )
150         END_3D
151         !                    !* upstream tracer flux in the k direction *!
152         DO_3D_11_11( 2, jpkm1 )
153            zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
154            zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
155            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)
156         END_3D
157         IF( ln_linssh ) THEN    ! top ocean value (only in linear free surface as zwz has been w-masked)
158            IF( ln_isfcav ) THEN             ! top of the ice-shelf cavities and at the ocean surface
159               DO_2D_11_11
160                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface
161               END_2D
162            ELSE                             ! no cavities: only at the ocean surface
163               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
164            ENDIF
165         ENDIF
166         !               
167         DO_3D_00_00( 1, jpkm1 )
168            !                             ! total intermediate advective trends
169            ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
170               &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
171               &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
172            !                             ! update and guess with monotonic sheme
173            pt(ji,jj,jk,jn,Krhs) =                     pt(ji,jj,jk,jn,Krhs) +        ztra   / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
174            zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
175         END_3D
176         
177         IF ( ll_zAimp ) THEN
178            CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
179            !
180            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
181            DO_3D_00_00( 2, jpkm1 )
182               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
183               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
184               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)
185               zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes
186            END_3D
187            DO_3D_00_00( 1, jpkm1 )
188               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
189                  &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
190            END_3D
191            !
192         END IF
193         !               
194         IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes)
195            ztrdx(:,:,:) = zwx(:,:,:)   ;   ztrdy(:,:,:) = zwy(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:)
196         END IF
197         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes)
198         IF( l_ptr )   zptry(:,:,:) = zwy(:,:,:) 
199         !
200         !        !==  anti-diffusive flux : high order minus low order  ==!
201         !
202         SELECT CASE( kn_fct_h )    !* horizontal anti-diffusive fluxes
203         !
204         CASE(  2  )                   !- 2nd order centered
205            DO_3D_10_10( 1, jpkm1 )
206               zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx(ji,jj,jk)
207               zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy(ji,jj,jk)
208            END_3D
209            !
210         CASE(  4  )                   !- 4th order centered
211            zltu(:,:,jpk) = 0._wp            ! Bottom value : flux set to zero
212            zltv(:,:,jpk) = 0._wp
213            DO jk = 1, jpkm1                 ! Laplacian
214               DO_2D_10_10
215                  ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
216                  ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
217               END_2D
218               DO_2D_00_00
219                  zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6
220                  zltv(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6
221               END_2D
222            END DO
223            CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1. , zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
224            !
225            DO_3D_10_10( 1, jpkm1 )
226               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
227               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
228               !                                                  ! C4 minus upstream advective fluxes
229               zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
230               zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
231            END_3D
232            !
233         CASE(  41 )                   !- 4th order centered       ==>>   !!gm coding attempt   need to be tested
234            ztu(:,:,jpk) = 0._wp             ! Bottom value : flux set to zero
235            ztv(:,:,jpk) = 0._wp
236            DO_3D_10_10( 1, jpkm1 )
237               ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
238               ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
239            END_3D
240            CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1. , ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn)
241            !
242            DO_3D_00_00( 1, jpkm1 )
243               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)
244               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
245               !                                                  ! C4 interpolation of T at u- & v-points (x2)
246               zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj  ,jk) - ztu(ji+1,jj  ,jk) )
247               zC4t_v =  zC2t_v + r1_6 * ( ztv(ji  ,jj-1,jk) - ztv(ji  ,jj+1,jk) )
248               !                                                  ! C4 minus upstream advective fluxes
249               zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
250               zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
251            END_3D
252            !
253         END SELECT
254         !                     
255         SELECT CASE( kn_fct_v )    !* vertical anti-diffusive fluxes (w-masked interior values)
256         !
257         CASE(  2  )                   !- 2nd order centered
258            DO_3D_00_00( 2, jpkm1 )
259               zwz(ji,jj,jk) =  (  pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
260                  &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk)
261            END_3D
262            !
263         CASE(  4  )                   !- 4th order COMPACT
264            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )   ! zwt = COMPACT interpolation of T at w-point
265            DO_3D_00_00( 2, jpkm1 )
266               zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
267            END_3D
268            !
269         END SELECT
270         IF( ln_linssh ) THEN    ! top ocean value: high order = upstream  ==>>  zwz=0
271            zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked
272         ENDIF
273         !         
274         IF ( ll_zAimp ) THEN
275            DO_3D_00_00( 1, jpkm1 )
276               !                             ! total intermediate advective trends
277               ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
278                  &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
279                  &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
280               ztw(ji,jj,jk)  = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
281            END_3D
282            !
283            CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
284            !
285            DO_3D_00_00( 2, jpkm1 )
286               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
287               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
288               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)
289            END_3D
290         END IF
291         !
292         CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1.,  zwz, 'W',  1. )
293         !
294         !        !==  monotonicity algorithm  ==!
295         !
296         CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx, zwy, zwz, zwi, p2dt )
297         !
298         !        !==  final trend with corrected fluxes  ==!
299         !
300         DO_3D_00_00( 1, jpkm1 )
301            ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
302               &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
303               &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
304            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm)
305            zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
306         END_3D
307         !
308         IF ( ll_zAimp ) THEN
309            !
310            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp
311            DO_3D_00_00( 2, jpkm1 )
312               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
313               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
314               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)
315               zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic
316            END_3D
317            DO_3D_00_00( 1, jpkm1 )
318               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
319                  &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
320            END_3D
321         END IF         
322         !
323         IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport
324            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< add anti-diffusive fluxes
325            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  !     to upstream fluxes
326            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  !
327            !
328            IF( l_trd ) THEN              ! trend diagnostics
329               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx, pU, pt(:,:,:,jn,Kmm) )
330               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy, pV, pt(:,:,:,jn,Kmm) )
331               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz, pW, pt(:,:,:,jn,Kmm) )
332            ENDIF
333            !                             ! heat/salt transport
334            IF( l_hst )   CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
335            !
336         ENDIF
337         IF( l_ptr ) THEN              ! "Poleward" transports
338            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< add anti-diffusive fluxes
339            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
340         ENDIF
341         !
342      END DO                     ! end of tracer loop
343      !
344      IF ( ll_zAimp ) THEN
345         DEALLOCATE( zwdia, zwinf, zwsup )
346      ENDIF
347      IF( l_trd .OR. l_hst ) THEN
348         DEALLOCATE( ztrdx, ztrdy, ztrdz )
349      ENDIF
350      IF( l_ptr ) THEN
351         DEALLOCATE( zptry )
352      ENDIF
353      !
354   END SUBROUTINE tra_adv_fct
355
356
357   SUBROUTINE nonosc( Kmm, pbef, paa, pbb, pcc, paft, p2dt )
358      !!---------------------------------------------------------------------
359      !!                    ***  ROUTINE nonosc  ***
360      !!     
361      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
362      !!       scheme and the before field by a nonoscillatory algorithm
363      !!
364      !! **  Method  :   ... ???
365      !!       warning : pbef and paft must be masked, but the boundaries
366      !!       conditions on the fluxes are not necessary zalezak (1979)
367      !!       drange (1995) multi-dimensional forward-in-time and upstream-
368      !!       in-space based differencing for fluid
369      !!----------------------------------------------------------------------
370      INTEGER                          , INTENT(in   ) ::   Kmm             ! time level index
371      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! tracer time-step
372      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
373      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb, pcc   ! monotonic fluxes in the 3 directions
374      !
375      INTEGER  ::   ji, jj, jk   ! dummy loop indices
376      INTEGER  ::   ikm1         ! local integer
377      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars
378      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
379      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo
380      !!----------------------------------------------------------------------
381      !
382      zbig  = 1.e+40_wp
383      zrtrn = 1.e-15_wp
384      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
385
386      ! Search local extrema
387      ! --------------------
388      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
389      zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ),   &
390         &        paft * tmask - zbig * ( 1._wp - tmask )  )
391      zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ),   &
392         &        paft * tmask + zbig * ( 1._wp - tmask )  )
393
394      DO jk = 1, jpkm1
395         ikm1 = MAX(jk-1,1)
396         DO_2D_00_00
397
398            ! search maximum in neighbourhood
399            zup = MAX(  zbup(ji  ,jj  ,jk  ),   &
400               &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   &
401               &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ),   &
402               &        zbup(ji  ,jj  ,ikm1), zbup(ji  ,jj  ,jk+1)  )
403
404            ! search minimum in neighbourhood
405            zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   &
406               &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   &
407               &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ),   &
408               &        zbdo(ji  ,jj  ,ikm1), zbdo(ji  ,jj  ,jk+1)  )
409
410            ! positive part of the flux
411            zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
412               & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   &
413               & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
414
415            ! negative part of the flux
416            zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
417               & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   &
418               & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
419
420            ! up & down beta terms
421            zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt
422            zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
423            zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt
424         END_2D
425      END DO
426      CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1. , zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
427
428      ! 3. monotonic flux in the i & j direction (paa & pbb)
429      ! ----------------------------------------
430      DO_3D_00_00( 1, jpkm1 )
431         zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
432         zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
433         zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) )
434         paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
435
436         zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
437         zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
438         zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) )
439         pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
440
441! monotonic flux in the k direction, i.e. pcc
442! -------------------------------------------
443         za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
444         zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
445         zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
446         pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
447      END_3D
448      CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1. , pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
449      !
450   END SUBROUTINE nonosc
451
452
453   SUBROUTINE interp_4th_cpt_org( pt_in, pt_out )
454      !!----------------------------------------------------------------------
455      !!                  ***  ROUTINE interp_4th_cpt_org  ***
456      !!
457      !! **  Purpose :   Compute the interpolation of tracer at w-point
458      !!
459      !! **  Method  :   4th order compact interpolation
460      !!----------------------------------------------------------------------
461      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt_in    ! now tracer fields
462      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pt_out   ! now tracer field interpolated at w-pts
463      !
464      INTEGER :: ji, jj, jk   ! dummy loop integers
465      REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
466      !!----------------------------------------------------------------------
467     
468      DO_3D_11_11( 3, jpkm1 )
469         zwd (ji,jj,jk) = 4._wp
470         zwi (ji,jj,jk) = 1._wp
471         zws (ji,jj,jk) = 1._wp
472         zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
473         !
474         IF( tmask(ji,jj,jk+1) == 0._wp) THEN   ! Switch to second order centered at bottom
475            zwd (ji,jj,jk) = 1._wp
476            zwi (ji,jj,jk) = 0._wp
477            zws (ji,jj,jk) = 0._wp
478            zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )   
479         ENDIF
480      END_3D
481      !
482      jk = 2                                          ! Switch to second order centered at top
483      DO_2D_11_11
484         zwd (ji,jj,jk) = 1._wp
485         zwi (ji,jj,jk) = 0._wp
486         zws (ji,jj,jk) = 0._wp
487         zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
488      END_2D
489      !
490      !                       !==  tridiagonal solve  ==!
491      DO_2D_11_11
492         zwt(ji,jj,2) = zwd(ji,jj,2)
493      END_2D
494      DO_3D_11_11( 3, jpkm1 )
495         zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
496      END_3D
497      !
498      DO_2D_11_11
499         pt_out(ji,jj,2) = zwrm(ji,jj,2)
500      END_2D
501      DO_3D_11_11( 3, jpkm1 )
502         pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
503      END_3D
504
505      DO_2D_11_11
506         pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
507      END_2D
508      DO_3DS_11_11( jpk-2, 2, -1 )
509         pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
510      END_3D
511      !   
512   END SUBROUTINE interp_4th_cpt_org
513   
514
515   SUBROUTINE interp_4th_cpt( pt_in, pt_out )
516      !!----------------------------------------------------------------------
517      !!                  ***  ROUTINE interp_4th_cpt  ***
518      !!
519      !! **  Purpose :   Compute the interpolation of tracer at w-point
520      !!
521      !! **  Method  :   4th order compact interpolation
522      !!----------------------------------------------------------------------
523      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt_in    ! field at t-point
524      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pt_out   ! field interpolated at w-point
525      !
526      INTEGER ::   ji, jj, jk   ! dummy loop integers
527      INTEGER ::   ikt, ikb     ! local integers
528      REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
529      !!----------------------------------------------------------------------
530      !
531      !                      !==  build the three diagonal matrix & the RHS  ==!
532      !
533      DO_3D_00_00( 3, jpkm1 )
534         zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp                 !       diagonal
535         zwi (ji,jj,jk) =         wmask(ji,jj,jk)                         ! lower diagonal
536         zws (ji,jj,jk) =         wmask(ji,jj,jk)                         ! upper diagonal
537         zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk)                     &   ! RHS
538            &           *       ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) )
539      END_3D
540      !
541!!gm
542!      SELECT CASE( kbc )               !* boundary condition
543!      CASE( np_NH   )   ! Neumann homogeneous at top & bottom
544!      CASE( np_CEN2 )   ! 2nd order centered  at top & bottom
545!      END SELECT
546!!gm 
547      !
548      IF ( ln_isfcav ) THEN            ! set level two values which may not be set in ISF case
549         zwd(:,:,2) = 1._wp  ;  zwi(:,:,2) = 0._wp  ;  zws(:,:,2) = 0._wp  ;  zwrm(:,:,2) = 0._wp
550      END IF
551      !
552      DO_2D_00_00
553         ikt = mikt(ji,jj) + 1            ! w-point below the 1st  wet point
554         ikb = MAX(mbkt(ji,jj), 2)        !     -   above the last wet point
555         !
556         zwd (ji,jj,ikt) = 1._wp          ! top
557         zwi (ji,jj,ikt) = 0._wp
558         zws (ji,jj,ikt) = 0._wp
559         zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,ikt-1) + pt_in(ji,jj,ikt) )
560         !
561         zwd (ji,jj,ikb) = 1._wp          ! bottom
562         zwi (ji,jj,ikb) = 0._wp
563         zws (ji,jj,ikb) = 0._wp
564         zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) )           
565      END_2D
566      !
567      !                       !==  tridiagonal solver  ==!
568      !
569      DO_2D_00_00
570         zwt(ji,jj,2) = zwd(ji,jj,2)
571      END_2D
572      DO_3D_00_00( 3, jpkm1 )
573         zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
574      END_3D
575      !
576      DO_2D_00_00
577         pt_out(ji,jj,2) = zwrm(ji,jj,2)
578      END_2D
579      DO_3D_00_00( 3, jpkm1 )
580         pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
581      END_3D
582
583      DO_2D_00_00
584         pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
585      END_2D
586      DO_3DS_00_00( jpk-2, 2, -1 )
587         pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
588      END_3D
589      !   
590   END SUBROUTINE interp_4th_cpt
591
592
593   SUBROUTINE tridia_solver( pD, pU, pL, pRHS, pt_out , klev )
594      !!----------------------------------------------------------------------
595      !!                  ***  ROUTINE tridia_solver  ***
596      !!
597      !! **  Purpose :   solve a symmetric 3diagonal system
598      !!
599      !! **  Method  :   solve M.t_out = RHS(t)  where M is a tri diagonal matrix ( jpk*jpk )
600      !!     
601      !!             ( D_1 U_1  0   0   0  )( t_1 )   ( RHS_1 )
602      !!             ( L_2 D_2 U_2  0   0  )( t_2 )   ( RHS_2 )
603      !!             (  0  L_3 D_3 U_3  0  )( t_3 ) = ( RHS_3 )
604      !!             (        ...          )( ... )   ( ...  )
605      !!             (  0   0   0  L_k D_k )( t_k )   ( RHS_k )
606      !!     
607      !!        M is decomposed in the product of an upper and lower triangular matrix.
608      !!        The tri-diagonals matrix is given as input 3D arrays:   pD, pU, pL
609      !!        (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal).
610      !!        The solution is pta.
611      !!        The 3d array zwt is used as a work space array.
612      !!----------------------------------------------------------------------
613      REAL(wp),DIMENSION(:,:,:), INTENT(in   ) ::   pD, pU, PL    ! 3-diagonal matrix
614      REAL(wp),DIMENSION(:,:,:), INTENT(in   ) ::   pRHS          ! Right-Hand-Side
615      REAL(wp),DIMENSION(:,:,:), INTENT(  out) ::   pt_out        !!gm field at level=F(klev)
616      INTEGER                  , INTENT(in   ) ::   klev          ! =1 pt_out at w-level
617      !                                                           ! =0 pt at t-level
618      INTEGER ::   ji, jj, jk   ! dummy loop integers
619      INTEGER ::   kstart       ! local indices
620      REAL(wp),DIMENSION(jpi,jpj,jpk) ::   zwt   ! 3D work array
621      !!----------------------------------------------------------------------
622      !
623      kstart =  1  + klev
624      !
625      DO_2D_00_00
626         zwt(ji,jj,kstart) = pD(ji,jj,kstart)
627      END_2D
628      DO_3D_00_00( kstart+1, jpkm1 )
629         zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1)
630      END_3D
631      !
632      DO_2D_00_00
633         pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart)
634      END_2D
635      DO_3D_00_00( kstart+1, jpkm1 )
636         pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
637      END_3D
638
639      DO_2D_00_00
640         pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
641      END_2D
642      DO_3DS_00_00( jpk-2, kstart, -1 )
643         pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
644      END_3D
645      !
646   END SUBROUTINE tridia_solver
647
648   !!======================================================================
649END MODULE traadv_fct
Note: See TracBrowser for help on using the repository browser.