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

Last change on this file since 13633 was 13633, checked in by clem, 7 months ago

trunk: repro issue(?) for UMx advection scheme. One loop index was incorrect.

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