source: NEMO/trunk/src/ICE/icedyn_adv_umx.F90 @ 11612

Last change on this file since 11612 was 11612, checked in by clem, 16 months ago

make Prather the default advection scheme for ice simulations and drastically reduce the number of communications in Prather routine (by a factor of 10)

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