source: NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90 @ 10910

Last change on this file since 10910 was 10910, checked in by clem, 20 months ago

Major change: the advection scheme UMx has been revisited to clean all the unphysical values which occured. Minor change: the ref config SPITZ12 has been slightly modified to test more parameterizations that are available in the code (melt ponds and landfast)

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