source: NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_GPU/src/ICE/icedyn_adv_umx.F90 @ 13149

Last change on this file since 13149 was 13149, checked in by andmirek, 15 months ago

Ticket #2482: Dissable restart and use allocatable arrays

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