source: NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icedyn_adv_umx.F90 @ 10888

Last change on this file since 10888 was 10888, checked in by davestorkey, 20 months ago

branches/UKMO/NEMO_4.0_mirror : clear SVN keywords

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