source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE/icedyn_adv_umx.F90

Last change on this file was 12546, checked in by orioltp, 13 months ago

Adding precision specification in hardcoded reals and other modifications to allow compilation without forcing reals without precision specification to a certain value through compiler flags

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