source: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_adv_umx.F90 @ 13466

Last change on this file since 13466 was 13466, checked in by smasson, 7 months ago

r4_trunk: merge r4 13280:13310, see #2523

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