source: NEMO/trunk/src/ICE/icedyn_adv_umx.F90 @ 10519

Last change on this file since 10519 was 10519, checked in by clem, 19 months ago

solve a compilation issue with agrif. It looks like agrif does not like 2 subroutines to have the same name eventhough they are not in the same module and they are not declared as public. I would like to say wtf but lets remain polite and say whats going on? Does agrif consider all the subroutines as public?

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