New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_adv_umx.F90 in NEMO/trunk/src/ICE – NEMO

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

Last change on this file since 11229 was 10945, checked in by clem, 5 years ago

finish repairing the code when only ice dynamics is activated

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