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/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE – NEMO

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icedyn_adv_umx.F90 @ 13553

Last change on this file since 13553 was 13553, checked in by hadcv, 3 years ago

Merge in trunk up to [13550]

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