source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ICE/icedyn_adv_umx.F90 @ 12808

Last change on this file since 12808 was 12489, checked in by davestorkey, 9 months ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp → rn_atfp (namelist parameter used everywhere) and rau0 → rho0.

This version gives bit-comparable results to the previous version of the trunk.

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