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.
icedyn_adv_umx.F90 in NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/MY_SRC – NEMO

source: NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/MY_SRC/icedyn_adv_umx.F90 @ 10281

Last change on this file since 10281 was 10281, checked in by clem, 5 years ago

advection routine containing all the tests implemented (Zalesak, Devore, etc)
File size: 60.1 KB
Line 
1MODULE icedyn_adv_umx
2   !!==============================================================================
3   !!                       ***  MODULE  icedyn_adv_umx  ***
4   !! sea-ice : advection using the ULTIMATE-MACHO scheme
5   !!==============================================================================
6   !! History :  3.6  !  2014-11  (C. Rousset, G. Madec)  Original code
7   !!            4.0  !  2018     (many people)           SI3 [aka Sea Ice cube]
8   !!----------------------------------------------------------------------
9#if defined key_si3
10   !!----------------------------------------------------------------------
11   !!   'key_si3'                                       SI3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!   ice_dyn_adv_umx   : update the tracer trend with the 3D advection trends using a TVD scheme
14   !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders
15   !!   macho             : ???
16   !!   nonosc            : compute monotonic tracer fluxes by a non-oscillatory algorithm
17   !!----------------------------------------------------------------------
18   USE phycst         ! physical constant
19   USE dom_oce        ! ocean domain
20   USE sbc_oce , ONLY : nn_fsbc   ! update frequency of surface boundary condition
21   USE ice            ! sea-ice variables
22   USE icevar         ! sea-ice: operations
23   !
24   USE in_out_manager ! I/O manager
25   USE lib_mpp        ! MPP library
26   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
27   USE lbclnk         ! lateral boundary conditions (or mpp links)
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   ice_dyn_adv_umx   ! called by icedyn_adv.F90
33     
34   REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6
35   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
36
37   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_a
38   
39   ! alternate directions for upstream
40   ! clem: it gives results for Lipscomb test that are the same as "ll_upsxy=false"
41   LOGICAL :: ll_upsxy = .FALSE.
42   
43   ! prelimiter
44   LOGICAL :: ll_prelimiter = .TRUE.
45   LOGICAL :: ll_prelimiter_zalesak = .TRUE.   ! from: Zalesak(1979) eq. 14 => better than Devore
46   LOGICAL :: ll_prelimiter_devore  = .FALSE.  ! from: Devore eq. 11
47
48   ! iterate on the limiter (only nonosc)
49   LOGICAL :: ll_limiter_it2 = .TRUE.
50   
51
52   !! * Substitutions
53#  include "vectopt_loop_substitute.h90"
54   !!----------------------------------------------------------------------
55   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
56   !! $Id: icedyn_adv_umx.F90 4499 2014-02-18 15:14:31Z timgraham $
57   !! Software governed by the CeCILL licence     (./LICENSE)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice,  &
62      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
63      !!----------------------------------------------------------------------
64      !!                  ***  ROUTINE ice_dyn_adv_umx  ***
65      !!
66      !! **  Purpose :   Compute the now trend due to total advection of
67      !!                 tracers and add it to the general trend of tracer equations
68      !!                 using an "Ultimate-Macho" scheme
69      !!
70      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
71      !!----------------------------------------------------------------------
72      INTEGER                     , INTENT(in   ) ::   kn_umx     ! order of the scheme (1-5=UM or 20=CEN2)
73      INTEGER                     , INTENT(in   ) ::   kt         ! time step
74      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
75      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
76      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
77      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
78      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
79      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
80      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
81      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
82      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
83      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
84      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
85      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
86      !
87      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
88      INTEGER  ::   icycle                  ! number of sub-timestep for the advection
89      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers
90      REAL(wp) ::   zcfl , zdt
91      REAL(wp) ::   zeps = 0.0_wp           ! shift in concentration to avoid division by 0
92      !                                     !    must be >= 0.01 and the best seems to be 0.1
93      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx, zcu_box, zcv_box, zfu, zfv
94      REAL(wp), DIMENSION(jpi,jpj) ::   zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, zh_ip 
95      !!----------------------------------------------------------------------
96      !
97      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
98      !
99      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !       
100      zcfl  =            MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
101      zcfl  = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
102      IF( lk_mpp )   CALL mpp_max( zcfl )
103
104      IF( zcfl > 0.5 ) THEN   ;   icycle = 2 
105      ELSE                    ;   icycle = 1 
106      ENDIF
107
108      zdt = rdt_ice / REAL(icycle)
109
110      ! --- transport --- !
111      zudy(:,:) = pu_ice(:,:) * e2u(:,:)
112      zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
113
114      ! --- define velocity for advection: u*grad(H) --- !
115      DO jj = 2, jpjm1
116         DO ji = fs_2, fs_jpim1
117            IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
118            ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj)
119            ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj)
120            ENDIF
121
122            IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
123            ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1)
124            ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  )
125            ENDIF
126         END DO
127      END DO
128
129      IF(.NOT. ALLOCATED(z1_a )) ALLOCATE(z1_a (jpi,jpj))
130      !---------------!
131      !== advection ==!
132      !---------------!
133      DO jt = 1, icycle
134
135         zamsk = 1._wp
136         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) )   ! Open water area
137
138         DO jl = 1, jpl
139            ! to avoid a problem with the limiter nonosc when A gets close to 0
140            pa_i(:,:,jl) = pa_i(:,:,jl) + zeps * tmask(:,:,1)
141            !
142            WHERE( pa_i(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_i(:,:,jl)
143            ELSEWHERE                        ;   z1_a(:,:) = 0.   !!;   pa_i(:,:,jl) = 0. ; pv_i(:,:,jl) = 0.
144            END WHERE
145            !
146            zamsk = 1._wp
147            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zfu, zfv ) ! Ice area
148            !!zfu = zudy ; zfv = zvdx
149            !!CALL adv_umx( kn_umx, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, pv_i(:,:,jl), pv_i(:,:,jl) )
150            zamsk = 0._wp ; zh_i(:,:) = pv_i (:,:,jl) * z1_a(:,:)
151            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_i(:,:), pv_i (:,:,jl) )              ! Ice volume
152            zamsk = 0._wp ; zh_s(:,:) = pv_s (:,:,jl) * z1_a(:,:)
153            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_s(:,:), pv_s (:,:,jl) )              ! Snw volume
154            zamsk = 0._wp ; zs_i(:,:) = psv_i(:,:,jl) * z1_a(:,:)
155            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zs_i(:,:), psv_i(:,:,jl) )              ! Salt content
156            zamsk = 0._wp ; zo_i(:,:) = poa_i(:,:,jl) * z1_a(:,:)
157            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zo_i(:,:), poa_i(:,:,jl) )              ! Age content
158            DO jk = 1, nlay_i
159               zamsk = 0._wp ; ze_i(:,:) = pe_i(:,:,jk,jl) * z1_a(:,:)
160               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu, zfv, zcu_box, zcv_box, ze_i(:,:), pe_i(:,:,jk,jl) )           ! Ice heat content
161            END DO
162            DO jk = 1, nlay_s
163               zamsk = 0._wp ; ze_s(:,:) = pe_s(:,:,jk,jl) * z1_a(:,:)
164               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu, zfv, zcu_box, zcv_box, ze_s(:,:), pe_s(:,:,jk,jl) )           ! Snw heat content
165            END DO
166            !
167            IF ( ln_pnd_H12 ) THEN   ! melt ponds (must be the last ones to be advected because of z1_a...)
168               ! to avoid a problem with the limiter nonosc when A gets close to 0
169               pa_ip(:,:,jl) = pa_ip(:,:,jl) + zeps * tmask(:,:,1)
170               !
171               WHERE( pa_ip(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_ip(:,:,jl)
172               ELSEWHERE                         ;   z1_a(:,:) = 0.
173               END WHERE
174               !
175               zamsk = 1._wp
176               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zfu, zfv ) ! mp fraction
177               zamsk = 0._wp ; zh_ip(:,:) = pv_ip(:,:,jl) * z1_a(:,:)
178               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_ip(:,:), pv_ip(:,:,jl) )              ! mp volume
179            ENDIF
180            !
181            ! to avoid a problem with the limiter nonosc when A gets close to 0
182            DO jj = 2, jpjm1
183               DO ji = fs_2, fs_jpim1
184                  !pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps ) * tmask(ji,jj,1)
185                  pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps &
186                     &             + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt &
187                     &             ) * tmask(ji,jj,1)
188                  IF ( ln_pnd_H12 ) THEN   ! melt ponds
189                     pa_ip(ji,jj,jl) = ( pa_ip(ji,jj,jl) - zeps &
190                        &              + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt &
191                        &              ) * tmask(ji,jj,1)
192                  ENDIF
193              END DO
194            END DO
195            CALL lbc_lnk( pa_i(:,:,jl), 'T',  1. )
196            !
197         END DO
198
199      END DO
200      !
201   END SUBROUTINE ice_dyn_adv_umx
202
203   
204   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pfu, pfv )
205      !!----------------------------------------------------------------------
206      !!                  ***  ROUTINE adv_umx  ***
207      !!
208      !! **  Purpose :   Compute the now trend due to total advection of
209      !!       tracers and add it to the general trend of tracer equations
210      !!
211      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
212      !!       corrected flux (monotonic correction)
213      !!       note: - this advection scheme needs a leap-frog time scheme
214      !!
215      !! ** Action : - pt  the after advective tracer
216      !!----------------------------------------------------------------------
217      REAL(wp)                    , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0)
218      INTEGER                     , INTENT(in   )           ::   kn_umx         ! order of the scheme (1-5=UM or 20=CEN2)
219      INTEGER                     , INTENT(in   )           ::   jt             ! number of sub-iteration
220      INTEGER                     , INTENT(in   )           ::   kt             ! number of iteration
221      REAL(wp)                    , INTENT(in   )           ::   pdt            ! tracer time-step
222      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2
223      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u
224      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity
225      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pt             ! tracer field
226      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptc            ! tracer content field
227      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pfu, pfv       ! high order fluxes
228      !
229      INTEGER  ::   ji, jj           ! dummy loop indices 
230      REAL(wp) ::   ztra             ! local scalar
231      INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid)
232      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt, zz1_a
233      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups, zt_ups   ! only for nonosc
234      !!----------------------------------------------------------------------
235      !
236      !  upstream (_ups) advection with initial mass fluxes
237      ! ---------------------------------------------------
238      IF( .NOT. ll_upsxy ) THEN
239
240         ! fluxes
241         DO jj = 1, jpjm1
242            DO ji = 1, fs_jpim1
243               zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj)
244               zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1)
245            END DO
246         END DO
247
248      ELSE
249         !
250         IF    ( pamsk == 1. ) THEN ; zz1_a(:,:) = 1._wp     ;        ! 1   if advection of A
251         ELSEIF( pamsk == 0. ) THEN ; zz1_a(:,:) = z1_a(:,:) ; ENDIF  ! 1/A if advection of V
252         !
253         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
254            ! flux in x-direction
255            DO jj = 1, jpjm1
256               DO ji = 1, fs_jpim1
257                  zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj)
258               END DO
259            END DO
260            ! first guess of tracer content from u-flux
261            DO jj = 2, jpjm1
262               DO ji = fs_2, fs_jpim1   ! vector opt.
263                  zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * zz1_a(ji,jj)
264               END DO
265            END DO
266            CALL lbc_lnk( zpt, 'T', 1. )
267            !
268            ! flux in y-direction
269            DO jj = 1, jpjm1
270               DO ji = 1, fs_jpim1
271                  zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1)
272               END DO
273            END DO
274            !
275         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
276            ! flux in y-direction
277            DO jj = 1, jpjm1
278               DO ji = 1, fs_jpim1
279                  zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1)
280               END DO
281            END DO
282            ! first guess of tracer content from v-flux
283            DO jj = 2, jpjm1
284               DO ji = fs_2, fs_jpim1   ! vector opt.
285                  zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * zz1_a(ji,jj) 
286               END DO
287            END DO
288            CALL lbc_lnk( zpt, 'T', 1. )
289            !
290            ! flux in x-direction
291            DO jj = 1, jpjm1
292               DO ji = 1, fs_jpim1
293                  zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj)
294               END DO
295            END DO
296            !
297         ENDIF
298         
299      ENDIF
300      ! guess after content field with upstream scheme
301      DO jj = 2, jpjm1
302         DO ji = fs_2, fs_jpim1
303            ztra          = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   &
304               &                + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1) ) * r1_e1e2t(ji,jj)
305            zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
306         END DO
307      END DO
308      CALL lbc_lnk( zt_ups, 'T', 1. )
309     
310      ! High order (_ho) fluxes
311      ! -----------------------
312      SELECT CASE( kn_umx )
313         !
314      CASE ( 20 )                          !== centered second order ==!
315         !
316         CALL cen2( kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  &
317            &       zt_ups, zfu_ups, zfv_ups )
318         !
319      CASE ( 1:5 )                         !== 1st to 5th order ULTIMATE-MACHO scheme ==!
320         !
321         CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho,  &
322            &        zt_ups, zfu_ups, zfv_ups )
323         !
324      END SELECT
325         
326      ! output high order fluxes
327      ! ------------------------
328      IF( PRESENT(pfu) ) THEN
329         DO jj = 1, jpjm1
330            DO ji = 1, fs_jpim1
331               pfu(ji,jj) =  zfu_ho(ji,jj)
332               pfv(ji,jj) =  zfv_ho(ji,jj)
333            END DO
334         END DO
335         !!CALL lbc_lnk( pfu, 'U',  -1. ) ! clem: not needed I think
336         !!CALL lbc_lnk( pfv, 'V',  -1. )
337      ENDIF
338     
339      ! final trend with corrected fluxes
340      ! ------------------------------------
341      DO jj = 2, jpjm1
342         DO ji = fs_2, fs_jpim1 
343            ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  &  !        Div(uaH) or        Div(ua)
344               &         ) * r1_e1e2t(ji,jj) 
345            ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
346         END DO
347      END DO
348      CALL lbc_lnk( ptc, 'T',  1. )
349      !
350   END SUBROUTINE adv_umx
351
352   SUBROUTINE cen2( kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, &
353      &             pt_ups, pfu_ups, pfv_ups )
354      !!---------------------------------------------------------------------
355      !!                    ***  ROUTINE macho  ***
356      !!     
357      !! **  Purpose :   compute 
358      !!
359      !! **  Method  :   ... ???
360      !!                 TIM = transient interpolation Modeling
361      !!
362      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
363      !!----------------------------------------------------------------------
364      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter
365      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration
366      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration
367      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step
368      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt               ! tracer fields
369      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
370      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components
371      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc              ! tracer content at before time step
372      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes
373      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content
374      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes
375      !
376      INTEGER  ::   ji, jj    ! dummy loop indices
377      LOGICAL  ::   ll_xy = .TRUE.
378      REAL(wp), DIMENSION(jpi,jpj) ::   zzt
379      !!----------------------------------------------------------------------
380      !
381      IF( .NOT.ll_xy ) THEN   !-- no alternate directions --!
382         !
383         DO jj = 1, jpjm1
384            DO ji = 1, fs_jpim1
385               pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
386               pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
387            END DO
388         END DO
389         IF    ( kn_limiter == 1 ) THEN
390            CALL nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
391         ELSEIF( kn_limiter == 2 ) THEN
392            CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
393            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
394         ELSEIF( kn_limiter == 3 ) THEN
395            CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
396            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
397         ENDIF
398         !
399      ELSE                    !-- alternate directions --!
400         !
401         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
402            !
403            ! flux in x-direction
404            DO jj = 1, jpjm1
405               DO ji = 1, fs_jpim1
406                  pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
407               END DO
408            END DO
409            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
410            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
411
412            ! first guess of tracer content from u-flux
413            DO jj = 2, jpjm1
414               DO ji = fs_2, fs_jpim1   ! vector opt.
415                  zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
416               END DO
417            END DO
418            CALL lbc_lnk( zzt, 'T', 1. )
419
420            ! flux in y-direction
421            DO jj = 1, jpjm1
422               DO ji = 1, fs_jpim1
423                  pfv_ho(ji,jj) = 0.5 * pv(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) )
424               END DO
425            END DO
426            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
427            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
428
429         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
430            !
431            ! flux in y-direction
432            DO jj = 1, jpjm1
433               DO ji = 1, fs_jpim1
434                  pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
435               END DO
436            END DO
437            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
438            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
439            !
440            ! first guess of tracer content from v-flux
441            DO jj = 2, jpjm1
442               DO ji = fs_2, fs_jpim1   ! vector opt.
443                  zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
444               END DO
445            END DO
446            CALL lbc_lnk( zzt, 'T', 1. )
447            !
448            ! flux in x-direction
449            DO jj = 1, jpjm1
450               DO ji = 1, fs_jpim1
451                  pfu_ho(ji,jj) = 0.5 * pu(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) )
452               END DO
453            END DO
454            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
455            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
456
457         ENDIF
458         IF( kn_limiter == 1 )   CALL nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
459         
460      ENDIF
461   
462   END SUBROUTINE cen2
463
464   
465   SUBROUTINE macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, pt_u, pt_v, pfu_ho, pfv_ho, &
466      &              pt_ups, pfu_ups, pfv_ups )
467      !!---------------------------------------------------------------------
468      !!                    ***  ROUTINE macho  ***
469      !!     
470      !! **  Purpose :   compute 
471      !!
472      !! **  Method  :   ... ???
473      !!                 TIM = transient interpolation Modeling
474      !!
475      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
476      !!----------------------------------------------------------------------
477      REAL(wp)                    , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
478      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter
479      INTEGER                     , INTENT(in   ) ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2)
480      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration
481      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration
482      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step
483      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt               ! tracer fields
484      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
485      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components
486      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox     ! upstream velocity
487      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc              ! tracer content at before time step
488      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u, pt_v       ! tracer at u- and v-points
489      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes
490      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content
491      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes
492      !
493      INTEGER  ::   ji, jj    ! dummy loop indices
494      REAL(wp) ::   ztra
495      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zzfu_ho, zzfv_ho
496      !!----------------------------------------------------------------------
497      !
498      IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
499         !
500         !                                                        !--  ultimate interpolation of pt at u-point  --!
501         CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
502         !                                                        !--  limiter in x --!
503         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
504         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
505         !                                                        !--  advective form update in zzt  --!
506         DO jj = 2, jpjm1
507            DO ji = fs_2, fs_jpim1   ! vector opt.
508               zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  &
509                  &                   - pt   (ji,jj) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk
510               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
511            END DO
512         END DO
513         CALL lbc_lnk( zzt, 'T', 1. )
514         !                                                        !--  ultimate interpolation of pt at v-point  --!
515         CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho )
516         !                                                        !--  limiter in y --!
517         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
518         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
519         !
520      ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
521         !
522         !                                                        !--  ultimate interpolation of pt at v-point  --!
523         CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
524         !                                                        !--  limiter in y --!
525         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
526         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
527         !                                                        !--  advective form update in zzt  --!
528         DO jj = 2, jpjm1
529            DO ji = fs_2, fs_jpim1
530               zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  &
531                  &                    - pt  (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk
532               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
533            END DO
534         END DO
535         CALL lbc_lnk( zzt, 'T', 1. )
536         !                                                        !--  ultimate interpolation of pt at u-point  --!
537         CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho )
538         !                                                        !--  limiter in x --!
539         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
540         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
541         !
542      ENDIF
543      IF( kn_limiter == 1 ) THEN
544         IF( .NOT. ll_limiter_it2 ) THEN
545            CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
546         ELSE
547            zzfu_ho(:,:) = pfu_ho(:,:)
548            zzfv_ho(:,:) = pfv_ho(:,:)
549            ! 1st iteration of nonosc (limit the flux with the upstream solution)
550            CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
551            ! guess after content field with high order
552            DO jj = 2, jpjm1
553               DO ji = fs_2, fs_jpim1
554                  ztra       = - ( zzfu_ho(ji,jj) - zzfu_ho(ji-1,jj) + zzfv_ho(ji,jj) - zzfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj)
555                  zzt(ji,jj) =   ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
556               END DO
557            END DO
558            CALL lbc_lnk( zzt, 'T', 1. )
559            ! 2nd iteration of nonosc (limit the flux with the limited high order solution)
560            CALL nonosc_2d ( pdt, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
561         ENDIF
562      ENDIF
563      !
564   END SUBROUTINE macho
565
566
567   SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
568      !!---------------------------------------------------------------------
569      !!                    ***  ROUTINE ultimate_x  ***
570      !!     
571      !! **  Purpose :   compute 
572      !!
573      !! **  Method  :   ... ???
574      !!                 TIM = transient interpolation Modeling
575      !!
576      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
577      !!----------------------------------------------------------------------
578      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
579      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
580      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu        ! ice i-velocity component
581      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc       ! ice i-velocity * A component
582      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
583      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u      ! tracer at u-point
584      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho    ! high order flux
585      !
586      INTEGER  ::   ji, jj             ! dummy loop indices
587      REAL(wp) ::   zcu, zdx2, zdx4    !   -      -
588      REAL(wp), DIMENSION(jpi,jpj) ::   ztu1, ztu2, ztu3, ztu4
589      !!----------------------------------------------------------------------
590      !
591      !                                                     !--  Laplacian in i-direction  --!
592      DO jj = 2, jpjm1         ! First derivative (gradient)
593         DO ji = 1, fs_jpim1
594            ztu1(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
595         END DO
596         !                     ! Second derivative (Laplacian)
597         DO ji = fs_2, fs_jpim1
598            ztu2(ji,jj) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj)
599         END DO
600      END DO
601      CALL lbc_lnk( ztu2, 'T', 1. )
602      !
603      !                                                     !--  BiLaplacian in i-direction  --!
604      DO jj = 2, jpjm1         ! Third derivative
605         DO ji = 1, fs_jpim1
606            ztu3(ji,jj) = ( ztu2(ji+1,jj) - ztu2(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
607         END DO
608         !                     ! Fourth derivative
609         DO ji = fs_2, fs_jpim1
610            ztu4(ji,jj) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj)
611         END DO
612      END DO
613      CALL lbc_lnk( ztu4, 'T', 1. )
614      !
615      !
616      SELECT CASE (kn_umx )
617      !
618      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
619         !       
620         DO jj = 2, jpjm1
621            DO ji = 1, fs_jpim1   ! vector opt.
622               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                              pt(ji+1,jj) + pt(ji,jj)   &
623                  &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )
624            END DO
625         END DO
626         !
627      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
628         !
629         DO jj = 2, jpjm1
630            DO ji = 1, fs_jpim1   ! vector opt.
631               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
632               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                                   pt(ji+1,jj) + pt(ji,jj)   &
633                  &                                               -              zcu   * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
634            END DO
635         END DO
636         ! 
637      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
638         !
639         DO jj = 2, jpjm1
640            DO ji = 1, fs_jpim1   ! vector opt.
641               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
642               zdx2 = e1u(ji,jj) * e1u(ji,jj)
643!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
644               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                         pt  (ji+1,jj) + pt  (ji,jj)        &
645                  &                                               -              zcu   * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   &
646                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj) + ztu2(ji,jj)        &
647                  &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
648            END DO
649         END DO
650         !
651      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
652         !
653         DO jj = 2, jpjm1
654            DO ji = 1, fs_jpim1   ! vector opt.
655               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
656               zdx2 = e1u(ji,jj) * e1u(ji,jj)
657!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
658               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                   pt  (ji+1,jj) + pt  (ji,jj)        &
659                  &                                               -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   &
660                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj) + ztu2(ji,jj)        &
661                  &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
662            END DO
663         END DO
664         !
665      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
666         !
667         DO jj = 2, jpjm1
668            DO ji = 1, fs_jpim1   ! vector opt.
669               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
670               zdx2 = e1u(ji,jj) * e1u(ji,jj)
671!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
672               zdx4 = zdx2 * zdx2
673               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (               (                   pt  (ji+1,jj) + pt  (ji,jj)       &
674                  &                                                     -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) ) )   &
675                  &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj) + ztu2(ji,jj)       &
676                  &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) )   &
677                  &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj) + ztu4(ji,jj)       &
678                  &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj) - ztu4(ji,jj) ) ) )
679            END DO
680         END DO
681         !
682      END SELECT
683      !                                                     !-- High order flux in i-direction  --!
684      DO jj = 1, jpjm1
685         DO ji = 1, fs_jpim1   ! vector opt.
686            pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj)
687         END DO
688      END DO
689      !
690   END SUBROUTINE ultimate_x
691   
692 
693   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
694      !!---------------------------------------------------------------------
695      !!                    ***  ROUTINE ultimate_y  ***
696      !!     
697      !! **  Purpose :   compute 
698      !!
699      !! **  Method  :   ... ???
700      !!                 TIM = transient interpolation Modeling
701      !!
702      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
703      !!----------------------------------------------------------------------
704      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
705      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
706      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pv        ! ice j-velocity component
707      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvc       ! ice j-velocity*A component
708      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
709      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_v      ! tracer at v-point
710      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfv_ho    ! high order flux
711      !
712      INTEGER  ::   ji, jj       ! dummy loop indices
713      REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
714      REAL(wp), DIMENSION(jpi,jpj) ::   ztv1, ztv2, ztv3, ztv4
715      !!----------------------------------------------------------------------
716      !
717      !                                                     !--  Laplacian in j-direction  --!
718      DO jj = 1, jpjm1         ! First derivative (gradient)
719         DO ji = fs_2, fs_jpim1
720            ztv1(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
721         END DO
722      END DO
723      DO jj = 2, jpjm1         ! Second derivative (Laplacian)
724         DO ji = fs_2, fs_jpim1
725            ztv2(ji,jj) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj)
726         END DO
727      END DO
728      CALL lbc_lnk( ztv2, 'T', 1. )
729      !
730      !                                                     !--  BiLaplacian in j-direction  --!
731      DO jj = 1, jpjm1         ! First derivative
732         DO ji = fs_2, fs_jpim1
733            ztv3(ji,jj) = ( ztv2(ji,jj+1) - ztv2(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
734         END DO
735      END DO
736      DO jj = 2, jpjm1         ! Second derivative
737         DO ji = fs_2, fs_jpim1
738            ztv4(ji,jj) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj)
739         END DO
740      END DO
741      CALL lbc_lnk( ztv4, 'T', 1. )
742      !
743      !
744      SELECT CASE (kn_umx )
745      !
746      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
747         DO jj = 1, jpjm1
748            DO ji = fs_2, fs_jpim1
749               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  &
750                  &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
751            END DO
752         END DO
753         !
754      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
755         DO jj = 1, jpjm1
756            DO ji = fs_2, fs_jpim1
757               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
758               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (        ( pt(ji,jj+1) + pt(ji,jj) )  &
759                  &                                     - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) )
760            END DO
761         END DO
762         CALL lbc_lnk( pt_v, 'V',  1. )
763         !
764      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
765         DO jj = 1, jpjm1
766            DO ji = fs_2, fs_jpim1
767               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
768               zdy2 = e2v(ji,jj) * e2v(ji,jj)
769!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
770               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)       &
771                  &                                     -                        zcv   * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   &
772                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1) + ztv2(ji,jj)       &
773                  &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
774            END DO
775         END DO
776         !
777      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
778         DO jj = 1, jpjm1
779            DO ji = fs_2, fs_jpim1
780               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
781               zdy2 = e2v(ji,jj) * e2v(ji,jj)
782!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
783               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                           ( pt  (ji,jj+1) + pt  (ji,jj)       &
784                  &                                               -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   &
785                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                   ztv2(ji,jj+1) + ztv2(ji,jj)       &
786                  &                                               - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
787            END DO
788         END DO
789         !
790      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
791         DO jj = 1, jpjm1
792            DO ji = fs_2, fs_jpim1
793               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
794               zdy2 = e2v(ji,jj) * e2v(ji,jj)
795!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
796               zdy4 = zdy2 * zdy2
797               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)      &
798                  &                                                     -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )  &
799                  &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) *     (                   ztv2(ji,jj+1) + ztv2(ji,jj)      &
800                  &                                                     - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) )  &
801                  &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1) + ztv4(ji,jj)      &
802                  &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1) - ztv4(ji,jj) ) ) )
803            END DO
804         END DO
805         !
806      END SELECT
807      !                                                     !-- High order flux in j-direction  --!
808      DO jj = 1, jpjm1
809         DO ji = 1, fs_jpim1   ! vector opt.
810            pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj)
811         END DO
812      END DO
813      !
814   END SUBROUTINE ultimate_y
815     
816
817   SUBROUTINE nonosc_2d( pdt, ptc, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho )
818      !!---------------------------------------------------------------------
819      !!                    ***  ROUTINE nonosc  ***
820      !!     
821       !! **  Purpose :   compute monotonic tracer fluxes from the upstream
822      !!       scheme and the before field by a nonoscillatory algorithm
823      !!
824      !! **  Method  :   ... ???
825      !!       warning : ptc and pt_low must be masked, but the boundaries
826      !!       conditions on the fluxes are not necessary zalezak (1979)
827      !!       drange (1995) multi-dimensional forward-in-time and upstream-
828      !!       in-space based differencing for fluid
829      !!----------------------------------------------------------------------
830      REAL(wp)                     , INTENT(in   ) ::   pdt              ! tracer time-step
831      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt_low      ! before field & upstream guess of after field
832      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pfv_low, pfu_low ! upstream flux
833      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux
834      !
835      INTEGER  ::   ji, jj    ! dummy loop indices
836      REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt   ! local scalars
837      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign        !   -      -
838      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zdiv, zti_low, ztj_low
839      !!----------------------------------------------------------------------
840      zbig = 1.e+40_wp
841      zsml = 1.e-15_wp
842
843      ! test on divergence
844      DO jj = 2, jpjm1
845         DO ji = fs_2, fs_jpim1   ! vector opt. 
846            zdiv(ji,jj) =  - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) + pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) 
847         END DO
848      END DO
849      CALL lbc_lnk( zdiv, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign)
850
851      ! antidiffusive flux : high order minus low order
852      ! --------------------------------------------------
853      DO jj = 1, jpjm1
854         DO ji = 1, fs_jpim1   ! vector opt.
855            pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_low(ji,jj)
856            pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_low(ji,jj)
857         END DO
858      END DO
859
860      ! extreme case where pfu_ho has to be zero
861      ! ----------------------------------------
862      !                                    pfu_ho
863      !                           *         --->
864      !                        |      |  *   |        |
865      !                        |      |      |    *   |   
866      !                        |      |      |        |    *
867      !            t_low :       i-1     i       i+1       i+2   
868      IF( ll_prelimiter ) THEN
869
870         DO jj = 2, jpjm1
871            DO ji = fs_2, fs_jpim1 
872               zti_low(ji,jj)= pt_low(ji+1,jj  )
873               ztj_low(ji,jj)= pt_low(ji  ,jj+1)
874            END DO
875         END DO
876         CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. )
877
878         IF( ll_prelimiter_zalesak ) THEN
879            DO jj = 2, jpjm1
880               DO ji = fs_2, fs_jpim1
881                  IF(    SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji+1,jj) - pt_low (ji  ,jj) ) ) THEN
882                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj) - zti_low(ji  ,jj) ) )   pfu_ho(ji,jj) = 0.
883                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji-1,jj) ) )   pfu_ho(ji,jj) = 0.
884                  ENDIF
885                  IF(    SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji,jj+1) - pt_low (ji,jj  ) ) ) THEN
886                     IF( SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1) - ztj_low(ji,jj  ) ) )   pfv_ho(ji,jj) = 0.
887                     IF( SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji,jj  ) - pt_low (ji,jj-1) ) )   pfv_ho(ji,jj) = 0.
888                  ENDIF
889               END DO
890            END DO
891            CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
892
893         ELSEIF( ll_prelimiter_devore ) THEN
894            z1_dt = 1._wp / pdt
895            DO jj = 2, jpjm1
896               DO ji = fs_2, fs_jpim1
897                  zsign = SIGN( 1., pt_low(ji+1,jj) - pt_low(ji,jj) )
898                  pfu_ho(ji,jj) =  zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj)) , &
899                     &                          zsign * ( pt_low (ji  ,jj) - pt_low (ji-1,jj) ) * e1e2t(ji  ,jj) * z1_dt , &
900                     &                          zsign * ( zti_low(ji+1,jj) - zti_low(ji  ,jj) ) * e1e2t(ji+1,jj) * z1_dt ) )
901                 
902                  zsign = SIGN( 1., pt_low(ji,jj+1) - pt_low(ji,jj) )
903                  pfv_ho(ji,jj) =  zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj)) , &
904                     &                          zsign * ( pt_low (ji,jj  ) - pt_low (ji,jj-1) ) * e1e2t(ji,jj  ) * z1_dt , &
905                     &                          zsign * ( ztj_low(ji,jj+1) - ztj_low(ji,jj  ) ) * e1e2t(ji,jj+1) * z1_dt ) )
906               END DO
907            END DO
908            CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
909           
910         ENDIF
911         
912      ENDIF
913     
914      ! Search local extrema
915      ! --------------------
916      ! max/min of ptc & pt_low with large negative/positive value (-/+zbig) outside ice cover
917      DO jj = 1, jpj
918         DO ji = 1, jpi
919            IF( ptc(ji,jj) == 0._wp .AND. pt_low(ji,jj) == 0._wp .AND. zdiv(ji,jj) == 0._wp ) THEN
920               zbup(ji,jj) = -zbig
921               zbdo(ji,jj) =  zbig
922            ELSE               
923               zbup(ji,jj) = MAX( ptc(ji,jj) , pt_low(ji,jj) )
924               zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_low(ji,jj) )
925            ENDIF
926         END DO
927      END DO
928
929      z1_dt = 1._wp / pdt
930      DO jj = 2, jpjm1
931         DO ji = fs_2, fs_jpim1   ! vector opt.
932            !
933            !zup  = MAX( zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ), zbup(ji  ,jj-1), zbup(ji  ,jj+1) )  ! search max/min in neighbourhood
934            !zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1) )
935            !
936            zup  = MAX( zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ), zbup(ji  ,jj-1), zbup(ji  ,jj+1), &  ! search max/min in neighbourhood
937               &                     zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1)  )
938            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1), &
939               &                     zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1)  )
940            !
941            zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji  ,jj) ) &  ! positive/negative part of the flux
942               & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj  ) )
943            zneg = MAX( 0., pfu_ho(ji  ,jj) ) - MIN( 0., pfu_ho(ji-1,jj) ) &
944               & + MAX( 0., pfv_ho(ji,jj  ) ) - MIN( 0., pfv_ho(ji,jj-1) )
945            !
946            !                                  ! up & down beta terms
947!!clem            zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt
948!!clem            zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt
949            IF( zpos >= epsi20 ) THEN ; zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt
950            ELSE                      ; zbetup(ji,jj) = zbig
951            ENDIF
952            !
953            IF( zneg >= epsi20 ) THEN ; zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
954            ELSE                      ; zbetdo(ji,jj) = zbig
955            ENDIF
956            !
957         END DO
958      END DO
959      CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
960
961      ! monotonic flux in the y direction
962      ! ---------------------------------
963      DO jj = 1, jpjm1
964         DO ji = 1, fs_jpim1   ! vector opt.
965            zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) )
966            zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) )
967            zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj) )
968            !
969            pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_low(ji,jj)
970         END DO
971      END DO
972
973      DO jj = 1, jpjm1
974         DO ji = 1, fs_jpim1   ! vector opt.
975            zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) )
976            zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) )
977            zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj) )
978            !
979            pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_low(ji,jj)
980         END DO
981      END DO
982      !
983   END SUBROUTINE nonosc_2d
984
985   SUBROUTINE limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
986      !!---------------------------------------------------------------------
987      !!                    ***  ROUTINE limiter_x  ***
988      !!     
989      !! **  Purpose :   compute flux limiter
990      !!----------------------------------------------------------------------
991      REAL(wp)                     , INTENT(in   )           ::   pdt          ! tracer time-step
992      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pu           ! ice i-velocity => u*e2
993      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   puc          ! ice i-velocity *A => u*e2*a
994      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pt           ! ice tracer
995      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout)           ::   pfu_ho       ! high order flux
996      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ), OPTIONAL ::   pfu_ups      ! upstream flux
997      !
998      REAL(wp) ::   Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr
999      INTEGER  ::   ji, jj    ! dummy loop indices
1000      REAL(wp), DIMENSION (jpi,jpj) ::   zslpx       ! tracer slopes
1001      !!----------------------------------------------------------------------
1002      !
1003      DO jj = 2, jpjm1
1004         DO ji = fs_2, fs_jpim1   ! vector opt.
1005            zslpx(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * umask(ji,jj,1)
1006         END DO
1007      END DO
1008      CALL lbc_lnk( zslpx, 'U', -1.)   ! lateral boundary cond.
1009     
1010      DO jj = 2, jpjm1
1011         DO ji = fs_2, fs_jpim1   ! vector opt.
1012            uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
1013
1014            Rjm = zslpx(ji-1,jj)
1015            Rj  = zslpx(ji  ,jj)
1016            Rjp = zslpx(ji+1,jj)
1017
1018            IF( PRESENT(pfu_ups) ) THEN
1019
1020               IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm
1021               ELSE                        ;   Rr = Rjp
1022               ENDIF
1023
1024               zh3 = pfu_ho(ji,jj) - pfu_ups(ji,jj)     
1025               IF( Rj > 0. ) THEN
1026                  zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(puc(ji,jj)),  &
1027                     &        MIN( 2. * Rr * 0.5 * ABS(puc(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) )
1028               ELSE
1029                  zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(puc(ji,jj)),  &
1030                     &        MIN(-2. * Rr * 0.5 * ABS(puc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) )
1031               ENDIF
1032               pfu_ho(ji,jj) = pfu_ups(ji,jj) + zlimiter
1033               
1034            ELSE
1035               IF( Rj /= 0. ) THEN
1036                  IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
1037                  ELSE                        ;   Cr = Rjp / Rj
1038                  ENDIF
1039               ELSE
1040                  Cr = 0.
1041                  !IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20
1042                  !ELSE                        ;   Cr = Rjp * 1.e20
1043                  !ENDIF
1044               ENDIF
1045               
1046               ! -- superbee --
1047               zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1048               ! -- van albada 2 --
1049               !!zpsi = 2.*Cr / (Cr*Cr+1.)
1050
1051               ! -- sweby (with beta=1) --
1052               !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1053               ! -- van Leer --
1054               !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1055               ! -- ospre --
1056               !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
1057               ! -- koren --
1058               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
1059               ! -- charm --
1060               !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
1061               !ELSE                 ;   zpsi = 0.
1062               !ENDIF
1063               ! -- van albada 1 --
1064               !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
1065               ! -- smart --
1066               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
1067               ! -- umist --
1068               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1069
1070               ! high order flux corrected by the limiter
1071               pfu_ho(ji,jj) = pfu_ho(ji,jj) - ABS( puc(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5
1072
1073            ENDIF
1074         END DO
1075      END DO
1076      CALL lbc_lnk( pfu_ho, 'U', -1.)   ! lateral boundary cond.
1077      !
1078   END SUBROUTINE limiter_x
1079
1080   SUBROUTINE limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
1081      !!---------------------------------------------------------------------
1082      !!                    ***  ROUTINE limiter_y  ***
1083      !!     
1084      !! **  Purpose :   compute flux limiter
1085      !!----------------------------------------------------------------------
1086      REAL(wp)                     , INTENT(in   )           ::   pdt          ! tracer time-step
1087      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pv           ! ice i-velocity => u*e2
1088      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pvc          ! ice i-velocity *A => u*e2*a
1089      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pt           ! ice tracer
1090      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout)           ::   pfv_ho       ! high order flux
1091      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ), OPTIONAL ::   pfv_ups      ! upstream flux
1092      !
1093      REAL(wp) ::   Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr
1094      INTEGER  ::   ji, jj    ! dummy loop indices
1095      REAL(wp), DIMENSION (jpi,jpj) ::   zslpy       ! tracer slopes
1096      !!----------------------------------------------------------------------
1097      !
1098      DO jj = 2, jpjm1
1099         DO ji = fs_2, fs_jpim1   ! vector opt.
1100            zslpy(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * vmask(ji,jj,1)
1101         END DO
1102      END DO
1103      CALL lbc_lnk( zslpy, 'V', -1.)   ! lateral boundary cond.
1104     
1105      DO jj = 2, jpjm1
1106         DO ji = fs_2, fs_jpim1   ! vector opt.
1107            vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
1108
1109            Rjm = zslpy(ji,jj-1)
1110            Rj  = zslpy(ji,jj  )
1111            Rjp = zslpy(ji,jj+1)
1112
1113            IF( PRESENT(pfv_ups) ) THEN
1114
1115               IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm
1116               ELSE                        ;   Rr = Rjp
1117               ENDIF
1118
1119               zh3 = pfv_ho(ji,jj) - pfv_ups(ji,jj)     
1120               IF( Rj > 0. ) THEN
1121                  zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pvc(ji,jj)),  &
1122                     &        MIN( 2. * Rr * 0.5 * ABS(pvc(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) )
1123               ELSE
1124                  zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pvc(ji,jj)),  &
1125                     &        MIN(-2. * Rr * 0.5 * ABS(pvc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) )
1126               ENDIF
1127               pfv_ho(ji,jj) = pfv_ups(ji,jj) + zlimiter
1128
1129            ELSE
1130
1131               IF( Rj /= 0. ) THEN
1132                  IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
1133                  ELSE                        ;   Cr = Rjp / Rj
1134                  ENDIF
1135               ELSE
1136                  Cr = 0.
1137                  !IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20
1138                  !ELSE                        ;   Cr = Rjp * 1.e20
1139                  !ENDIF
1140               ENDIF
1141
1142               ! -- superbee --
1143               zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1144               ! -- van albada 2 --
1145               !!zpsi = 2.*Cr / (Cr*Cr+1.)
1146
1147               ! -- sweby (with beta=1) --
1148               !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1149               ! -- van Leer --
1150               !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1151               ! -- ospre --
1152               !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
1153               ! -- koren --
1154               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
1155               ! -- charm --
1156               !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
1157               !ELSE                 ;   zpsi = 0.
1158               !ENDIF
1159               ! -- van albada 1 --
1160               !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
1161               ! -- smart --
1162               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
1163               ! -- umist --
1164               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1165
1166               ! high order flux corrected by the limiter
1167               pfv_ho(ji,jj) = pfv_ho(ji,jj) - ABS( pvc(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5
1168
1169            ENDIF
1170         END DO
1171      END DO
1172      CALL lbc_lnk( pfv_ho, 'V', -1.)   ! lateral boundary cond.
1173      !
1174   END SUBROUTINE limiter_y
1175
1176#else
1177   !!----------------------------------------------------------------------
1178   !!   Default option           Dummy module         NO SI3 sea-ice model
1179   !!----------------------------------------------------------------------
1180#endif
1181
1182   !!======================================================================
1183END MODULE icedyn_adv_umx
Note: See TracBrowser for help on using the repository browser.