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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_fct.F90 @ 10946

Last change on this file since 10946 was 10946, checked in by acc, 5 years ago

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