New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_adv_umx.F90 in NEMO/trunk/src/ICE – NEMO

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

Last change on this file since 10446 was 10446, checked in by clem, 4 years ago

clean the routine icedyn_adv_umx

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