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 @ 10475

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

some light cleaning

  • Property svn:keywords set to Id
File size: 67.7 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_aip(:,:,:)
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, 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, pubox, pvbox, 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._wp ) 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, 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   ) ::   pt_ups           ! upstream guess of tracer
465      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes
466      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes
467      !
468      INTEGER  ::   ji, jj, jl    ! dummy loop indices
469      REAL(wp) ::   ztra          ! local scalar
470      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zpt
471      !!----------------------------------------------------------------------
472      !
473      IF( .NOT.ll_hoxy ) THEN           !** no alternate directions **!
474         !
475         DO jl = 1, jpl
476            DO jj = 1, jpjm1
477               DO ji = 1, fs_jpim1
478                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) )
479                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) )
480               END DO
481            END DO
482         END DO
483         !
484         IF    ( kn_limiter == 1 ) THEN
485            CALL nonosc( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
486         ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN
487            CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
488            CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
489         ENDIF
490         !
491      ELSE                              !** alternate directions **!
492         !
493         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
494            !
495            DO jl = 1, jpl              !-- flux in x-direction
496               DO jj = 1, jpjm1
497                  DO ji = 1, fs_jpim1
498                     pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
499                  END DO
500               END DO
501            END DO
502            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
503
504            DO jl = 1, jpl              !-- first guess of tracer from u-flux
505               DO jj = 2, jpjm1
506                  DO ji = fs_2, fs_jpim1
507                     ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) )              &
508                        &   + ( pu    (ji,jj   ) - pu    (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk)
509                     !
510                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
511                  END DO
512               END DO
513            END DO
514            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
515
516            DO jl = 1, jpl              !-- flux in y-direction
517               DO jj = 1, jpjm1
518                  DO ji = 1, fs_jpim1
519                     pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) )
520                  END DO
521               END DO
522            END DO
523            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
524
525         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
526            !
527            DO jl = 1, jpl              !-- flux in y-direction
528               DO jj = 1, jpjm1
529                  DO ji = 1, fs_jpim1
530                     pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
531                  END DO
532               END DO
533            END DO
534            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
535            !
536            DO jl = 1, jpl              !-- first guess of tracer from v-flux
537               DO jj = 2, jpjm1
538                  DO ji = fs_2, fs_jpim1
539                     ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) )  &
540                        &   + ( pv    (ji,jj   ) - pv    (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
541                     !
542                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
543                  END DO
544               END DO
545            END DO
546            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
547            !
548            DO jl = 1, jpl              !-- flux in x-direction
549               DO jj = 1, jpjm1
550                  DO ji = 1, fs_jpim1
551                     pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) )
552                  END DO
553               END DO
554            END DO
555            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
556
557         ENDIF
558         IF( kn_limiter == 1 )   CALL nonosc( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
559         
560      ENDIF
561   
562   END SUBROUTINE cen2
563
564   
565   SUBROUTINE macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, pubox, pvbox, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
566      !!---------------------------------------------------------------------
567      !!                    ***  ROUTINE macho  ***
568      !!     
569      !! **  Purpose :   compute the high order fluxes using Ultimate-Macho scheme 
570      !!
571      !! **  Method  :   ...
572      !!
573      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
574      !!----------------------------------------------------------------------
575      REAL(wp)                        , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
576      INTEGER                         , INTENT(in   ) ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2)
577      INTEGER                         , INTENT(in   ) ::   jt               ! number of sub-iteration
578      INTEGER                         , INTENT(in   ) ::   kt               ! number of iteration
579      REAL(wp)                        , INTENT(in   ) ::   pdt              ! tracer time-step
580      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt               ! tracer fields
581      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
582      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pubox, pvbox     ! upstream velocity
583      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt_ups           ! upstream guess of tracer
584      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes
585      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes
586      !
587      INTEGER  ::   ji, jj, jl    ! dummy loop indices
588      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zt_u, zt_v, zpt
589      !!----------------------------------------------------------------------
590      !
591      IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
592         !
593         !                                                        !--  ultimate interpolation of pt at u-point  --!
594         CALL ultimate_x( kn_umx, pdt, pt, pu, zt_u, pfu_ho )
595         !                                                        !--  limiter in x --!
596         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
597         !                                                        !--  advective form update in zpt  --!
598         DO jl = 1, jpl
599            DO jj = 2, jpjm1
600               DO ji = fs_2, fs_jpim1
601                  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) &
602                     &                              + pt   (ji,jj,jl) * ( pu  (ji,jj   ) - pu  (ji-1,jj   ) ) * r1_e1e2t(ji,jj) &
603                     &                                                                                        * pamsk           &
604                     &                             ) * pdt ) * tmask(ji,jj,1) 
605               END DO
606            END DO
607         END DO
608         CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
609         !
610         !                                                        !--  ultimate interpolation of pt at v-point  --!
611         IF( ll_hoxy ) THEN
612            CALL ultimate_y( kn_umx, pdt, zpt, pv, zt_v, pfv_ho )
613         ELSE
614            CALL ultimate_y( kn_umx, pdt, pt , pv, zt_v, pfv_ho )
615         ENDIF
616         !                                                        !--  limiter in y --!
617         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
618         !         
619         !
620      ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
621         !
622         !                                                        !--  ultimate interpolation of pt at v-point  --!
623         CALL ultimate_y( kn_umx, pdt, pt, pv, zt_v, pfv_ho )
624         !                                                        !--  limiter in y --!
625         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
626         !                                                        !--  advective form update in zpt  --!
627         DO jl = 1, jpl
628            DO jj = 2, jpjm1
629               DO ji = fs_2, fs_jpim1
630                  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) &
631                     &                              + pt   (ji,jj,jl) * ( pv  (ji,jj   ) - pv  (ji,jj-1   ) ) * r1_e1e2t(ji,jj) &
632                     &                                                                                        * pamsk           &
633                     &                             ) * pdt ) * tmask(ji,jj,1) 
634               END DO
635            END DO
636         END DO
637         CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
638         !
639         !                                                        !--  ultimate interpolation of pt at u-point  --!
640         IF( ll_hoxy ) THEN
641            CALL ultimate_x( kn_umx, pdt, zpt, pu, zt_u, pfu_ho )
642         ELSE
643            CALL ultimate_x( kn_umx, pdt, pt , pu, zt_u, pfu_ho )
644         ENDIF
645         !                                                        !--  limiter in x --!
646         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
647         !
648      ENDIF
649
650      IF( kn_limiter == 1 )   CALL nonosc( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
651      !
652   END SUBROUTINE macho
653
654
655   SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho )
656      !!---------------------------------------------------------------------
657      !!                    ***  ROUTINE ultimate_x  ***
658      !!     
659      !! **  Purpose :   compute tracer at u-points
660      !!
661      !! **  Method  :   ...
662      !!
663      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
664      !!----------------------------------------------------------------------
665      INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
666      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step
667      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pu        ! ice i-velocity component
668      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt        ! tracer fields
669      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_u      ! tracer at u-point
670      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho    ! high order flux
671      !
672      INTEGER  ::   ji, jj, jl             ! dummy loop indices
673      REAL(wp) ::   zcu, zdx2, zdx4        !   -      -
674      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztu1, ztu2, ztu3, ztu4
675      !!----------------------------------------------------------------------
676      !
677      !                                                     !--  Laplacian in i-direction  --!
678      DO jl = 1, jpl
679         DO jj = 2, jpjm1         ! First derivative (gradient)
680            DO ji = 1, fs_jpim1
681               ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
682            END DO
683            !                     ! Second derivative (Laplacian)
684            DO ji = fs_2, fs_jpim1
685               ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj)
686            END DO
687         END DO
688      END DO
689      CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1. )
690      !
691      !                                                     !--  BiLaplacian in i-direction  --!
692      DO jl = 1, jpl
693         DO jj = 2, jpjm1         ! Third derivative
694            DO ji = 1, fs_jpim1
695               ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
696            END DO
697            !                     ! Fourth derivative
698            DO ji = fs_2, fs_jpim1
699               ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj)
700            END DO
701         END DO
702      END DO
703      CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1. )
704      !
705      !
706      SELECT CASE (kn_umx )
707      !
708      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
709         !       
710         DO jl = 1, jpl
711            DO jj = 1, jpjm1
712               DO ji = 1, fs_jpim1   ! vector opt.
713                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
714                     &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
715               END DO
716            END DO
717         END DO
718         !
719      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
720         !
721         DO jl = 1, jpl
722            DO jj = 1, jpjm1
723               DO ji = 1, fs_jpim1   ! vector opt.
724                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
725                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
726                     &                                                            - zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
727               END DO
728            END DO
729         END DO
730         
731      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
732         !
733         DO jl = 1, jpl
734            DO jj = 1, jpjm1
735               DO ji = 1, fs_jpim1   ! vector opt.
736                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
737                  zdx2 = e1u(ji,jj) * e1u(ji,jj)
738!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj)
739                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
740                     &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
741                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
742                     &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
743               END DO
744            END DO
745         END DO
746         !
747      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
748         !
749         DO jl = 1, jpl
750            DO jj = 1, jpjm1
751               DO ji = 1, fs_jpim1   ! vector opt.
752                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
753                  zdx2 = e1u(ji,jj) * e1u(ji,jj)
754!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj)
755                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
756                     &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
757                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
758                     &                                                   - 0.5_wp * zcu   * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
759               END DO
760            END DO
761         END DO
762         !
763      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
764         !
765         DO jl = 1, jpl
766            DO jj = 1, jpjm1
767               DO ji = 1, fs_jpim1   ! vector opt.
768                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
769                  zdx2 = e1u(ji,jj) * e1u(ji,jj)
770!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj)
771                  zdx4 = zdx2 * zdx2
772                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (        (                       pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
773                     &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
774                     &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) * (                       ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
775                     &                                                   - 0.5_wp * zcu   * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) &
776                     &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)     &
777                     &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
778               END DO
779            END DO
780         END DO
781         !
782      END SELECT
783      !
784      ! if pt at u-point is negative then use the upstream value
785      !    this should not be necessary if a proper sea-ice mask is set in Ultimate
786      !    to degrade the order of the scheme when necessary (for ex. at the ice edge)
787      IF( ll_neg ) THEN
788         DO jl = 1, jpl
789            DO jj = 1, jpjm1
790               DO ji = 1, fs_jpim1
791                  IF( pt_u(ji,jj,jl) < 0._wp ) THEN
792                     pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
793                        &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
794                  ENDIF
795               END DO
796            END DO
797         END DO
798      ENDIF
799      !                                                     !-- High order flux in i-direction  --!
800      DO jl = 1, jpl
801         DO jj = 1, jpjm1
802            DO ji = 1, fs_jpim1   ! vector opt.
803               pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl)
804            END DO
805         END DO
806      END DO
807      !
808   END SUBROUTINE ultimate_x
809   
810 
811   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho )
812      !!---------------------------------------------------------------------
813      !!                    ***  ROUTINE ultimate_y  ***
814      !!     
815      !! **  Purpose :   compute tracer at v-points
816      !!
817      !! **  Method  :   ...
818      !!
819      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
820      !!----------------------------------------------------------------------
821      INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
822      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step
823      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pv        ! ice j-velocity component
824      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt        ! tracer fields
825      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_v      ! tracer at v-point
826      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfv_ho    ! high order flux
827      !
828      INTEGER  ::   ji, jj, jl         ! dummy loop indices
829      REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
830      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztv1, ztv2, ztv3, ztv4
831      !!----------------------------------------------------------------------
832      !
833      !                                                     !--  Laplacian in j-direction  --!
834      DO jl = 1, jpl
835         DO jj = 1, jpjm1         ! First derivative (gradient)
836            DO ji = fs_2, fs_jpim1
837               ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
838            END DO
839         END DO
840         DO jj = 2, jpjm1         ! Second derivative (Laplacian)
841            DO ji = fs_2, fs_jpim1
842               ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj)
843            END DO
844         END DO
845      END DO
846      CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1. )
847      !
848      !                                                     !--  BiLaplacian in j-direction  --!
849      DO jl = 1, jpl
850         DO jj = 1, jpjm1         ! First derivative
851            DO ji = fs_2, fs_jpim1
852               ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
853            END DO
854         END DO
855         DO jj = 2, jpjm1         ! Second derivative
856            DO ji = fs_2, fs_jpim1
857               ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj)
858            END DO
859         END DO
860      END DO
861      CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1. )
862      !
863      !
864      SELECT CASE (kn_umx )
865         !
866      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
867         DO jl = 1, jpl
868            DO jj = 1, jpjm1
869               DO ji = 1, fs_jpim1
870                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   &
871                     &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
872               END DO
873            END DO
874         END DO
875         !
876      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
877         DO jl = 1, jpl
878            DO jj = 1, jpjm1
879               DO ji = 1, fs_jpim1
880                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
881                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   &
882                     &                                                            - zcv *   ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
883               END DO
884            END DO
885         END DO
886         !
887      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
888         DO jl = 1, jpl
889            DO jj = 1, jpjm1
890               DO ji = 1, fs_jpim1
891                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
892                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
893!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj)
894                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
895                     &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
896                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
897                     &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
898               END DO
899            END DO
900         END DO
901         !
902      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
903         DO jl = 1, jpl
904            DO jj = 1, jpjm1
905               DO ji = 1, fs_jpim1
906                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
907                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
908!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj)
909                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
910                     &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
911                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
912                     &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
913               END DO
914            END DO
915         END DO
916         !
917      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
918         DO jl = 1, jpl
919            DO jj = 1, jpjm1
920               DO ji = 1, fs_jpim1
921                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
922                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
923!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj)
924                  zdy4 = zdy2 * zdy2
925                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
926                     &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
927                     &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) * (                       ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
928                     &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) &
929                     &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl)     &
930                     &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
931               END DO
932            END DO
933         END DO
934         !
935      END SELECT
936      !
937      ! if pt at v-point is negative then use the upstream value
938      !    this should not be necessary if a proper sea-ice mask is set in Ultimate
939      !    to degrade the order of the scheme when necessary (for ex. at the ice edge)
940      IF( ll_neg ) THEN
941         DO jl = 1, jpl
942            DO jj = 1, jpjm1
943               DO ji = 1, fs_jpim1
944                  IF( pt_v(ji,jj,jl) < 0._wp ) THEN
945                     pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
946                        &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
947                  ENDIF
948               END DO
949            END DO
950         END DO
951      ENDIF
952      !                                                     !-- High order flux in j-direction  --!
953      DO jl = 1, jpl
954         DO jj = 1, jpjm1
955            DO ji = 1, fs_jpim1   ! vector opt.
956               pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
957            END DO
958         END DO
959      END DO
960      !
961   END SUBROUTINE ultimate_y
962     
963
964   SUBROUTINE nonosc( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
965      !!---------------------------------------------------------------------
966      !!                    ***  ROUTINE nonosc  ***
967      !!     
968      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
969      !!       scheme and the before field by a nonoscillatory algorithm
970      !!
971      !! **  Method  :   ...
972      !!----------------------------------------------------------------------
973      REAL(wp)                   , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
974      REAL(wp)                   , INTENT(in   ) ::   pdt              ! tracer time-step
975      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pu               ! ice i-velocity => u*e2
976      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pv               ! ice j-velocity => v*e1
977      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pt, pt_ups       ! before field & upstream guess of after field
978      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pfv_ups, pfu_ups ! upstream flux
979      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux
980      !
981      INTEGER  ::   ji, jj, jl    ! dummy loop indices
982      REAL(wp) ::   zpos, zneg, zbig, zup, zdo, z1_dt              ! local scalars
983      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zcoef, zzt       !   -      -
984      REAL(wp), DIMENSION(jpi,jpj    ) :: zbup, zbdo
985      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_ups, ztj_ups
986      !!----------------------------------------------------------------------
987      zbig = 1.e+40_wp
988     
989      ! antidiffusive flux : high order minus low order
990      ! --------------------------------------------------
991      DO jl = 1, jpl
992         DO jj = 1, jpjm1
993            DO ji = 1, fs_jpim1   ! vector opt.
994               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
995               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
996            END DO
997         END DO
998      END DO
999
1000      ! extreme case where pfu_ho has to be zero
1001      ! ----------------------------------------
1002      !                                    pfu_ho
1003      !                           *         --->
1004      !                        |      |  *   |        |
1005      !                        |      |      |    *   |   
1006      !                        |      |      |        |    *
1007      !            t_ups :       i-1     i       i+1       i+2   
1008      IF( ll_prelimiter_zalesak ) THEN
1009         
1010         DO jl = 1, jpl
1011            DO jj = 2, jpjm1
1012               DO ji = fs_2, fs_jpim1 
1013                  zti_ups(ji,jj,jl)= pt_ups(ji+1,jj  ,jl)
1014                  ztj_ups(ji,jj,jl)= pt_ups(ji  ,jj+1,jl)
1015               END DO
1016            END DO
1017         END DO
1018         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. )
1019
1020         DO jl = 1, jpl
1021            DO jj = 2, jpjm1
1022               DO ji = fs_2, fs_jpim1
1023                  IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj  ,jl) - pt_ups(ji,jj,jl) ) <= 0._wp .AND.  &
1024                     & pfv_ho(ji,jj,jl) * ( pt_ups(ji  ,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0._wp ) THEN
1025                     !
1026                     IF(  pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj  ,jl) - zti_ups(ji,jj,jl) ) <= 0._wp .AND.  &
1027                        & pfv_ho(ji,jj,jl) * ( ztj_ups(ji  ,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0._wp ) THEN
1028                        pfu_ho(ji,jj,jl)=0._wp
1029                        pfv_ho(ji,jj,jl)=0._wp
1030                     ENDIF
1031                     !
1032                     IF(  pfu_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji-1,jj  ,jl) ) <= 0._wp .AND.  &
1033                        & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji  ,jj-1,jl) ) <= 0._wp ) THEN
1034                        pfu_ho(ji,jj,jl)=0._wp
1035                        pfv_ho(ji,jj,jl)=0._wp
1036                     ENDIF
1037                     !
1038                  ENDIF
1039               END DO
1040            END DO
1041         END DO
1042         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
1043
1044      ENDIF
1045
1046      ! Search local extrema
1047      ! --------------------
1048      ! max/min of pt & pt_ups with large negative/positive value (-/+zbig) outside ice cover
1049      z1_dt = 1._wp / pdt
1050      DO jl = 1, jpl
1051         
1052         DO jj = 1, jpj
1053            DO ji = 1, jpi
1054               IF    ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
1055                  zbup(ji,jj) = -zbig
1056                  zbdo(ji,jj) =  zbig
1057               ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN
1058                  zbup(ji,jj) = pt_ups(ji,jj,jl)
1059                  zbdo(ji,jj) = pt_ups(ji,jj,jl)
1060               ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
1061                  zbup(ji,jj) = pt(ji,jj,jl)
1062                  zbdo(ji,jj) = pt(ji,jj,jl)
1063               ELSE
1064                  zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
1065                  zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
1066               ENDIF
1067            END DO
1068         END DO
1069
1070         DO jj = 2, jpjm1
1071            DO ji = fs_2, fs_jpim1   ! vector opt.
1072               !
1073               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
1074               zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) )
1075               !
1076               zpos = MAX( 0._wp, pfu_ho(ji-1,jj  ,jl) ) - MIN( 0._wp, pfu_ho(ji  ,jj  ,jl) ) &  ! positive/negative part of the flux
1077                  & + MAX( 0._wp, pfv_ho(ji  ,jj-1,jl) ) - MIN( 0._wp, pfv_ho(ji  ,jj  ,jl) )
1078               zneg = MAX( 0._wp, pfu_ho(ji  ,jj  ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj  ,jl) ) &
1079                  & + MAX( 0._wp, pfv_ho(ji  ,jj  ,jl) ) - MIN( 0._wp, pfv_ho(ji  ,jj-1,jl) )
1080               !
1081               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) ) &
1082                  &          ) * ( 1. - pamsk )
1083               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) ) &
1084                  &          ) * ( 1. - pamsk )
1085               !
1086               !                                  ! up & down beta terms
1087               IF( zpos > 0._wp ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt
1088               ELSE                    ; zbetup(ji,jj,jl) = 0._wp ! zbig
1089               ENDIF
1090               !
1091               IF( zneg > 0._wp ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
1092               ELSE                    ; zbetdo(ji,jj,jl) = 0._wp ! zbig
1093               ENDIF
1094               !
1095               ! if all the points are outside ice cover
1096               IF( zup == -zbig )   zbetup(ji,jj,jl) = 0._wp ! zbig
1097               IF( zdo ==  zbig )   zbetdo(ji,jj,jl) = 0._wp ! zbig           
1098               !
1099            END DO
1100         END DO
1101      END DO
1102      CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
1103
1104     
1105      ! monotonic flux in the y direction
1106      ! ---------------------------------
1107      DO jl = 1, jpl
1108         DO jj = 1, jpjm1
1109            DO ji = 1, fs_jpim1   ! vector opt.
1110               zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) )
1111               zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) )
1112               zcu = 0.5_wp + SIGN( 0.5_wp , pfu_ho(ji,jj,jl) )
1113               !
1114               zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
1115               !
1116               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl)
1117               !
1118            END DO
1119         END DO
1120
1121         DO jj = 1, jpjm1
1122            DO ji = 1, fs_jpim1   ! vector opt.
1123               zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) )
1124               zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) )
1125               zcv = 0.5_wp + SIGN( 0.5_wp , pfv_ho(ji,jj,jl) )
1126               !
1127               zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
1128               !
1129               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl)
1130               !
1131            END DO
1132         END DO
1133
1134         ! clem test
1135!!         DO jj = 2, jpjm1
1136!!            DO ji = 2, fs_jpim1   ! vector opt.
1137!!               zzt = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1138!!                  &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1139!!                  &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1140!!                  &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1141!!                  &         ) * tmask(ji,jj,1)
1142!!               IF( zzt < -epsi20 ) THEN
1143!!                  WRITE(numout,*) 'T<0 nonosc',zzt
1144!!               ENDIF
1145!!            END DO
1146!!         END DO
1147
1148      END DO
1149      !
1150   END SUBROUTINE nonosc
1151
1152   
1153   SUBROUTINE limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
1154      !!---------------------------------------------------------------------
1155      !!                    ***  ROUTINE limiter_x  ***
1156      !!     
1157      !! **  Purpose :   compute flux limiter
1158      !!----------------------------------------------------------------------
1159      REAL(wp)                  , INTENT(in   ) ::   pdt          ! tracer time-step
1160      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu           ! ice i-velocity => u*e2
1161      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt           ! ice tracer
1162      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pfu_ups      ! upstream flux
1163      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pfu_ho       ! high order flux
1164      !
1165      REAL(wp) ::   Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr
1166      INTEGER  ::   ji, jj, jl    ! dummy loop indices
1167      REAL(wp), DIMENSION (jpi,jpj,jpl) ::   zslpx       ! tracer slopes
1168      !!----------------------------------------------------------------------
1169      !
1170      DO jl = 1, jpl
1171         DO jj = 2, jpjm1
1172            DO ji = fs_2, fs_jpim1   ! vector opt.
1173               zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1)
1174            END DO
1175         END DO
1176      END DO
1177      CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.)   ! lateral boundary cond.
1178     
1179      DO jl = 1, jpl
1180         DO jj = 2, jpjm1
1181            DO ji = fs_2, fs_jpim1   ! vector opt.
1182               uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
1183               
1184               Rjm = zslpx(ji-1,jj,jl)
1185               Rj  = zslpx(ji  ,jj,jl)
1186               Rjp = zslpx(ji+1,jj,jl)
1187
1188               IF( kn_limiter == 3 ) THEN
1189
1190                  IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm
1191                  ELSE                        ;   Rr = Rjp
1192                  ENDIF
1193
1194                  zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)     
1195                  IF( Rj > 0. ) THEN
1196                     zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)),  &
1197                        &        MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
1198                  ELSE
1199                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)),  &
1200                        &        MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
1201                  ENDIF
1202                  pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter
1203
1204               ELSEIF( kn_limiter == 2 ) THEN
1205                  IF( Rj /= 0. ) THEN
1206                     IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
1207                     ELSE                        ;   Cr = Rjp / Rj
1208                     ENDIF
1209                  ELSE
1210                     Cr = 0.
1211                  ENDIF
1212
1213                  ! -- superbee --
1214                  zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1215                  ! -- van albada 2 --
1216                  !!zpsi = 2.*Cr / (Cr*Cr+1.)
1217                  ! -- sweby (with beta=1) --
1218                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1219                  ! -- van Leer --
1220                  !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1221                  ! -- ospre --
1222                  !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
1223                  ! -- koren --
1224                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
1225                  ! -- charm --
1226                  !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
1227                  !ELSE                 ;   zpsi = 0.
1228                  !ENDIF
1229                  ! -- van albada 1 --
1230                  !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
1231                  ! -- smart --
1232                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
1233                  ! -- umist --
1234                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1235
1236                  ! high order flux corrected by the limiter
1237                  pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5
1238
1239               ENDIF
1240            END DO
1241         END DO
1242      END DO
1243      CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.)   ! lateral boundary cond.
1244      !
1245   END SUBROUTINE limiter_x
1246
1247   
1248   SUBROUTINE limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
1249      !!---------------------------------------------------------------------
1250      !!                    ***  ROUTINE limiter_y  ***
1251      !!     
1252      !! **  Purpose :   compute flux limiter
1253      !!----------------------------------------------------------------------
1254      REAL(wp)                   , INTENT(in   ) ::   pdt          ! tracer time-step
1255      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pv           ! ice i-velocity => u*e2
1256      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pt           ! ice tracer
1257      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pfv_ups      ! upstream flux
1258      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ho       ! high order flux
1259      !
1260      REAL(wp) ::   Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr
1261      INTEGER  ::   ji, jj, jl    ! dummy loop indices
1262      REAL(wp), DIMENSION (jpi,jpj,jpl) ::   zslpy       ! tracer slopes
1263      !!----------------------------------------------------------------------
1264      !
1265      DO jl = 1, jpl
1266         DO jj = 2, jpjm1
1267            DO ji = fs_2, fs_jpim1   ! vector opt.
1268               zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1)
1269            END DO
1270         END DO
1271      END DO
1272      CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.)   ! lateral boundary cond.
1273
1274      DO jl = 1, jpl
1275         DO jj = 2, jpjm1
1276            DO ji = fs_2, fs_jpim1   ! vector opt.
1277               vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
1278
1279               Rjm = zslpy(ji,jj-1,jl)
1280               Rj  = zslpy(ji,jj  ,jl)
1281               Rjp = zslpy(ji,jj+1,jl)
1282
1283               IF( kn_limiter == 3 ) THEN
1284
1285                  IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm
1286                  ELSE                        ;   Rr = Rjp
1287                  ENDIF
1288
1289                  zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)     
1290                  IF( Rj > 0. ) THEN
1291                     zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)),  &
1292                        &        MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
1293                  ELSE
1294                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)),  &
1295                        &        MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
1296                  ENDIF
1297                  pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter
1298
1299               ELSEIF( kn_limiter == 2 ) THEN
1300
1301                  IF( Rj /= 0. ) THEN
1302                     IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
1303                     ELSE                        ;   Cr = Rjp / Rj
1304                     ENDIF
1305                  ELSE
1306                     Cr = 0.
1307                  ENDIF
1308
1309                  ! -- superbee --
1310                  zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1311                  ! -- van albada 2 --
1312                  !!zpsi = 2.*Cr / (Cr*Cr+1.)
1313                  ! -- sweby (with beta=1) --
1314                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1315                  ! -- van Leer --
1316                  !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1317                  ! -- ospre --
1318                  !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
1319                  ! -- koren --
1320                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
1321                  ! -- charm --
1322                  !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
1323                  !ELSE                 ;   zpsi = 0.
1324                  !ENDIF
1325                  ! -- van albada 1 --
1326                  !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
1327                  ! -- smart --
1328                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
1329                  ! -- umist --
1330                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1331
1332                  ! high order flux corrected by the limiter
1333                  pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5
1334
1335               ENDIF
1336            END DO
1337         END DO
1338      END DO
1339      CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.)   ! lateral boundary cond.
1340      !
1341   END SUBROUTINE limiter_y
1342
1343#else
1344   !!----------------------------------------------------------------------
1345   !!   Default option           Dummy module         NO SI3 sea-ice model
1346   !!----------------------------------------------------------------------
1347#endif
1348
1349   !!======================================================================
1350END MODULE icedyn_adv_umx
Note: See TracBrowser for help on using the repository browser.