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

Last change on this file since 10579 was 10579, checked in by clem, 21 months ago

change an option in the advection scheme to avoid weird ice temperature values at very low ice concentration in ORCA2. However this option breaks the test case ICE_ADV2D. I chose ORCA2 over ICE_ADV2D

  • 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_ice        : compute monotonic tracer fluxes by a non-oscillatory algorithm
17   !!----------------------------------------------------------------------
18   USE phycst         ! physical constant
19   USE dom_oce        ! ocean domain
20   USE sbc_oce , ONLY : nn_fsbc   ! update frequency of surface boundary condition
21   USE ice            ! sea-ice variables
22   USE icevar         ! sea-ice: operations
23   !
24   USE in_out_manager ! I/O manager
25   USE lib_mpp        ! MPP library
26   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
27   USE lbclnk         ! lateral boundary conditions (or mpp links)
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   ice_dyn_adv_umx   ! called by icedyn_adv.F90
33     
34   REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6
35   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
36   
37   ! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid)
38   INTEGER ::   kn_limiter = 1
39
40   ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme
41   !   clem: if set to true, the 2D test case "diagonal advection" does not work (I do not understand why)
42   !         but in realistic cases, it avoids having very negative ice temperature (-50) at low ice concentration
43   LOGICAL ::   ll_neg = .TRUE.
44   
45   ! alternate directions for upstream
46   LOGICAL ::   ll_upsxy = .TRUE.
47
48   ! alternate directions for high order
49   LOGICAL ::   ll_hoxy = .TRUE.
50   
51   ! prelimiter: use it to avoid overshoot in H
52   !   clem: if set to true, the 2D test case "diagnoal advection" does not work (I do not understand why)
53   LOGICAL ::   ll_prelimiter_zalesak = .FALSE.  ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D
54
55
56   !! * Substitutions
57#  include "vectopt_loop_substitute.h90"
58   !!----------------------------------------------------------------------
59   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
60   !! $Id$
61   !! Software governed by the CeCILL licence     (./LICENSE)
62   !!----------------------------------------------------------------------
63CONTAINS
64
65   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice,  &
66      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
67      !!----------------------------------------------------------------------
68      !!                  ***  ROUTINE ice_dyn_adv_umx  ***
69      !!
70      !! **  Purpose :   Compute the now trend due to total advection of
71      !!                 tracers and add it to the general trend of tracer equations
72      !!                 using an "Ultimate-Macho" scheme
73      !!
74      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
75      !!----------------------------------------------------------------------
76      INTEGER                     , INTENT(in   ) ::   kn_umx     ! order of the scheme (1-5=UM or 20=CEN2)
77      INTEGER                     , INTENT(in   ) ::   kt         ! time step
78      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
79      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
80      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
81      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
82      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
83      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
84      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
85      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
86      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
87      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
88      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
89      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
90      !
91      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
92      INTEGER  ::   icycle                  ! number of sub-timestep for the advection
93      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers
94      REAL(wp) ::   zdt
95      REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! send zcflnow and receive zcflprv
96      REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box 
97      REAL(wp), DIMENSION(jpi,jpj)     ::   zati1, zati2
98      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho
99      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip
100      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhvar
101      !!----------------------------------------------------------------------
102      !
103      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
104      !
105      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !
106      !     When needed, the advection split is applied at the next time-step in order to avoid blocking global comm.
107      !     ...this should not affect too much the stability... Was ok on the tests we did...
108      zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
109      zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
110     
111      ! non-blocking global communication send zcflnow and receive zcflprv
112      CALL mpp_delay_max( 'icedyn_adv_umx', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
113
114      IF( zcflprv(1) > .5 ) THEN   ;   icycle = 2
115      ELSE                         ;   icycle = 1
116      ENDIF
117     
118      zdt = rdt_ice / REAL(icycle)
119
120      ! --- transport --- !
121      zudy(:,:) = pu_ice(:,:) * e2u(:,:)
122      zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
123
124      ! --- define velocity for advection: u*grad(H) --- !
125      DO jj = 2, jpjm1
126         DO ji = fs_2, fs_jpim1
127            IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
128            ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj)
129            ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj)
130            ENDIF
131
132            IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
133            ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1)
134            ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  )
135            ENDIF
136         END DO
137      END DO
138
139      !---------------!
140      !== advection ==!
141      !---------------!
142      DO jt = 1, icycle
143
144         ! record at_i before advection (for open water)
145         zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
146         
147         ! inverse of A and Ap
148         WHERE( pa_i(:,:,:) >= epsi20 )   ;   z1_ai(:,:,:) = 1._wp / pa_i(:,:,:)
149         ELSEWHERE                        ;   z1_ai(:,:,:) = 0.
150         END WHERE
151         WHERE( pa_ip(:,:,:) >= epsi20 )  ;   z1_aip(:,:,:) = 1._wp / pa_ip(:,:,:)
152         ELSEWHERE                        ;   z1_aip(:,:,:) = 0.
153         END WHERE
154         !
155         ! set u*a=u for advection of A only
156         DO jl = 1, jpl
157            zua_ho(:,:,jl) = zudy(:,:)
158            zva_ho(:,:,jl) = zvdx(:,:)
159         END DO
160         
161         zamsk = 1._wp
162         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
163         zamsk = 0._wp
164         !
165         zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
166         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_i )                ! Ice volume
167         !
168         zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
169         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_s )                ! Snw volume
170         !
171         zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
172         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, psv_i )               ! Salt content
173         !
174         zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:)
175         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i )               ! Age content
176         !
177         DO jk = 1, nlay_i
178            zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
179            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
180         END DO
181         !
182         DO jk = 1, nlay_s
183            zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
184            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
185         END DO
186         !
187         IF ( ln_pnd_H12 ) THEN
188            ! set u*a=u for advection of Ap only
189            DO jl = 1, jpl
190               zua_ho(:,:,jl) = zudy(:,:)
191               zva_ho(:,:,jl) = zvdx(:,:)
192            END DO
193           
194            zamsk = 1._wp
195            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
196            zamsk = 0._wp
197            !
198            zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:)
199            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_ip )                 ! mp volume
200         ENDIF
201         !
202         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
203         DO jj = 2, jpjm1
204            DO ji = fs_2, fs_jpim1
205               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                                                   ! Open water area
206                  &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
207            END DO
208         END DO
209         CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T',  1. )
210         !
211      END DO
212      !
213   END SUBROUTINE ice_dyn_adv_umx
214
215   
216   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho )
217      !!----------------------------------------------------------------------
218      !!                  ***  ROUTINE adv_umx  ***
219      !!
220      !! **  Purpose :   Compute the now trend due to total advection of
221      !!                 tracers and add it to the general trend of tracer equations
222      !!
223      !! **  Method  :   - calculate upstream fluxes and upstream solution for tracer H
224      !!                 - calculate tracer H at u and v points (Ultimate)
225      !!                 - calculate the high order fluxes using alterning directions (Macho?)
226      !!                 - apply a limiter on the fluxes (nonosc_ice)
227      !!                 - convert this tracer flux to a tracer content flux (uH -> uV)
228      !!                 - calculate the high order solution for tracer content V
229      !!
230      !! ** Action : solve 2 equations => a) da/dt = -div(ua)
231      !!                                  b) dV/dt = -div(uV) using dH/dt = -u.grad(H)
232      !!             in eq. b), - fluxes uH are evaluated (with UMx) and limited (with nonosc_ice). This step is necessary to get a good H.
233      !!                        - then we convert this flux to a "volume" flux this way => uH*ua/u
234      !!                             where ua is the flux from eq. a)
235      !!                        - at last we estimate dV/dt = -div(uH*ua/u)
236      !!
237      !! ** Note : - this method can lead to small negative V (since we only limit H) => corrected in icedyn_adv.F90 conserving mass etc.
238      !!           - negative tracers at u-v points can also occur from the Ultimate scheme (usually at the ice edge) and the solution for now
239      !!             is to apply an upstream scheme when it occurs. A better solution would be to degrade the order of
240      !!             the scheme automatically by applying a mask of the ice cover inside Ultimate (not done).
241      !!           - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good.
242      !!             Large values of H can appear for very small ice concentration, and when it does it messes the things up since we
243      !!             work on H (and not V). It probably comes from the prelimiter of zalesak which is coded for 1D and not 2D.
244      !!             Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice
245      !!             concentration is small).
246      !! To-do: expand the prelimiter from zalesak to make it work in 2D
247      !!----------------------------------------------------------------------
248      REAL(wp)                        , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0)
249      INTEGER                         , INTENT(in   )           ::   kn_umx         ! order of the scheme (1-5=UM or 20=CEN2)
250      INTEGER                         , INTENT(in   )           ::   jt             ! number of sub-iteration
251      INTEGER                         , INTENT(in   )           ::   kt             ! number of iteration
252      REAL(wp)                        , INTENT(in   )           ::   pdt            ! tracer time-step
253      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2
254      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u
255      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pubox, pvbox   ! upstream velocity
256      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   pt             ! tracer field
257      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   ptc            ! tracer content field
258      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes
259      !
260      INTEGER  ::   ji, jj, jl       ! dummy loop indices 
261      REAL(wp) ::   ztra             ! local scalar
262      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ho , zfv_ho , zpt
263      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ups, zfv_ups, zt_ups
264      !!----------------------------------------------------------------------
265      !
266      ! Upstream (_ups) fluxes
267      ! -----------------------
268      CALL upstream( pamsk, jt, kt, pdt, pt, pu, pv, zt_ups, zfu_ups, zfv_ups )
269     
270      ! High order (_ho) fluxes
271      ! -----------------------
272      SELECT CASE( kn_umx )
273         !
274      CASE ( 20 )                          !== centered second order ==!
275         !
276         CALL cen2( pamsk, jt, kt, pdt, pt, pu, pv, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
277         !
278      CASE ( 1:5 )                         !== 1st to 5th order ULTIMATE-MACHO scheme ==!
279         !
280         CALL macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, pubox, pvbox, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
281         !
282      END SELECT
283      !
284      !              --ho    --ho
285      ! new fluxes = u*H  *  u*a / u
286      ! ----------------------------
287      IF( pamsk == 0._wp ) THEN
288         DO jl = 1, jpl
289            DO jj = 1, jpjm1
290               DO ji = 1, fs_jpim1
291                  IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN
292                     zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
293                     zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
294                  ELSE
295                     zfu_ho (ji,jj,jl) = 0._wp
296                     zfu_ups(ji,jj,jl) = 0._wp
297                  ENDIF
298                  !
299                  IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN
300                     zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
301                     zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
302                  ELSE
303                     zfv_ho (ji,jj,jl) = 0._wp 
304                     zfv_ups(ji,jj,jl) = 0._wp 
305                  ENDIF
306               END DO
307            END DO
308         END DO
309      ENDIF
310      !                                   --ho
311      ! in case of advection of A: output u*a
312      ! -------------------------------------
313      IF( PRESENT( pua_ho ) ) THEN
314         DO jl = 1, jpl
315            DO jj = 1, jpjm1
316               DO ji = 1, fs_jpim1
317                  pua_ho(ji,jj,jl) = zfu_ho(ji,jj,jl)
318                  pva_ho(ji,jj,jl) = zfv_ho(ji,jj,jl)
319               END DO
320            END DO
321         END DO
322      ENDIF
323      !
324      ! final trend with corrected fluxes
325      ! ---------------------------------
326      DO jl = 1, jpl
327         DO jj = 2, jpjm1
328            DO ji = fs_2, fs_jpim1 
329               ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) 
330               !
331               ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)               
332            END DO
333         END DO
334      END DO
335      CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1. )
336      !
337   END SUBROUTINE adv_umx
338
339
340   SUBROUTINE upstream( pamsk, jt, kt, pdt, pt, pu, pv, pt_ups, pfu_ups, pfv_ups )
341      !!---------------------------------------------------------------------
342      !!                    ***  ROUTINE upstream  ***
343      !!     
344      !! **  Purpose :   compute the upstream fluxes and upstream guess of tracer
345      !!----------------------------------------------------------------------
346      REAL(wp)                        , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
347      INTEGER                         , INTENT(in   ) ::   jt               ! number of sub-iteration
348      INTEGER                         , INTENT(in   ) ::   kt               ! number of iteration
349      REAL(wp)                        , INTENT(in   ) ::   pdt              ! tracer time-step
350      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt               ! tracer fields
351      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
352      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_ups           ! upstream guess of tracer
353      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ups, pfv_ups ! upstream fluxes
354      !
355      INTEGER  ::   ji, jj, jl    ! dummy loop indices
356      REAL(wp) ::   ztra          ! local scalar
357      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zpt
358      !!----------------------------------------------------------------------
359
360      IF( .NOT. ll_upsxy ) THEN         !** no alternate directions **!
361         !
362         DO jl = 1, jpl
363            DO jj = 1, jpjm1
364               DO ji = 1, fs_jpim1
365                  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)
366                  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)
367               END DO
368            END DO
369         END DO
370         !
371      ELSE                              !** alternate directions **!
372         !
373         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
374            !
375            DO jl = 1, jpl              !-- flux in x-direction
376               DO jj = 1, jpjm1
377                  DO ji = 1, fs_jpim1
378                     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)
379                  END DO
380               END DO
381            END DO
382            !
383            DO jl = 1, jpl              !-- first guess of tracer from u-flux
384               DO jj = 2, jpjm1
385                  DO ji = fs_2, fs_jpim1
386                     ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) )              &
387                        &   + ( pu     (ji,jj   ) - pu     (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk)
388                     !
389                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
390                  END DO
391               END DO
392            END DO
393            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
394            !
395            DO jl = 1, jpl              !-- flux in y-direction
396               DO jj = 1, jpjm1
397                  DO ji = 1, fs_jpim1
398                     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)
399                  END DO
400               END DO
401            END DO
402            !
403         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
404            !
405            DO jl = 1, jpl              !-- flux in y-direction
406               DO jj = 1, jpjm1
407                  DO ji = 1, fs_jpim1
408                     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)
409                  END DO
410               END DO
411            END DO
412            !
413            DO jl = 1, jpl              !-- first guess of tracer from v-flux
414               DO jj = 2, jpjm1
415                  DO ji = fs_2, fs_jpim1
416                     ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) )  &
417                        &   + ( pv     (ji,jj   ) - pv     (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
418                     !
419                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
420                  END DO
421               END DO
422            END DO
423            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
424            !
425            DO jl = 1, jpl              !-- flux in x-direction
426               DO jj = 1, jpjm1
427                  DO ji = 1, fs_jpim1
428                     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)
429                  END DO
430               END DO
431            END DO
432            !
433         ENDIF
434         
435      ENDIF
436      !
437      DO jl = 1, jpl                    !-- after tracer with upstream scheme
438         DO jj = 2, jpjm1
439            DO ji = fs_2, fs_jpim1
440               ztra = - (   pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj  ,jl)   &
441                  &       + pfv_ups(ji,jj,jl) - pfv_ups(ji  ,jj-1,jl) ) &
442                  &   + (   pu     (ji,jj   ) - pu     (ji-1,jj     )   &
443                  &       + pv     (ji,jj   ) - pv     (ji  ,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
444               !
445               pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
446            END DO
447         END DO
448      END DO
449      CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1. )
450
451   END SUBROUTINE upstream
452
453   
454   SUBROUTINE cen2( pamsk, jt, kt, pdt, pt, pu, pv, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
455      !!---------------------------------------------------------------------
456      !!                    ***  ROUTINE cen2  ***
457      !!     
458      !! **  Purpose :   compute the high order fluxes using a centered
459      !!                 second order scheme
460      !!----------------------------------------------------------------------
461      REAL(wp)                        , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
462      INTEGER                         , INTENT(in   ) ::   jt               ! number of sub-iteration
463      INTEGER                         , INTENT(in   ) ::   kt               ! number of iteration
464      REAL(wp)                        , INTENT(in   ) ::   pdt              ! tracer time-step
465      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt               ! tracer fields
466      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
467      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt_ups           ! upstream guess of tracer
468      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes
469      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes
470      !
471      INTEGER  ::   ji, jj, jl    ! dummy loop indices
472      REAL(wp) ::   ztra          ! local scalar
473      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zpt
474      !!----------------------------------------------------------------------
475      !
476      IF( .NOT.ll_hoxy ) THEN           !** no alternate directions **!
477         !
478         DO jl = 1, jpl
479            DO jj = 1, jpjm1
480               DO ji = 1, fs_jpim1
481                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) )
482                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) )
483               END DO
484            END DO
485         END DO
486         !
487         IF    ( kn_limiter == 1 ) THEN
488            CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
489         ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN
490            CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
491            CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
492         ENDIF
493         !
494      ELSE                              !** alternate directions **!
495         !
496         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
497            !
498            DO jl = 1, jpl              !-- flux in x-direction
499               DO jj = 1, jpjm1
500                  DO ji = 1, fs_jpim1
501                     pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
502                  END DO
503               END DO
504            END DO
505            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
506
507            DO jl = 1, jpl              !-- first guess of tracer from u-flux
508               DO jj = 2, jpjm1
509                  DO ji = fs_2, fs_jpim1
510                     ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) )              &
511                        &   + ( pu    (ji,jj   ) - pu    (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk)
512                     !
513                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
514                  END DO
515               END DO
516            END DO
517            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
518
519            DO jl = 1, jpl              !-- flux in y-direction
520               DO jj = 1, jpjm1
521                  DO ji = 1, fs_jpim1
522                     pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) )
523                  END DO
524               END DO
525            END DO
526            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
527
528         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
529            !
530            DO jl = 1, jpl              !-- flux in y-direction
531               DO jj = 1, jpjm1
532                  DO ji = 1, fs_jpim1
533                     pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
534                  END DO
535               END DO
536            END DO
537            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
538            !
539            DO jl = 1, jpl              !-- first guess of tracer from v-flux
540               DO jj = 2, jpjm1
541                  DO ji = fs_2, fs_jpim1
542                     ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) )  &
543                        &   + ( pv    (ji,jj   ) - pv    (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
544                     !
545                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
546                  END DO
547               END DO
548            END DO
549            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
550            !
551            DO jl = 1, jpl              !-- flux in x-direction
552               DO jj = 1, jpjm1
553                  DO ji = 1, fs_jpim1
554                     pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) )
555                  END DO
556               END DO
557            END DO
558            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
559
560         ENDIF
561         IF( kn_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
562         
563      ENDIF
564   
565   END SUBROUTINE cen2
566
567   
568   SUBROUTINE macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, pubox, pvbox, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
569      !!---------------------------------------------------------------------
570      !!                    ***  ROUTINE macho  ***
571      !!     
572      !! **  Purpose :   compute the high order fluxes using Ultimate-Macho scheme 
573      !!
574      !! **  Method  :   ...
575      !!
576      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
577      !!----------------------------------------------------------------------
578      REAL(wp)                        , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
579      INTEGER                         , INTENT(in   ) ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2)
580      INTEGER                         , INTENT(in   ) ::   jt               ! number of sub-iteration
581      INTEGER                         , INTENT(in   ) ::   kt               ! number of iteration
582      REAL(wp)                        , INTENT(in   ) ::   pdt              ! tracer time-step
583      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt               ! tracer fields
584      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
585      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pubox, pvbox     ! upstream velocity
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_ice( 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         !
890      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
891         DO jl = 1, jpl
892            DO jj = 1, jpjm1
893               DO ji = 1, fs_jpim1
894                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
895                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
896!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj)
897                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
898                     &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
899                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
900                     &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
901               END DO
902            END DO
903         END DO
904         !
905      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
906         DO jl = 1, jpl
907            DO jj = 1, jpjm1
908               DO ji = 1, fs_jpim1
909                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
910                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
911!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj)
912                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
913                     &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
914                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
915                     &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
916               END DO
917            END DO
918         END DO
919         !
920      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
921         DO jl = 1, jpl
922            DO jj = 1, jpjm1
923               DO ji = 1, fs_jpim1
924                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
925                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
926!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj)
927                  zdy4 = zdy2 * zdy2
928                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
929                     &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
930                     &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) * (                       ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
931                     &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) &
932                     &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl)     &
933                     &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
934               END DO
935            END DO
936         END DO
937         !
938      END SELECT
939      !
940      ! if pt at v-point is negative then use the upstream value
941      !    this should not be necessary if a proper sea-ice mask is set in Ultimate
942      !    to degrade the order of the scheme when necessary (for ex. at the ice edge)
943      IF( ll_neg ) THEN
944         DO jl = 1, jpl
945            DO jj = 1, jpjm1
946               DO ji = 1, fs_jpim1
947                  IF( pt_v(ji,jj,jl) < 0._wp ) THEN
948                     pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
949                        &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
950                  ENDIF
951               END DO
952            END DO
953         END DO
954      ENDIF
955      !                                                     !-- High order flux in j-direction  --!
956      DO jl = 1, jpl
957         DO jj = 1, jpjm1
958            DO ji = 1, fs_jpim1   ! vector opt.
959               pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
960            END DO
961         END DO
962      END DO
963      !
964   END SUBROUTINE ultimate_y
965     
966
967   SUBROUTINE nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
968      !!---------------------------------------------------------------------
969      !!                    ***  ROUTINE nonosc_ice  ***
970      !!     
971      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
972      !!       scheme and the before field by a non-oscillatory algorithm
973      !!
974      !! **  Method  :   ...
975      !!----------------------------------------------------------------------
976      REAL(wp)                   , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
977      REAL(wp)                   , INTENT(in   ) ::   pdt              ! tracer time-step
978      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pu               ! ice i-velocity => u*e2
979      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pv               ! ice j-velocity => v*e1
980      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pt, pt_ups       ! before field & upstream guess of after field
981      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pfv_ups, pfu_ups ! upstream flux
982      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux
983      !
984      INTEGER  ::   ji, jj, jl    ! dummy loop indices
985      REAL(wp) ::   zpos, zneg, zbig, zup, zdo, z1_dt              ! local scalars
986      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zcoef, zzt       !   -      -
987      REAL(wp), DIMENSION(jpi,jpj    ) :: zbup, zbdo
988      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_ups, ztj_ups
989      !!----------------------------------------------------------------------
990      zbig = 1.e+40_wp
991     
992      ! antidiffusive flux : high order minus low order
993      ! --------------------------------------------------
994      DO jl = 1, jpl
995         DO jj = 1, jpjm1
996            DO ji = 1, fs_jpim1   ! vector opt.
997               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
998               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
999            END DO
1000         END DO
1001      END DO
1002
1003      ! extreme case where pfu_ho has to be zero
1004      ! ----------------------------------------
1005      !                                    pfu_ho
1006      !                           *         --->
1007      !                        |      |  *   |        |
1008      !                        |      |      |    *   |   
1009      !                        |      |      |        |    *
1010      !            t_ups :       i-1     i       i+1       i+2   
1011      IF( ll_prelimiter_zalesak ) THEN
1012         
1013         DO jl = 1, jpl
1014            DO jj = 2, jpjm1
1015               DO ji = fs_2, fs_jpim1 
1016                  zti_ups(ji,jj,jl)= pt_ups(ji+1,jj  ,jl)
1017                  ztj_ups(ji,jj,jl)= pt_ups(ji  ,jj+1,jl)
1018               END DO
1019            END DO
1020         END DO
1021         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. )
1022
1023         DO jl = 1, jpl
1024            DO jj = 2, jpjm1
1025               DO ji = fs_2, fs_jpim1
1026                  IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj  ,jl) - pt_ups(ji,jj,jl) ) <= 0._wp .AND.  &
1027                     & pfv_ho(ji,jj,jl) * ( pt_ups(ji  ,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0._wp ) THEN
1028                     !
1029                     IF(  pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj  ,jl) - zti_ups(ji,jj,jl) ) <= 0._wp .AND.  &
1030                        & pfv_ho(ji,jj,jl) * ( ztj_ups(ji  ,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0._wp ) THEN
1031                        pfu_ho(ji,jj,jl)=0._wp
1032                        pfv_ho(ji,jj,jl)=0._wp
1033                     ENDIF
1034                     !
1035                     IF(  pfu_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji-1,jj  ,jl) ) <= 0._wp .AND.  &
1036                        & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji  ,jj-1,jl) ) <= 0._wp ) THEN
1037                        pfu_ho(ji,jj,jl)=0._wp
1038                        pfv_ho(ji,jj,jl)=0._wp
1039                     ENDIF
1040                     !
1041                  ENDIF
1042               END DO
1043            END DO
1044         END DO
1045         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
1046
1047      ENDIF
1048
1049      ! Search local extrema
1050      ! --------------------
1051      ! max/min of pt & pt_ups with large negative/positive value (-/+zbig) outside ice cover
1052      z1_dt = 1._wp / pdt
1053      DO jl = 1, jpl
1054         
1055         DO jj = 1, jpj
1056            DO ji = 1, jpi
1057               IF    ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
1058                  zbup(ji,jj) = -zbig
1059                  zbdo(ji,jj) =  zbig
1060               ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN
1061                  zbup(ji,jj) = pt_ups(ji,jj,jl)
1062                  zbdo(ji,jj) = pt_ups(ji,jj,jl)
1063               ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
1064                  zbup(ji,jj) = pt(ji,jj,jl)
1065                  zbdo(ji,jj) = pt(ji,jj,jl)
1066               ELSE
1067                  zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
1068                  zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
1069               ENDIF
1070            END DO
1071         END DO
1072
1073         DO jj = 2, jpjm1
1074            DO ji = fs_2, fs_jpim1   ! vector opt.
1075               !
1076               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
1077               zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) )
1078               !
1079               zpos = MAX( 0._wp, pfu_ho(ji-1,jj  ,jl) ) - MIN( 0._wp, pfu_ho(ji  ,jj  ,jl) ) &  ! positive/negative part of the flux
1080                  & + MAX( 0._wp, pfv_ho(ji  ,jj-1,jl) ) - MIN( 0._wp, pfv_ho(ji  ,jj  ,jl) )
1081               zneg = MAX( 0._wp, pfu_ho(ji  ,jj  ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj  ,jl) ) &
1082                  & + MAX( 0._wp, pfv_ho(ji  ,jj  ,jl) ) - MIN( 0._wp, pfv_ho(ji  ,jj-1,jl) )
1083               !
1084               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) ) &
1085                  &          ) * ( 1. - pamsk )
1086               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) ) &
1087                  &          ) * ( 1. - pamsk )
1088               !
1089               !                                  ! up & down beta terms
1090               IF( zpos > 0._wp ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt
1091               ELSE                    ; zbetup(ji,jj,jl) = 0._wp ! zbig
1092               ENDIF
1093               !
1094               IF( zneg > 0._wp ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
1095               ELSE                    ; zbetdo(ji,jj,jl) = 0._wp ! zbig
1096               ENDIF
1097               !
1098               ! if all the points are outside ice cover
1099               IF( zup == -zbig )   zbetup(ji,jj,jl) = 0._wp ! zbig
1100               IF( zdo ==  zbig )   zbetdo(ji,jj,jl) = 0._wp ! zbig           
1101               !
1102            END DO
1103         END DO
1104      END DO
1105      CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
1106
1107     
1108      ! monotonic flux in the y direction
1109      ! ---------------------------------
1110      DO jl = 1, jpl
1111         DO jj = 1, jpjm1
1112            DO ji = 1, fs_jpim1   ! vector opt.
1113               zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) )
1114               zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) )
1115               zcu = 0.5_wp + SIGN( 0.5_wp , pfu_ho(ji,jj,jl) )
1116               !
1117               zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
1118               !
1119               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl)
1120               !
1121            END DO
1122         END DO
1123
1124         DO jj = 1, jpjm1
1125            DO ji = 1, fs_jpim1   ! vector opt.
1126               zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) )
1127               zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) )
1128               zcv = 0.5_wp + SIGN( 0.5_wp , pfv_ho(ji,jj,jl) )
1129               !
1130               zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
1131               !
1132               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl)
1133               !
1134            END DO
1135         END DO
1136
1137         ! clem test
1138!!         DO jj = 2, jpjm1
1139!!            DO ji = 2, fs_jpim1   ! vector opt.
1140!!               zzt = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1141!!                  &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1142!!                  &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1143!!                  &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1144!!                  &         ) * tmask(ji,jj,1)
1145!!               IF( zzt < -epsi20 ) THEN
1146!!                  WRITE(numout,*) 'T<0 nonosc_ice',zzt
1147!!               ENDIF
1148!!            END DO
1149!!         END DO
1150
1151      END DO
1152      !
1153   END SUBROUTINE nonosc_ice
1154
1155   
1156   SUBROUTINE limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
1157      !!---------------------------------------------------------------------
1158      !!                    ***  ROUTINE limiter_x  ***
1159      !!     
1160      !! **  Purpose :   compute flux limiter
1161      !!----------------------------------------------------------------------
1162      REAL(wp)                  , INTENT(in   ) ::   pdt          ! tracer time-step
1163      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu           ! ice i-velocity => u*e2
1164      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt           ! ice tracer
1165      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pfu_ups      ! upstream flux
1166      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pfu_ho       ! high order flux
1167      !
1168      REAL(wp) ::   Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr
1169      INTEGER  ::   ji, jj, jl    ! dummy loop indices
1170      REAL(wp), DIMENSION (jpi,jpj,jpl) ::   zslpx       ! tracer slopes
1171      !!----------------------------------------------------------------------
1172      !
1173      DO jl = 1, jpl
1174         DO jj = 2, jpjm1
1175            DO ji = fs_2, fs_jpim1   ! vector opt.
1176               zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1)
1177            END DO
1178         END DO
1179      END DO
1180      CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.)   ! lateral boundary cond.
1181     
1182      DO jl = 1, jpl
1183         DO jj = 2, jpjm1
1184            DO ji = fs_2, fs_jpim1   ! vector opt.
1185               uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
1186               
1187               Rjm = zslpx(ji-1,jj,jl)
1188               Rj  = zslpx(ji  ,jj,jl)
1189               Rjp = zslpx(ji+1,jj,jl)
1190
1191               IF( kn_limiter == 3 ) THEN
1192
1193                  IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm
1194                  ELSE                        ;   Rr = Rjp
1195                  ENDIF
1196
1197                  zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)     
1198                  IF( Rj > 0. ) THEN
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                  ELSE
1202                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)),  &
1203                        &        MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
1204                  ENDIF
1205                  pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter
1206
1207               ELSEIF( kn_limiter == 2 ) THEN
1208                  IF( Rj /= 0. ) THEN
1209                     IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
1210                     ELSE                        ;   Cr = Rjp / Rj
1211                     ENDIF
1212                  ELSE
1213                     Cr = 0.
1214                  ENDIF
1215
1216                  ! -- superbee --
1217                  zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1218                  ! -- van albada 2 --
1219                  !!zpsi = 2.*Cr / (Cr*Cr+1.)
1220                  ! -- sweby (with beta=1) --
1221                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1222                  ! -- van Leer --
1223                  !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1224                  ! -- ospre --
1225                  !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
1226                  ! -- koren --
1227                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
1228                  ! -- charm --
1229                  !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
1230                  !ELSE                 ;   zpsi = 0.
1231                  !ENDIF
1232                  ! -- van albada 1 --
1233                  !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
1234                  ! -- smart --
1235                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
1236                  ! -- umist --
1237                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1238
1239                  ! high order flux corrected by the limiter
1240                  pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5
1241
1242               ENDIF
1243            END DO
1244         END DO
1245      END DO
1246      CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.)   ! lateral boundary cond.
1247      !
1248   END SUBROUTINE limiter_x
1249
1250   
1251   SUBROUTINE limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
1252      !!---------------------------------------------------------------------
1253      !!                    ***  ROUTINE limiter_y  ***
1254      !!     
1255      !! **  Purpose :   compute flux limiter
1256      !!----------------------------------------------------------------------
1257      REAL(wp)                   , INTENT(in   ) ::   pdt          ! tracer time-step
1258      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pv           ! ice i-velocity => u*e2
1259      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pt           ! ice tracer
1260      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pfv_ups      ! upstream flux
1261      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ho       ! high order flux
1262      !
1263      REAL(wp) ::   Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr
1264      INTEGER  ::   ji, jj, jl    ! dummy loop indices
1265      REAL(wp), DIMENSION (jpi,jpj,jpl) ::   zslpy       ! tracer slopes
1266      !!----------------------------------------------------------------------
1267      !
1268      DO jl = 1, jpl
1269         DO jj = 2, jpjm1
1270            DO ji = fs_2, fs_jpim1   ! vector opt.
1271               zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1)
1272            END DO
1273         END DO
1274      END DO
1275      CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.)   ! lateral boundary cond.
1276
1277      DO jl = 1, jpl
1278         DO jj = 2, jpjm1
1279            DO ji = fs_2, fs_jpim1   ! vector opt.
1280               vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
1281
1282               Rjm = zslpy(ji,jj-1,jl)
1283               Rj  = zslpy(ji,jj  ,jl)
1284               Rjp = zslpy(ji,jj+1,jl)
1285
1286               IF( kn_limiter == 3 ) THEN
1287
1288                  IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm
1289                  ELSE                        ;   Rr = Rjp
1290                  ENDIF
1291
1292                  zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)     
1293                  IF( Rj > 0. ) THEN
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                  ELSE
1297                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)),  &
1298                        &        MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
1299                  ENDIF
1300                  pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter
1301
1302               ELSEIF( kn_limiter == 2 ) THEN
1303
1304                  IF( Rj /= 0. ) THEN
1305                     IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
1306                     ELSE                        ;   Cr = Rjp / Rj
1307                     ENDIF
1308                  ELSE
1309                     Cr = 0.
1310                  ENDIF
1311
1312                  ! -- superbee --
1313                  zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1314                  ! -- van albada 2 --
1315                  !!zpsi = 2.*Cr / (Cr*Cr+1.)
1316                  ! -- sweby (with beta=1) --
1317                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1318                  ! -- van Leer --
1319                  !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1320                  ! -- ospre --
1321                  !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
1322                  ! -- koren --
1323                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
1324                  ! -- charm --
1325                  !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
1326                  !ELSE                 ;   zpsi = 0.
1327                  !ENDIF
1328                  ! -- van albada 1 --
1329                  !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
1330                  ! -- smart --
1331                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
1332                  ! -- umist --
1333                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1334
1335                  ! high order flux corrected by the limiter
1336                  pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5
1337
1338               ENDIF
1339            END DO
1340         END DO
1341      END DO
1342      CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.)   ! lateral boundary cond.
1343      !
1344   END SUBROUTINE limiter_y
1345
1346#else
1347   !!----------------------------------------------------------------------
1348   !!   Default option           Dummy module         NO SI3 sea-ice model
1349   !!----------------------------------------------------------------------
1350#endif
1351
1352   !!======================================================================
1353END MODULE icedyn_adv_umx
Note: See TracBrowser for help on using the repository browser.