source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn_adv_umx.F90 @ 12475

Last change on this file since 12475 was 12475, checked in by dancopsey, 14 months ago
  • Add more print statements.
  • Move away from using snow to ice diagnostics and use a new snow to pond one instead.
File size: 86.1 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 fields
14   !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders
15   !!   macho             : compute the fluxes
16   !!   nonosc_ice        : limit the fluxes using 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   INTEGER, PARAMETER ::   np_advS = 1         ! advection for S and T:    dVS/dt = -div(      uVS     ) => np_advS = 1
36   !                                                                    or dVS/dt = -div( uA * uHS / u ) => np_advS = 2
37   !                                                                    or dVS/dt = -div( uV * uS  / u ) => np_advS = 3
38   INTEGER, PARAMETER ::   np_limiter = 1      ! limiter: 1 = nonosc
39   !                                                      2 = superbee
40   !                                                      3 = h3
41   LOGICAL            ::   ll_upsxy  = .TRUE.   ! alternate directions for upstream
42   LOGICAL            ::   ll_hoxy   = .TRUE.   ! alternate directions for high order
43   LOGICAL            ::   ll_neg    = .TRUE.   ! if T interpolated at u/v points is negative or v_i < 1.e-6
44   !                                                 then interpolate T at u/v points using the upstream scheme
45   LOGICAL            ::   ll_prelim = .FALSE.  ! prelimiter from: Zalesak(1979) eq. 14 => not well defined in 2D
46   !
47   REAL(wp)           ::   z1_6   = 1._wp /   6._wp   ! =1/6
48   REAL(wp)           ::   z1_120 = 1._wp / 120._wp   ! =1/120
49   !
50   INTEGER, ALLOCATABLE, DIMENSION(:,:,:) ::   imsk_small, jmsk_small
51   !
52   !! * Substitutions
53#  include "vectopt_loop_substitute.h90"
54   !!----------------------------------------------------------------------
55   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
56   !! $Id$
57   !! Software governed by the CeCILL licence     (./LICENSE)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  &
62      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
63      !!----------------------------------------------------------------------
64      !!                  ***  ROUTINE ice_dyn_adv_umx  ***
65      !!
66      !! **  Purpose :   Compute the now trend due to total advection of
67      !!                 tracers and add it to the general trend of tracer equations
68      !!                 using an "Ultimate-Macho" scheme
69      !!
70      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
71      !!----------------------------------------------------------------------
72      INTEGER                     , INTENT(in   ) ::   kn_umx     ! order of the scheme (1-5=UM or 20=CEN2)
73      INTEGER                     , INTENT(in   ) ::   kt         ! time step
74      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
75      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
76      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_i       ! ice thickness
77      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_s       ! snw thickness
78      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_ip      ! ice pond thickness
79      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
80      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
81      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
82      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
83      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
84      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
85      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
86      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
87      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness
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, zvi_cen
95      REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! for global communication
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) ::   zu_cat, zv_cat
99      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho, zua_ups, zva_ups
100      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai , z1_aip, zhvar
101      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max
102      !
103      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs 
104      !!----------------------------------------------------------------------
105      !
106      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
107
108      write(numout,*)'ice_dyn_adv_umx 1:psv_i, pe_s = ',psv_i(3,4,1), ' ', pe_i (3,4,1,1)
109      !
110      ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- !
111      DO jl = 1, jpl
112         DO jj = 2, jpjm1
113            DO ji = fs_2, fs_jpim1
114               zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), &
115                  &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), &
116                  &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), &
117                  &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) )
118               zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), &
119                  &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), &
120                  &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), &
121                  &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) )
122               zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), &
123                  &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), &
124                  &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), &
125                  &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) )
126            END DO
127         END DO
128      END DO
129      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. )
130      !
131      !
132      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- !
133      !        Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
134      !              this should not affect too much the stability
135      zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
136      zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
137     
138      ! non-blocking global communication send zcflnow and receive zcflprv
139      CALL mpp_delay_max( 'icedyn_adv_umx', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
140
141      IF( zcflprv(1) > .5 ) THEN   ;   icycle = 2
142      ELSE                         ;   icycle = 1
143      ENDIF
144      zdt = rdt_ice / REAL(icycle)
145
146      write(numout,*)'ice_dyn_adv_umx 2: pu_ice, e2u = ',pu_ice(3,4), e2u(3,4)
147
148      ! --- transport --- !
149      zudy(:,:) = pu_ice(:,:) * e2u(:,:)
150      zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
151
152      write(numout,*)'ice_dyn_adv_umx 2: zudy, zvdx, zu_cat, zv_cat, zcu_box, zcv_box = ',zudy(3,4), zvdx(3,4), zu_cat(3,4,1), zv_cat(3,4,1), zcu_box(3,4), zcv_box(3,4)
153
154      !
155      ! setup transport for each ice cat
156      DO jl = 1, jpl
157         zu_cat(:,:,jl) = zudy(:,:)
158         zv_cat(:,:,jl) = zvdx(:,:)
159      END DO
160
161      write(numout,*)'ice_dyn_adv_umx 3: zudy, zvdx, zu_cat, zv_cat, zcu_box, zcv_box = ',zudy(3,4), zvdx(3,4), zu_cat(3,4,1), zv_cat(3,4,1), zcu_box(3,4), zcv_box(3,4)
162
163      !
164      ! --- define velocity for advection: u*grad(H) --- !
165      DO jj = 2, jpjm1
166         DO ji = fs_2, fs_jpim1
167            IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
168            ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj)
169            ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj)
170            ENDIF
171
172            IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
173            ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1)
174            ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  )
175            ENDIF
176         END DO
177      END DO
178
179      write(numout,*)'ice_dyn_adv_umx 4: zudy, zvdx, zu_cat, zv_cat, zcu_box, zcv_box = ',zudy(3,4), zvdx(3,4), zu_cat(3,4,1), zv_cat(3,4,1), zcu_box(3,4), zcv_box(3,4)
180
181      !---------------!
182      !== advection ==!
183      !---------------!
184      DO jt = 1, icycle
185
186         ! record at_i before advection (for open water)
187         zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
188         
189         ! inverse of A and Ap
190         WHERE( pa_i(:,:,:) >= epsi20 )   ;   z1_ai(:,:,:) = 1._wp / pa_i(:,:,:)
191         ELSEWHERE                        ;   z1_ai(:,:,:) = 0.
192         END WHERE
193         WHERE( pa_ip(:,:,:) >= epsi20 )  ;   z1_aip(:,:,:) = 1._wp / pa_ip(:,:,:)
194         ELSEWHERE                        ;   z1_aip(:,:,:) = 0.
195         END WHERE
196         !
197         ! setup a mask where advection will be upstream
198         IF( ll_neg ) THEN
199            IF( .NOT. ALLOCATED(imsk_small) )   ALLOCATE( imsk_small(jpi,jpj,jpl) ) 
200            IF( .NOT. ALLOCATED(jmsk_small) )   ALLOCATE( jmsk_small(jpi,jpj,jpl) ) 
201            DO jl = 1, jpl
202               DO jj = 1, jpjm1
203                  DO ji = 1, jpim1
204                     zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) )
205                     IF( zvi_cen < epsi06) THEN   ;   imsk_small(ji,jj,jl) = 0
206                     ELSE                         ;   imsk_small(ji,jj,jl) = 1   ;   ENDIF
207                     zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) )
208                     IF( zvi_cen < epsi06) THEN   ;   jmsk_small(ji,jj,jl) = 0
209                     ELSE                         ;   jmsk_small(ji,jj,jl) = 1   ;   ENDIF
210                  END DO
211               END DO
212            END DO
213         ENDIF
214         !
215         ! ----------------------- !
216         ! ==> start advection <== !
217         ! ----------------------- !
218         !
219         !== Ice area ==!
220         zamsk = 1._wp
221         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zu_cat , zv_cat , zcu_box, zcv_box, &
222            &                                      pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho )
223
224         write(numout,*)'ice_dyn_adv_umx 5: zudy, zvdx, zu_cat, zv_cat, zcu_box, zcv_box = ',zudy(3,4), zvdx(3,4), zu_cat(3,4,1), zv_cat(3,4,1), zcu_box(3,4), zcv_box(3,4)
225
226         !
227         !                             ! --------------------------------- !
228         IF( np_advS == 1 ) THEN       ! -- advection form: -div( uVS ) -- !
229            !                          ! --------------------------------- !
230            zamsk = 0._wp
231            !== Ice volume ==!
232            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
233            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
234               &                                      zhvar, pv_i, zua_ups, zva_ups )
235            !== Snw volume ==!         
236            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
237            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
238               &                                      zhvar, pv_s, zua_ups, zva_ups )
239            !
240            zamsk = 1._wp
241            !== Salt content ==!
242            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
243               &                                      psv_i, psv_i )
244            write(numout,*)'ice_dyn_adv_umx 6:psv_i, pe_s, jt = ',psv_i(3,4,1), ' ', pe_i (3,4,1,1), ' ', jt
245            !== Ice heat content ==!
246            DO jk = 1, nlay_i
247               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
248                  &                                      pe_i(:,:,jk,:), pe_i(:,:,jk,:) )
249            END DO
250            !== Snw heat content ==!
251            DO jk = 1, nlay_s
252               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
253                  &                                      pe_s(:,:,jk,:), pe_s(:,:,jk,:) )
254            END DO
255            !
256            !                          ! ------------------------------------------ !
257         ELSEIF( np_advS == 2 ) THEN   ! -- advection form: -div( uA * uHS / u ) -- !
258            !                          ! ------------------------------------------ !
259            zamsk = 0._wp
260            !== Ice volume ==!
261            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
262            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
263               &                                      zhvar, pv_i, zua_ups, zva_ups )
264            !== Snw volume ==!         
265            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
266            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
267               &                                      zhvar, pv_s, zua_ups, zva_ups )
268            !== Salt content ==!
269            zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
270            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
271               &                                      zhvar, psv_i, zua_ups, zva_ups )
272            write(numout,*)'ice_dyn_adv_umx 7:psv_i, pe_s, jt = ',psv_i(3,4,1), ' ', pe_i (3,4,1,1), ' ', jt
273            !== Ice heat content ==!
274            DO jk = 1, nlay_i
275               zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
276               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, &
277                  &                                      zhvar, pe_i(:,:,jk,:), zua_ups, zva_ups )
278            END DO
279            !== Snw heat content ==!
280            DO jk = 1, nlay_s
281               zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
282               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, &
283                  &                                      zhvar, pe_s(:,:,jk,:), zua_ups, zva_ups )
284            END DO
285            !
286            !                          ! ----------------------------------------- !
287         ELSEIF( np_advS == 3 ) THEN   ! -- advection form: -div( uV * uS / u ) -- !
288            !                          ! ----------------------------------------- !
289            zamsk = 0._wp
290            !
291            ALLOCATE( zuv_ho (jpi,jpj,jpl), zvv_ho (jpi,jpj,jpl),  &
292               &      zuv_ups(jpi,jpj,jpl), zvv_ups(jpi,jpj,jpl), z1_vi(jpi,jpj,jpl), z1_vs(jpi,jpj,jpl) )
293            !
294            ! inverse of Vi
295            WHERE( pv_i(:,:,:) >= epsi20 )   ;   z1_vi(:,:,:) = 1._wp / pv_i(:,:,:)
296            ELSEWHERE                        ;   z1_vi(:,:,:) = 0.
297            END WHERE
298            ! inverse of Vs
299            WHERE( pv_s(:,:,:) >= epsi20 )   ;   z1_vs(:,:,:) = 1._wp / pv_s(:,:,:)
300            ELSEWHERE                        ;   z1_vs(:,:,:) = 0.
301            END WHERE
302            !
303            ! It is important to first calculate the ice fields and then the snow fields (because we use the same arrays)
304            !
305            !== Ice volume ==!
306            zuv_ups = zua_ups
307            zvv_ups = zva_ups
308            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
309            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
310               &                                      zhvar, pv_i, zuv_ups, zvv_ups, zuv_ho , zvv_ho )
311            !== Salt content ==!
312            zhvar(:,:,:) = psv_i(:,:,:) * z1_vi(:,:,:)
313            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, &
314               &                                      zhvar, psv_i, zuv_ups, zvv_ups )
315
316            write(numout,*)'ice_dyn_adv_umx 8:psv_i, pe_s, jt = ',psv_i(3,4,1), ' ', pe_i (3,4,1,1), ' ', jt
317            !== Ice heat content ==!
318            DO jk = 1, nlay_i
319               zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_vi(:,:,:)
320               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, &
321                  &                                      zhvar, pe_i(:,:,jk,:), zuv_ups, zvv_ups )
322            END DO
323            !== Snow volume ==!         
324            zuv_ups = zua_ups
325            zvv_ups = zva_ups
326            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
327            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
328               &                                      zhvar, pv_s, zuv_ups, zvv_ups, zuv_ho , zvv_ho )
329            !== Snw heat content ==!
330            DO jk = 1, nlay_s
331               zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_vs(:,:,:)
332               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, &
333                  &                                      zhvar, pe_s(:,:,jk,:), zuv_ups, zvv_ups )
334            END DO
335            !
336            DEALLOCATE( zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs )
337            !
338         ENDIF
339         !
340         !== Ice age ==!
341         IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN
342            zamsk = 1._wp
343            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
344               &                                      poa_i, poa_i )
345         ENDIF
346         !
347         !== melt ponds ==!
348         IF ( ln_pnd_H12 ) THEN
349            ! fraction
350            zamsk = 1._wp
351            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat , zv_cat , zcu_box, zcv_box, &
352               &                                      pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho )
353            ! volume
354            zamsk = 0._wp
355            zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:)
356            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
357               &                                      zhvar, pv_ip, zua_ups, zva_ups )
358            ! lid thickness
359            zamsk = 0._wp
360            zhvar(:,:,:) = plh_ip(:,:,:) * z1_aip(:,:,:)
361            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
362               &                                      zhvar, plh_ip, zua_ups, zva_ups )
363           
364         ENDIF
365         !
366         !== Open water area ==!
367         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
368         DO jj = 2, jpjm1
369            DO ji = fs_2, fs_jpim1
370               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & 
371                  &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
372            END DO
373         END DO
374         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1. )
375         !
376         !
377         ! --- Ensure non-negative fields and in-bound thicknesses --- !
378         ! Remove negative values (conservation is ensured)
379         !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
380         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
381
382         write(numout,*)'ice_dyn_adv_umx 9:psv_i, pe_s, jt = ',psv_i(3,4,1), ' ', pe_i (3,4,1,1), ' ', jt
383         !
384         ! Make sure ice thickness is not too big
385         !    (because ice thickness can be too large where ice concentration is very small)
386         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
387
388         write(numout,*)'ice_dyn_adv_umx 10:psv_i, pe_s, jt = ',psv_i(3,4,1), ' ', pe_i (3,4,1,1), ' ', jt
389
390      END DO
391
392      write(numout,*)'ice_dyn_adv_umx 11:psv_i, pe_s = ',psv_i(3,4,1), ' ', pe_i (3,4,1,1)
393      !
394   END SUBROUTINE ice_dyn_adv_umx
395
396   
397   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox,  &
398      &                                            pt, ptc, pua_ups, pva_ups, pua_ho, pva_ho )
399      !!----------------------------------------------------------------------
400      !!                  ***  ROUTINE adv_umx  ***
401      !!
402      !! **  Purpose :   Compute the now trend due to total advection of
403      !!                 tracers and add it to the general trend of tracer equations
404      !!
405      !! **  Method  :   - calculate upstream fluxes and upstream solution for tracers V/A(=H) etc
406      !!                 - calculate tracer H at u and v points (Ultimate)
407      !!                 - calculate the high order fluxes using alterning directions (Macho)
408      !!                 - apply a limiter on the fluxes (nonosc_ice)
409      !!                 - convert this tracer flux to a "volume" flux (uH -> uV)
410      !!                 - apply a limiter a second time on the volumes fluxes (nonosc_ice)
411      !!                 - calculate the high order solution for V
412      !!
413      !! ** Action : solve 3 equations => a) dA/dt  = -div(uA)
414      !!                                  b) dV/dt  = -div(uV)  using dH/dt = -u.grad(H)
415      !!                                  c) dVS/dt = -div(uVS) using either dHS/dt = -u.grad(HS) or dS/dt = -u.grad(S)
416      !!
417      !!             in eq. b), - fluxes uH are evaluated (with UMx) and limited with nonosc_ice. This step is necessary to get a good H.
418      !!                        - then we convert this flux to a "volume" flux this way => uH * uA / u
419      !!                             where uA is the flux from eq. a)
420      !!                             this "volume" flux is also limited with nonosc_ice (otherwise overshoots can occur)
421      !!                        - at last we estimate dV/dt = -div(uH * uA / u)
422      !!
423      !!             in eq. c), one can solve the equation for  S (ln_advS=T), then dVS/dt = -div(uV * uS  / u)
424      !!                                                or for HS (ln_advS=F), then dVS/dt = -div(uA * uHS / u)
425      !!
426      !! ** Note : - this method can lead to tiny negative V (-1.e-20) => set it to 0 while conserving mass etc.
427      !!           - At the ice edge, Ultimate scheme can lead to:
428      !!                              1) negative interpolated tracers at u-v points
429      !!                              2) non-zero interpolated tracers at u-v points eventhough there is no ice and velocity is outward
430      !!                              Solution for 1): apply an upstream scheme when it occurs. A better solution would be to degrade the order of
431      !!                                               the scheme automatically by applying a mask of the ice cover inside Ultimate (not done).
432      !!                              Solution for 2): we set it to 0 in this case
433      !!           - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good.
434      !!             Large values of H can appear for very small ice concentration, and when it does it messes the things up since we
435      !!             work on H (and not V). It is partly related to the multi-category approach
436      !!             Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice
437      !!             concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter
438      !!             since sv_i and e_i are still good.
439      !!----------------------------------------------------------------------
440      REAL(wp)                        , INTENT(in   )           ::   pamsk            ! advection of concentration (1) or other tracers (0)
441      INTEGER                         , INTENT(in   )           ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2)
442      INTEGER                         , INTENT(in   )           ::   jt               ! number of sub-iteration
443      INTEGER                         , INTENT(in   )           ::   kt               ! number of iteration
444      REAL(wp)                        , INTENT(in   )           ::   pdt              ! tracer time-step
445      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pu   , pv        ! 2 ice velocity components => u*e2
446      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   )           ::   puc  , pvc       ! 2 ice velocity components => u*e2 or u*a*e2u
447      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pubox, pvbox     ! upstream velocity
448      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   pt               ! tracer field
449      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   ptc              ! tracer content field
450      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout), OPTIONAL ::   pua_ups, pva_ups ! upstream u*a fluxes
451      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho   ! high order u*a fluxes
452      !
453      INTEGER  ::   ji, jj, jl       ! dummy loop indices 
454      REAL(wp) ::   ztra             ! local scalar
455      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ho , zfv_ho , zpt
456      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ups, zfv_ups, zt_ups
457      !!----------------------------------------------------------------------
458      !
459      ! Upstream (_ups) fluxes
460      ! -----------------------
461      CALL upstream( pamsk, jt, kt, pdt, pt, pu, pv, zt_ups, zfu_ups, zfv_ups )
462     
463      ! High order (_ho) fluxes
464      ! -----------------------
465      SELECT CASE( kn_umx )
466         !
467      CASE ( 20 )                          !== centered second order ==!
468         !
469         CALL cen2( pamsk, jt, kt, pdt, pt, pu, pv, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
470         !
471      CASE ( 1:5 )                         !== 1st to 5th order ULTIMATE-MACHO scheme ==!
472         !
473         CALL macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, pubox, pvbox, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
474         !
475      END SELECT
476      !
477      !              --ho    --ho
478      ! new fluxes = u*H  *  u*a / u
479      ! ----------------------------
480      IF( pamsk == 0._wp ) THEN
481         DO jl = 1, jpl
482            DO jj = 1, jpjm1
483               DO ji = 1, fs_jpim1
484                  IF( ABS( pu(ji,jj) ) > epsi10 ) THEN
485                     zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj)
486                     zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj)
487                  ELSE
488                     zfu_ho (ji,jj,jl) = 0._wp
489                     zfu_ups(ji,jj,jl) = 0._wp
490                  ENDIF
491                  !
492                  IF( ABS( pv(ji,jj) ) > epsi10 ) THEN
493                     zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj)
494                     zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj)
495                  ELSE
496                     zfv_ho (ji,jj,jl) = 0._wp 
497                     zfv_ups(ji,jj,jl) = 0._wp 
498                  ENDIF
499               END DO
500            END DO
501         END DO
502
503         ! the new "volume" fluxes must also be "flux corrected"
504         ! thus we calculate the upstream solution and apply a limiter again
505         DO jl = 1, jpl
506            DO jj = 2, jpjm1
507               DO ji = fs_2, fs_jpim1
508                  ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) )
509                  !
510                  zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
511               END DO
512            END DO
513         END DO
514         CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T',  1. )
515         !
516         IF    ( np_limiter == 1 ) THEN
517            CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
518         ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN
519            CALL limiter_x( pdt, pu, ptc, zfu_ups, zfu_ho )
520            CALL limiter_y( pdt, pv, ptc, zfv_ups, zfv_ho )
521         ENDIF
522         !
523      ENDIF
524      !                                   --ho    --ups
525      ! in case of advection of A: output u*a and u*a
526      ! -----------------------------------------------
527      IF( PRESENT( pua_ho ) ) THEN
528         DO jl = 1, jpl
529            DO jj = 1, jpjm1
530               DO ji = 1, fs_jpim1
531                  pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl)
532                  pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl)
533              END DO
534            END DO
535         END DO
536      ENDIF
537      !
538      ! final trend with corrected fluxes
539      ! ---------------------------------
540      DO jl = 1, jpl
541         DO jj = 2, jpjm1
542            DO ji = fs_2, fs_jpim1 
543               ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) 
544               !
545               ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)               
546            END DO
547         END DO
548      END DO
549      CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1. )
550      !
551   END SUBROUTINE adv_umx
552
553
554   SUBROUTINE upstream( pamsk, jt, kt, pdt, pt, pu, pv, pt_ups, pfu_ups, pfv_ups )
555      !!---------------------------------------------------------------------
556      !!                    ***  ROUTINE upstream  ***
557      !!     
558      !! **  Purpose :   compute the upstream fluxes and upstream guess of tracer
559      !!----------------------------------------------------------------------
560      REAL(wp)                        , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
561      INTEGER                         , INTENT(in   ) ::   jt               ! number of sub-iteration
562      INTEGER                         , INTENT(in   ) ::   kt               ! number of iteration
563      REAL(wp)                        , INTENT(in   ) ::   pdt              ! tracer time-step
564      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt               ! tracer fields
565      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
566      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_ups           ! upstream guess of tracer
567      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ups, pfv_ups ! upstream fluxes
568      !
569      INTEGER  ::   ji, jj, jl    ! dummy loop indices
570      REAL(wp) ::   ztra          ! local scalar
571      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zpt
572      !!----------------------------------------------------------------------
573
574      IF( .NOT. ll_upsxy ) THEN         !** no alternate directions **!
575         !
576         DO jl = 1, jpl
577            DO jj = 1, jpjm1
578               DO ji = 1, fs_jpim1
579                  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)
580                  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)
581               END DO
582            END DO
583         END DO
584         !
585      ELSE                              !** alternate directions **!
586         !
587         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
588            !
589            DO jl = 1, jpl              !-- flux in x-direction
590               DO jj = 1, jpjm1
591                  DO ji = 1, fs_jpim1
592                     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)
593                  END DO
594               END DO
595            END DO
596            !
597            DO jl = 1, jpl              !-- first guess of tracer from u-flux
598               DO jj = 2, jpjm1
599                  DO ji = fs_2, fs_jpim1
600                     ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) )              &
601                        &   + ( pu     (ji,jj   ) - pu     (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk)
602                     !
603                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
604                  END DO
605               END DO
606            END DO
607            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
608            !
609            DO jl = 1, jpl              !-- flux in y-direction
610               DO jj = 1, jpjm1
611                  DO ji = 1, fs_jpim1
612                     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)
613                  END DO
614               END DO
615            END DO
616            !
617         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
618            !
619            DO jl = 1, jpl              !-- flux in y-direction
620               DO jj = 1, jpjm1
621                  DO ji = 1, fs_jpim1
622                     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)
623                  END DO
624               END DO
625            END DO
626            !
627            DO jl = 1, jpl              !-- first guess of tracer from v-flux
628               DO jj = 2, jpjm1
629                  DO ji = fs_2, fs_jpim1
630                     ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) )  &
631                        &   + ( pv     (ji,jj   ) - pv     (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
632                     !
633                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
634                  END DO
635               END DO
636            END DO
637            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
638            !
639            DO jl = 1, jpl              !-- flux in x-direction
640               DO jj = 1, jpjm1
641                  DO ji = 1, fs_jpim1
642                     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)
643                  END DO
644               END DO
645            END DO
646            !
647         ENDIF
648         
649      ENDIF
650      !
651      DO jl = 1, jpl                    !-- after tracer with upstream scheme
652         DO jj = 2, jpjm1
653            DO ji = fs_2, fs_jpim1
654               ztra = - (   pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj  ,jl)   &
655                  &       + pfv_ups(ji,jj,jl) - pfv_ups(ji  ,jj-1,jl) ) &
656                  &   + (   pu     (ji,jj   ) - pu     (ji-1,jj     )   &
657                  &       + pv     (ji,jj   ) - pv     (ji  ,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
658               !
659               pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
660            END DO
661         END DO
662      END DO
663      CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1. )
664
665   END SUBROUTINE upstream
666
667   
668   SUBROUTINE cen2( pamsk, jt, kt, pdt, pt, pu, pv, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
669      !!---------------------------------------------------------------------
670      !!                    ***  ROUTINE cen2  ***
671      !!     
672      !! **  Purpose :   compute the high order fluxes using a centered
673      !!                 second order scheme
674      !!----------------------------------------------------------------------
675      REAL(wp)                        , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
676      INTEGER                         , INTENT(in   ) ::   jt               ! number of sub-iteration
677      INTEGER                         , INTENT(in   ) ::   kt               ! number of iteration
678      REAL(wp)                        , INTENT(in   ) ::   pdt              ! tracer time-step
679      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt               ! tracer fields
680      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
681      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt_ups           ! upstream guess of tracer
682      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes
683      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes
684      !
685      INTEGER  ::   ji, jj, jl    ! dummy loop indices
686      REAL(wp) ::   ztra          ! local scalar
687      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zpt
688      !!----------------------------------------------------------------------
689      !
690      IF( .NOT.ll_hoxy ) THEN           !** no alternate directions **!
691         !
692         DO jl = 1, jpl
693            DO jj = 1, jpjm1
694               DO ji = 1, fs_jpim1
695                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) )
696                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) )
697               END DO
698            END DO
699         END DO
700         !
701         IF    ( np_limiter == 1 ) THEN
702            CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
703         ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN
704            CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
705            CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
706         ENDIF
707         !
708      ELSE                              !** alternate directions **!
709         !
710         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
711            !
712            DO jl = 1, jpl              !-- flux in x-direction
713               DO jj = 1, jpjm1
714                  DO ji = 1, fs_jpim1
715                     pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
716                  END DO
717               END DO
718            END DO
719            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
720
721            DO jl = 1, jpl              !-- first guess of tracer from u-flux
722               DO jj = 2, jpjm1
723                  DO ji = fs_2, fs_jpim1
724                     ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) )              &
725                        &   + ( pu    (ji,jj   ) - pu    (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk)
726                     !
727                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
728                  END DO
729               END DO
730            END DO
731            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
732
733            DO jl = 1, jpl              !-- flux in y-direction
734               DO jj = 1, jpjm1
735                  DO ji = 1, fs_jpim1
736                     pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) )
737                  END DO
738               END DO
739            END DO
740            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
741
742         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
743            !
744            DO jl = 1, jpl              !-- flux in y-direction
745               DO jj = 1, jpjm1
746                  DO ji = 1, fs_jpim1
747                     pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
748                  END DO
749               END DO
750            END DO
751            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
752            !
753            DO jl = 1, jpl              !-- first guess of tracer from v-flux
754               DO jj = 2, jpjm1
755                  DO ji = fs_2, fs_jpim1
756                     ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) )  &
757                        &   + ( pv    (ji,jj   ) - pv    (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
758                     !
759                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
760                  END DO
761               END DO
762            END DO
763            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
764            !
765            DO jl = 1, jpl              !-- flux in x-direction
766               DO jj = 1, jpjm1
767                  DO ji = 1, fs_jpim1
768                     pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) )
769                  END DO
770               END DO
771            END DO
772            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
773
774         ENDIF
775         IF( np_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
776         
777      ENDIF
778   
779   END SUBROUTINE cen2
780
781   
782   SUBROUTINE macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, pubox, pvbox, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
783      !!---------------------------------------------------------------------
784      !!                    ***  ROUTINE macho  ***
785      !!     
786      !! **  Purpose :   compute the high order fluxes using Ultimate-Macho scheme 
787      !!
788      !! **  Method  :   ...
789      !!
790      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
791      !!----------------------------------------------------------------------
792      REAL(wp)                        , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
793      INTEGER                         , INTENT(in   ) ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2)
794      INTEGER                         , INTENT(in   ) ::   jt               ! number of sub-iteration
795      INTEGER                         , INTENT(in   ) ::   kt               ! number of iteration
796      REAL(wp)                        , INTENT(in   ) ::   pdt              ! tracer time-step
797      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt               ! tracer fields
798      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
799      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pubox, pvbox     ! upstream velocity
800      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt_ups           ! upstream guess of tracer
801      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes
802      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes
803      !
804      INTEGER  ::   ji, jj, jl    ! dummy loop indices
805      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zt_u, zt_v, zpt
806      !!----------------------------------------------------------------------
807      !
808      IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
809         !
810         !                                                        !--  ultimate interpolation of pt at u-point  --!
811         CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho )
812         !                                                        !--  limiter in x --!
813         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
814         !                                                        !--  advective form update in zpt  --!
815         DO jl = 1, jpl
816            DO jj = 2, jpjm1
817               DO ji = fs_2, fs_jpim1
818                  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) &
819                     &                              + pt   (ji,jj,jl) * ( pu  (ji,jj   ) - pu  (ji-1,jj   ) ) * r1_e1e2t(ji,jj) &
820                     &                                                                                        * pamsk           &
821                     &                             ) * pdt ) * tmask(ji,jj,1)
822               END DO
823            END DO
824         END DO
825         CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
826         !
827         !                                                        !--  ultimate interpolation of pt at v-point  --!
828         IF( ll_hoxy ) THEN
829            CALL ultimate_y( pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho )
830         ELSE
831            CALL ultimate_y( pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho )
832         ENDIF
833         !                                                        !--  limiter in y --!
834         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
835         !         
836         !
837      ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
838         !
839         !                                                        !--  ultimate interpolation of pt at v-point  --!
840         CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho )
841         !                                                        !--  limiter in y --!
842         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
843         !                                                        !--  advective form update in zpt  --!
844         DO jl = 1, jpl
845            DO jj = 2, jpjm1
846               DO ji = fs_2, fs_jpim1
847                  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) &
848                     &                              + pt   (ji,jj,jl) * ( pv  (ji,jj   ) - pv  (ji,jj-1   ) ) * r1_e1e2t(ji,jj) &
849                     &                                                                                        * pamsk           &
850                     &                             ) * pdt ) * tmask(ji,jj,1) 
851               END DO
852            END DO
853         END DO
854         CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
855         !
856         !                                                        !--  ultimate interpolation of pt at u-point  --!
857         IF( ll_hoxy ) THEN
858            CALL ultimate_x( pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho )
859         ELSE
860            CALL ultimate_x( pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho )
861         ENDIF
862         !                                                        !--  limiter in x --!
863         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
864         !
865      ENDIF
866
867      IF( np_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
868      !
869   END SUBROUTINE macho
870
871
872   SUBROUTINE ultimate_x( pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho )
873      !!---------------------------------------------------------------------
874      !!                    ***  ROUTINE ultimate_x  ***
875      !!     
876      !! **  Purpose :   compute tracer at u-points
877      !!
878      !! **  Method  :   ...
879      !!
880      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
881      !!----------------------------------------------------------------------
882      REAL(wp)                        , INTENT(in   ) ::   pamsk     ! advection of concentration (1) or other tracers (0)
883      INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
884      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step
885      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pu        ! ice i-velocity component
886      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt        ! tracer fields
887      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_u      ! tracer at u-point
888      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho    ! high order flux
889      !
890      INTEGER  ::   ji, jj, jl             ! dummy loop indices
891      REAL(wp) ::   zcu, zdx2, zdx4        !   -      -
892      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztu1, ztu2, ztu3, ztu4
893      !!----------------------------------------------------------------------
894      !
895      !                                                     !--  Laplacian in i-direction  --!
896      DO jl = 1, jpl
897         DO jj = 2, jpjm1         ! First derivative (gradient)
898            DO ji = 1, fs_jpim1
899               ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
900            END DO
901            !                     ! Second derivative (Laplacian)
902            DO ji = fs_2, fs_jpim1
903               ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj)
904            END DO
905         END DO
906      END DO
907      CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1. )
908      !
909      !                                                     !--  BiLaplacian in i-direction  --!
910      DO jl = 1, jpl
911         DO jj = 2, jpjm1         ! Third derivative
912            DO ji = 1, fs_jpim1
913               ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
914            END DO
915            !                     ! Fourth derivative
916            DO ji = fs_2, fs_jpim1
917               ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj)
918            END DO
919         END DO
920      END DO
921      CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1. )
922      !
923      !
924      SELECT CASE (kn_umx )
925      !
926      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
927         !       
928         DO jl = 1, jpl
929            DO jj = 1, jpjm1
930               DO ji = 1, fs_jpim1   ! vector opt.
931                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
932                     &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
933               END DO
934            END DO
935         END DO
936         !
937      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
938         !
939         DO jl = 1, jpl
940            DO jj = 1, jpjm1
941               DO ji = 1, fs_jpim1   ! vector opt.
942                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
943                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
944                     &                                                            - zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
945               END DO
946            END DO
947         END DO
948         
949      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
950         !
951         DO jl = 1, jpl
952            DO jj = 1, jpjm1
953               DO ji = 1, fs_jpim1   ! vector opt.
954                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
955                  zdx2 = e1u(ji,jj) * e1u(ji,jj)
956!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj)
957                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
958                     &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
959                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
960                     &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
961               END DO
962            END DO
963         END DO
964         !
965      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
966         !
967         DO jl = 1, jpl
968            DO jj = 1, jpjm1
969               DO ji = 1, fs_jpim1   ! vector opt.
970                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
971                  zdx2 = e1u(ji,jj) * e1u(ji,jj)
972!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj)
973                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
974                     &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
975                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
976                     &                                                   - 0.5_wp * zcu   * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
977               END DO
978            END DO
979         END DO
980         !
981      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
982         !
983         DO jl = 1, jpl
984            DO jj = 1, jpjm1
985               DO ji = 1, fs_jpim1   ! vector opt.
986                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
987                  zdx2 = e1u(ji,jj) * e1u(ji,jj)
988!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj)
989                  zdx4 = zdx2 * zdx2
990                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (        (                       pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
991                     &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
992                     &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) * (                       ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
993                     &                                                   - 0.5_wp * zcu   * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) &
994                     &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)     &
995                     &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
996               END DO
997            END DO
998         END DO
999         !
1000      END SELECT
1001      !
1002      ! if pt at u-point is negative then use the upstream value
1003      !    this should not be necessary if a proper sea-ice mask is set in Ultimate
1004      !    to degrade the order of the scheme when necessary (for ex. at the ice edge)
1005      IF( ll_neg ) THEN
1006         DO jl = 1, jpl
1007            DO jj = 1, jpjm1
1008               DO ji = 1, fs_jpim1
1009                  IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
1010                     pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
1011                        &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
1012                  ENDIF
1013               END DO
1014            END DO
1015         END DO
1016      ENDIF
1017      !                                                     !-- High order flux in i-direction  --!
1018      DO jl = 1, jpl
1019         DO jj = 1, jpjm1
1020            DO ji = 1, fs_jpim1   ! vector opt.
1021               pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl)
1022            END DO
1023         END DO
1024      END DO
1025      !
1026   END SUBROUTINE ultimate_x
1027   
1028 
1029   SUBROUTINE ultimate_y( pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho )
1030      !!---------------------------------------------------------------------
1031      !!                    ***  ROUTINE ultimate_y  ***
1032      !!     
1033      !! **  Purpose :   compute tracer at v-points
1034      !!
1035      !! **  Method  :   ...
1036      !!
1037      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
1038      !!----------------------------------------------------------------------
1039      REAL(wp)                        , INTENT(in   ) ::   pamsk     ! advection of concentration (1) or other tracers (0)
1040      INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
1041      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step
1042      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pv        ! ice j-velocity component
1043      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt        ! tracer fields
1044      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_v      ! tracer at v-point
1045      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfv_ho    ! high order flux
1046      !
1047      INTEGER  ::   ji, jj, jl         ! dummy loop indices
1048      REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
1049      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztv1, ztv2, ztv3, ztv4
1050      !!----------------------------------------------------------------------
1051      !
1052      !                                                     !--  Laplacian in j-direction  --!
1053      DO jl = 1, jpl
1054         DO jj = 1, jpjm1         ! First derivative (gradient)
1055            DO ji = fs_2, fs_jpim1
1056               ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
1057            END DO
1058         END DO
1059         DO jj = 2, jpjm1         ! Second derivative (Laplacian)
1060            DO ji = fs_2, fs_jpim1
1061               ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj)
1062            END DO
1063         END DO
1064      END DO
1065      CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1. )
1066      !
1067      !                                                     !--  BiLaplacian in j-direction  --!
1068      DO jl = 1, jpl
1069         DO jj = 1, jpjm1         ! First derivative
1070            DO ji = fs_2, fs_jpim1
1071               ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
1072            END DO
1073         END DO
1074         DO jj = 2, jpjm1         ! Second derivative
1075            DO ji = fs_2, fs_jpim1
1076               ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj)
1077            END DO
1078         END DO
1079      END DO
1080      CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1. )
1081      !
1082      !
1083      SELECT CASE (kn_umx )
1084         !
1085      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
1086         DO jl = 1, jpl
1087            DO jj = 1, jpjm1
1088               DO ji = 1, fs_jpim1
1089                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   &
1090                     &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
1091               END DO
1092            END DO
1093         END DO
1094         !
1095      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
1096         DO jl = 1, jpl
1097            DO jj = 1, jpjm1
1098               DO ji = 1, fs_jpim1
1099                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1100                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   &
1101                     &                                                            - zcv *   ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
1102               END DO
1103            END DO
1104         END DO
1105         !
1106      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
1107         DO jl = 1, jpl
1108            DO jj = 1, jpjm1
1109               DO ji = 1, fs_jpim1
1110                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1111                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
1112!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj)
1113                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
1114                     &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
1115                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
1116                     &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
1117               END DO
1118            END DO
1119         END DO
1120         !
1121      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
1122         DO jl = 1, jpl
1123            DO jj = 1, jpjm1
1124               DO ji = 1, fs_jpim1
1125                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1126                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
1127!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj)
1128                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
1129                     &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
1130                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
1131                     &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
1132               END DO
1133            END DO
1134         END DO
1135         !
1136      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
1137         DO jl = 1, jpl
1138            DO jj = 1, jpjm1
1139               DO ji = 1, fs_jpim1
1140                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1141                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
1142!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj)
1143                  zdy4 = zdy2 * zdy2
1144                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
1145                     &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
1146                     &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) * (                       ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
1147                     &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) &
1148                     &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl)     &
1149                     &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
1150               END DO
1151            END DO
1152         END DO
1153         !
1154      END SELECT
1155      !
1156      ! if pt at v-point is negative then use the upstream value
1157      !    this should not be necessary if a proper sea-ice mask is set in Ultimate
1158      !    to degrade the order of the scheme when necessary (for ex. at the ice edge)
1159      IF( ll_neg ) THEN
1160         DO jl = 1, jpl
1161            DO jj = 1, jpjm1
1162               DO ji = 1, fs_jpim1
1163                  IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
1164                     pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
1165                        &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
1166                  ENDIF
1167               END DO
1168            END DO
1169         END DO
1170      ENDIF
1171      !                                                     !-- High order flux in j-direction  --!
1172      DO jl = 1, jpl
1173         DO jj = 1, jpjm1
1174            DO ji = 1, fs_jpim1   ! vector opt.
1175               pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
1176            END DO
1177         END DO
1178      END DO
1179      !
1180   END SUBROUTINE ultimate_y
1181     
1182
1183   SUBROUTINE nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
1184      !!---------------------------------------------------------------------
1185      !!                    ***  ROUTINE nonosc_ice  ***
1186      !!     
1187      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
1188      !!       scheme and the before field by a non-oscillatory algorithm
1189      !!
1190      !! **  Method  :   ...
1191      !!----------------------------------------------------------------------
1192      REAL(wp)                   , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
1193      REAL(wp)                   , INTENT(in   ) ::   pdt              ! tracer time-step
1194      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pu               ! ice i-velocity => u*e2
1195      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pv               ! ice j-velocity => v*e1
1196      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pt, pt_ups       ! before field & upstream guess of after field
1197      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pfv_ups, pfu_ups ! upstream flux
1198      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux
1199      !
1200      INTEGER  ::   ji, jj, jl    ! dummy loop indices
1201      REAL(wp) ::   zpos, zneg, zbig, zup, zdo, z1_dt              ! local scalars
1202      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zcoef, zzt       !   -      -
1203      REAL(wp), DIMENSION(jpi,jpj    ) :: zbup, zbdo
1204      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_ups, ztj_ups
1205      !!----------------------------------------------------------------------
1206      zbig = 1.e+40_wp
1207     
1208      ! antidiffusive flux : high order minus low order
1209      ! --------------------------------------------------
1210      DO jl = 1, jpl
1211         DO jj = 1, jpjm1
1212            DO ji = 1, fs_jpim1   ! vector opt.
1213               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
1214               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
1215            END DO
1216         END DO
1217      END DO
1218
1219      ! extreme case where pfu_ho has to be zero
1220      ! ----------------------------------------
1221      !                                    pfu_ho
1222      !                           *         --->
1223      !                        |      |  *   |        |
1224      !                        |      |      |    *   |   
1225      !                        |      |      |        |    *
1226      !            t_ups :       i-1     i       i+1       i+2   
1227      IF( ll_prelim ) THEN
1228         
1229         DO jl = 1, jpl
1230            DO jj = 2, jpjm1
1231               DO ji = fs_2, fs_jpim1 
1232                  zti_ups(ji,jj,jl)= pt_ups(ji+1,jj  ,jl)
1233                  ztj_ups(ji,jj,jl)= pt_ups(ji  ,jj+1,jl)
1234               END DO
1235            END DO
1236         END DO
1237         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. )
1238
1239         DO jl = 1, jpl
1240            DO jj = 2, jpjm1
1241               DO ji = fs_2, fs_jpim1
1242                  IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj  ,jl) - pt_ups(ji,jj,jl) ) <= 0._wp .AND.  &
1243                     & pfv_ho(ji,jj,jl) * ( pt_ups(ji  ,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0._wp ) THEN
1244                     !
1245                     IF(  pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj  ,jl) - zti_ups(ji,jj,jl) ) <= 0._wp .AND.  &
1246                        & pfv_ho(ji,jj,jl) * ( ztj_ups(ji  ,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0._wp ) THEN
1247                        pfu_ho(ji,jj,jl)=0._wp
1248                        pfv_ho(ji,jj,jl)=0._wp
1249                     ENDIF
1250                     !
1251                     IF(  pfu_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji-1,jj  ,jl) ) <= 0._wp .AND.  &
1252                        & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji  ,jj-1,jl) ) <= 0._wp ) THEN
1253                        pfu_ho(ji,jj,jl)=0._wp
1254                        pfv_ho(ji,jj,jl)=0._wp
1255                     ENDIF
1256                     !
1257                  ENDIF
1258               END DO
1259            END DO
1260         END DO
1261         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
1262
1263      ENDIF
1264
1265      ! Search local extrema
1266      ! --------------------
1267      ! max/min of pt & pt_ups with large negative/positive value (-/+zbig) outside ice cover
1268      z1_dt = 1._wp / pdt
1269      DO jl = 1, jpl
1270         
1271         DO jj = 1, jpj
1272            DO ji = 1, jpi
1273               IF    ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
1274                  zbup(ji,jj) = -zbig
1275                  zbdo(ji,jj) =  zbig
1276               ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN
1277                  zbup(ji,jj) = pt_ups(ji,jj,jl)
1278                  zbdo(ji,jj) = pt_ups(ji,jj,jl)
1279               ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
1280                  zbup(ji,jj) = pt(ji,jj,jl)
1281                  zbdo(ji,jj) = pt(ji,jj,jl)
1282               ELSE
1283                  zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
1284                  zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
1285               ENDIF
1286            END DO
1287         END DO
1288
1289         DO jj = 2, jpjm1
1290            DO ji = fs_2, fs_jpim1   ! vector opt.
1291               !
1292               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
1293               zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) )
1294               !
1295               zpos = MAX( 0._wp, pfu_ho(ji-1,jj  ,jl) ) - MIN( 0._wp, pfu_ho(ji  ,jj  ,jl) ) &  ! positive/negative part of the flux
1296                  & + MAX( 0._wp, pfv_ho(ji  ,jj-1,jl) ) - MIN( 0._wp, pfv_ho(ji  ,jj  ,jl) )
1297               zneg = MAX( 0._wp, pfu_ho(ji  ,jj  ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj  ,jl) ) &
1298                  & + MAX( 0._wp, pfv_ho(ji  ,jj  ,jl) ) - MIN( 0._wp, pfv_ho(ji  ,jj-1,jl) )
1299               !
1300               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) ) &
1301                  &          ) * ( 1. - pamsk )
1302               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) ) &
1303                  &          ) * ( 1. - pamsk )
1304               !
1305               !                                  ! up & down beta terms
1306               ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!)
1307               IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt
1308               ELSE                     ; zbetup(ji,jj,jl) = 0._wp ! zbig
1309               ENDIF
1310               !
1311               IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
1312               ELSE                     ; zbetdo(ji,jj,jl) = 0._wp ! zbig
1313               ENDIF
1314               !
1315               ! if all the points are outside ice cover
1316               IF( zup == -zbig )   zbetup(ji,jj,jl) = 0._wp ! zbig
1317               IF( zdo ==  zbig )   zbetdo(ji,jj,jl) = 0._wp ! zbig           
1318               !
1319            END DO
1320         END DO
1321      END DO
1322      CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
1323
1324     
1325      ! monotonic flux in the y direction
1326      ! ---------------------------------
1327      DO jl = 1, jpl
1328         DO jj = 1, jpjm1
1329            DO ji = 1, fs_jpim1   ! vector opt.
1330               zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) )
1331               zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) )
1332               zcu = 0.5_wp + SIGN( 0.5_wp , pfu_ho(ji,jj,jl) )
1333               !
1334               zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
1335               !
1336               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl)
1337               !
1338            END DO
1339         END DO
1340
1341         DO jj = 1, jpjm1
1342            DO ji = 1, fs_jpim1   ! vector opt.
1343               zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) )
1344               zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) )
1345               zcv = 0.5_wp + SIGN( 0.5_wp , pfv_ho(ji,jj,jl) )
1346               !
1347               zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
1348               !
1349               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl)
1350               !
1351            END DO
1352         END DO
1353
1354      END DO
1355      !
1356   END SUBROUTINE nonosc_ice
1357
1358   
1359   SUBROUTINE limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
1360      !!---------------------------------------------------------------------
1361      !!                    ***  ROUTINE limiter_x  ***
1362      !!     
1363      !! **  Purpose :   compute flux limiter
1364      !!----------------------------------------------------------------------
1365      REAL(wp)                  , INTENT(in   ) ::   pdt          ! tracer time-step
1366      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu           ! ice i-velocity => u*e2
1367      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt           ! ice tracer
1368      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pfu_ups      ! upstream flux
1369      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pfu_ho       ! high order flux
1370      !
1371      REAL(wp) ::   Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr
1372      INTEGER  ::   ji, jj, jl    ! dummy loop indices
1373      REAL(wp), DIMENSION (jpi,jpj,jpl) ::   zslpx       ! tracer slopes
1374      !!----------------------------------------------------------------------
1375      !
1376      DO jl = 1, jpl
1377         DO jj = 2, jpjm1
1378            DO ji = fs_2, fs_jpim1   ! vector opt.
1379               zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1)
1380            END DO
1381         END DO
1382      END DO
1383      CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.)   ! lateral boundary cond.
1384     
1385      DO jl = 1, jpl
1386         DO jj = 2, jpjm1
1387            DO ji = fs_2, fs_jpim1   ! vector opt.
1388               uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
1389               
1390               Rjm = zslpx(ji-1,jj,jl)
1391               Rj  = zslpx(ji  ,jj,jl)
1392               Rjp = zslpx(ji+1,jj,jl)
1393
1394               IF( np_limiter == 3 ) THEN
1395
1396                  IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm
1397                  ELSE                        ;   Rr = Rjp
1398                  ENDIF
1399
1400                  zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)     
1401                  IF( Rj > 0. ) THEN
1402                     zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)),  &
1403                        &        MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
1404                  ELSE
1405                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)),  &
1406                        &        MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
1407                  ENDIF
1408                  pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter
1409
1410               ELSEIF( np_limiter == 2 ) THEN
1411                  IF( Rj /= 0. ) THEN
1412                     IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
1413                     ELSE                        ;   Cr = Rjp / Rj
1414                     ENDIF
1415                  ELSE
1416                     Cr = 0.
1417                  ENDIF
1418
1419                  ! -- superbee --
1420                  zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1421                  ! -- van albada 2 --
1422                  !!zpsi = 2.*Cr / (Cr*Cr+1.)
1423                  ! -- sweby (with beta=1) --
1424                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1425                  ! -- van Leer --
1426                  !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1427                  ! -- ospre --
1428                  !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
1429                  ! -- koren --
1430                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
1431                  ! -- charm --
1432                  !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
1433                  !ELSE                 ;   zpsi = 0.
1434                  !ENDIF
1435                  ! -- van albada 1 --
1436                  !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
1437                  ! -- smart --
1438                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
1439                  ! -- umist --
1440                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1441
1442                  ! high order flux corrected by the limiter
1443                  pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5
1444
1445               ENDIF
1446            END DO
1447         END DO
1448      END DO
1449      CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.)   ! lateral boundary cond.
1450      !
1451   END SUBROUTINE limiter_x
1452
1453   
1454   SUBROUTINE limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
1455      !!---------------------------------------------------------------------
1456      !!                    ***  ROUTINE limiter_y  ***
1457      !!     
1458      !! **  Purpose :   compute flux limiter
1459      !!----------------------------------------------------------------------
1460      REAL(wp)                   , INTENT(in   ) ::   pdt          ! tracer time-step
1461      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pv           ! ice i-velocity => u*e2
1462      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pt           ! ice tracer
1463      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pfv_ups      ! upstream flux
1464      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ho       ! high order flux
1465      !
1466      REAL(wp) ::   Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr
1467      INTEGER  ::   ji, jj, jl    ! dummy loop indices
1468      REAL(wp), DIMENSION (jpi,jpj,jpl) ::   zslpy       ! tracer slopes
1469      !!----------------------------------------------------------------------
1470      !
1471      DO jl = 1, jpl
1472         DO jj = 2, jpjm1
1473            DO ji = fs_2, fs_jpim1   ! vector opt.
1474               zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1)
1475            END DO
1476         END DO
1477      END DO
1478      CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.)   ! lateral boundary cond.
1479
1480      DO jl = 1, jpl
1481         DO jj = 2, jpjm1
1482            DO ji = fs_2, fs_jpim1   ! vector opt.
1483               vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
1484
1485               Rjm = zslpy(ji,jj-1,jl)
1486               Rj  = zslpy(ji,jj  ,jl)
1487               Rjp = zslpy(ji,jj+1,jl)
1488
1489               IF( np_limiter == 3 ) THEN
1490
1491                  IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm
1492                  ELSE                        ;   Rr = Rjp
1493                  ENDIF
1494
1495                  zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)     
1496                  IF( Rj > 0. ) THEN
1497                     zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)),  &
1498                        &        MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
1499                  ELSE
1500                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)),  &
1501                        &        MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
1502                  ENDIF
1503                  pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter
1504
1505               ELSEIF( np_limiter == 2 ) THEN
1506
1507                  IF( Rj /= 0. ) THEN
1508                     IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
1509                     ELSE                        ;   Cr = Rjp / Rj
1510                     ENDIF
1511                  ELSE
1512                     Cr = 0.
1513                  ENDIF
1514
1515                  ! -- superbee --
1516                  zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1517                  ! -- van albada 2 --
1518                  !!zpsi = 2.*Cr / (Cr*Cr+1.)
1519                  ! -- sweby (with beta=1) --
1520                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1521                  ! -- van Leer --
1522                  !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1523                  ! -- ospre --
1524                  !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
1525                  ! -- koren --
1526                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
1527                  ! -- charm --
1528                  !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
1529                  !ELSE                 ;   zpsi = 0.
1530                  !ENDIF
1531                  ! -- van albada 1 --
1532                  !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
1533                  ! -- smart --
1534                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
1535                  ! -- umist --
1536                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1537
1538                  ! high order flux corrected by the limiter
1539                  pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5
1540
1541               ENDIF
1542            END DO
1543         END DO
1544      END DO
1545      CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.)   ! lateral boundary cond.
1546      !
1547   END SUBROUTINE limiter_y
1548
1549
1550   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
1551      !!-------------------------------------------------------------------
1552      !!                  ***  ROUTINE Hbig  ***
1553      !!
1554      !! ** Purpose : Thickness correction in case advection scheme creates
1555      !!              abnormally thick ice or snow
1556      !!
1557      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points
1558      !!                 (before advection) and reduce it by adapting ice concentration
1559      !!              2- check whether snow thickness is larger than the surrounding 9-points
1560      !!                 (before advection) and reduce it by sending the excess in the ocean
1561      !!              3- check whether snow load deplets the snow-ice interface below sea level$
1562      !!                 and reduce it by sending the excess in the ocean
1563      !!              4- correct pond fraction to avoid a_ip > a_i
1564      !!
1565      !! ** input   : Max thickness of the surrounding 9-points
1566      !!-------------------------------------------------------------------
1567      REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step
1568      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts
1569      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip
1570      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
1571      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i
1572      !
1573      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
1574      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zvs_excess, zfra
1575      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
1576      !!-------------------------------------------------------------------
1577      !
1578      z1_dt = 1._wp / pdt
1579      !
1580      DO jl = 1, jpl
1581
1582         DO jj = 1, jpj
1583            DO ji = 1, jpi
1584               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
1585                  !
1586                  !                               ! -- check h_ip -- !
1587                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
1588                  IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
1589                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
1590                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
1591                        pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl)
1592                     ENDIF
1593                  ENDIF
1594                  !
1595                  !                               ! -- check h_i -- !
1596                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
1597                  zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
1598                  IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
1599                     pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
1600                  ENDIF
1601                  !
1602                  !                               ! -- check h_s -- !
1603                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
1604                  zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
1605                  IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
1606                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
1607                     !
1608                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt
1609                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
1610                     !
1611                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
1612                     pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
1613                  ENDIF           
1614                  !
1615                  !                               ! -- check snow load -- !
1616                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean
1617                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)
1618                  !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically)
1619                  zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
1620                  IF( zvs_excess > 0._wp ) THEN
1621                     zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
1622                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
1623                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
1624                     !
1625                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
1626                     pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess
1627                  ENDIF
1628                 
1629               ENDIF
1630            END DO
1631         END DO
1632      END DO 
1633      !                                           !-- correct pond fraction to avoid a_ip > a_i
1634      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:)
1635      !
1636      !
1637   END SUBROUTINE Hbig
1638   
1639#else
1640   !!----------------------------------------------------------------------
1641   !!   Default option           Dummy module         NO SI3 sea-ice model
1642   !!----------------------------------------------------------------------
1643#endif
1644
1645   !!======================================================================
1646END MODULE icedyn_adv_umx
Note: See TracBrowser for help on using the repository browser.