source: NEMO/trunk/src/OCE/TRA/traadv_fct.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 13 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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