source: NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_GPU/src/ICE/icedyn_adv_umx.F90 @ 13153

Last change on this file since 13153 was 13153, checked in by andmirek, 15 months ago

Ticket #2482: fixes after testing with ifort

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