source: NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90 @ 10315

Last change on this file since 10315 was 10315, checked in by clem, 2 years ago

still playing around with nonosc limiter to get a good ice thickness at the ice edge. There is still a problem when the ice concentration is very small

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