source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_umx.F90 @ 12252

Last change on this file since 12252 was 12252, checked in by smasson, 11 months ago

rev12240_dev_r11943_MERGE_2019: same as [12251], merge trunk 12072:12248, all sette tests ok, GYRE_PISCES, AMM12, ISOMIP, VORTEX intentical to 12236

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