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

Last change on this file since 12379 was 12379, checked in by dancopsey, 15 months ago

Add meltpond lid thickness as a new prognostic.

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