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/2018/dev_r9947_SI3_advection/src/ICE – NEMO

source: NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90 @ 10267

Last change on this file since 10267 was 10267, checked in by clem, 6 years ago

improve ice thickness after advection (advection of H using u*a from advection of concentration)
File size: 53.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 trend with the 3D advection trends using a TVD scheme
14   !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders
15   !!   macho             : ???
16   !!   nonosc            : compute monotonic tracer fluxes by 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 lib_mpp        ! MPP library
26   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
27   USE lbclnk         ! lateral boundary conditions (or mpp links)
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   ice_dyn_adv_umx   ! called by icedyn_adv.F90
33     
34   REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6
35   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
36
37   !! * Substitutions
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
41   !! $Id: icedyn_adv_umx.F90 4499 2014-02-18 15:14:31Z timgraham $
42   !! Software governed by the CeCILL licence     (./LICENSE)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice,  &
47      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE ice_dyn_adv_umx  ***
50      !!
51      !! **  Purpose :   Compute the now trend due to total advection of
52      !!                 tracers and add it to the general trend of tracer equations
53      !!                 using an "Ultimate-Macho" scheme
54      !!
55      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
56      !!----------------------------------------------------------------------
57      INTEGER                     , INTENT(in   ) ::   kn_umx     ! order of the scheme (1-5=UM or 20=CEN2)
58      INTEGER                     , INTENT(in   ) ::   kt         ! time step
59      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
60      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
61      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
62      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
63      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
64      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
65      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
66      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
67      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
68      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
69      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
70      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
71      !
72      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
73      INTEGER  ::   icycle                  ! number of sub-timestep for the advection
74      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers
75      REAL(wp) ::   zcfl , zdt
76      REAL(wp) ::   zeps = 0.1_wp           ! shift in concentration to avoid division by 0
77      !                                     !    must be >= 0.01 and the best seems to be 0.1
78      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx, zcu_box, zcv_box, zfu, zfv
79      REAL(wp), DIMENSION(jpi,jpj) ::   z1_a, zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, z1_ap, zh_ip 
80      !!----------------------------------------------------------------------
81      !
82      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
83      !
84      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !       
85      zcfl  =            MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
86      zcfl  = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
87      IF( lk_mpp )   CALL mpp_max( zcfl )
88
89      IF( zcfl > 0.5 ) THEN   ;   icycle = 2 
90      ELSE                    ;   icycle = 1 
91      ENDIF
92
93      zdt = rdt_ice / REAL(icycle)
94
95      ! --- transport --- !
96      zudy(:,:) = pu_ice(:,:) * e2u(:,:)
97      zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
98
99      ! --- define velocity for advection: u*grad(H) --- !
100      DO jj = 2, jpjm1
101         DO ji = fs_2, fs_jpim1
102            IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
103            ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj)
104            ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj)
105            ENDIF
106
107            IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
108            ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1)
109            ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  )
110            ENDIF
111         END DO
112      END DO
113
114      !---------------!
115      !== advection ==!
116      !---------------!
117      DO jt = 1, icycle
118
119         zamsk = 1._wp
120         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) )   ! Open water area
121
122         DO jl = 1, jpl
123            ! to avoid a problem with the limiter nonosc when A gets close to 0
124            pa_i(:,:,jl) = pa_i(:,:,jl) + zeps * tmask(:,:,1)
125            !
126            WHERE( pa_i(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_i(:,:,jl)
127            ELSEWHERE                        ;   z1_a(:,:) = 0.   !!;   pa_i(:,:,jl) = 0. ; pv_i(:,:,jl) = 0.
128            END WHERE
129            !
130            zamsk = 1._wp
131            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zfu, zfv ) ! Ice area
132            !!zfu = zudy ; zfv = zvdx
133            !!CALL adv_umx( kn_umx, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, pv_i(:,:,jl), pv_i(:,:,jl) )
134            zamsk = 0._wp ; zh_i(:,:) = pv_i (:,:,jl) * z1_a(:,:)
135            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_i(:,:), pv_i (:,:,jl) )              ! Ice volume
136            zamsk = 0._wp ; zh_s(:,:) = pv_s (:,:,jl) * z1_a(:,:)
137            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_s(:,:), pv_s (:,:,jl) )              ! Snw volume
138            zamsk = 0._wp ; zs_i(:,:) = psv_i(:,:,jl) * z1_a(:,:)
139            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zs_i(:,:), psv_i(:,:,jl) )              ! Salt content
140            zamsk = 0._wp ; zo_i(:,:) = poa_i(:,:,jl) * z1_a(:,:)
141            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zo_i(:,:), poa_i(:,:,jl) )              ! Age content
142            DO jk = 1, nlay_i
143               zamsk = 0._wp ; ze_i(:,:) = pe_i(:,:,jk,jl) * z1_a(:,:)
144               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu, zfv, zcu_box, zcv_box, ze_i(:,:), pe_i(:,:,jk,jl) )           ! Ice heat content
145            END DO
146            DO jk = 1, nlay_s
147               zamsk = 0._wp ; ze_s(:,:) = pe_s(:,:,jk,jl) * z1_a(:,:)
148               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu, zfv, zcu_box, zcv_box, ze_s(:,:), pe_s(:,:,jk,jl) )           ! Snw heat content
149            END DO
150            !
151            IF ( ln_pnd_H12 ) THEN   ! melt ponds
152               ! to avoid a problem with the limiter nonosc when A gets close to 0
153               pa_ip(:,:,jl) = pa_ip(:,:,jl) + zeps * tmask(:,:,1)
154               !
155               WHERE( pa_ip(:,:,jl) > epsi20 )   ;   z1_ap(:,:) = 1._wp / pa_ip(:,:,jl)
156               ELSEWHERE                         ;   z1_ap(:,:) = 0.
157               END WHERE
158               !
159               zamsk = 1._wp
160               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zfu, zfv ) ! mp fraction
161               zamsk = 0._wp ; zh_ip(:,:) = pv_ip(:,:,jl) * z1_a(:,:)
162               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_ip(:,:), pv_ip(:,:,jl) )              ! mp volume
163            ENDIF
164            !
165            ! to avoid a problem with the limiter nonosc when A gets close to 0
166            DO jj = 2, jpjm1
167               DO ji = fs_2, fs_jpim1
168                  !pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps ) * tmask(ji,jj,1)
169                  pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps &
170                     &             + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt &
171                     &             ) * tmask(ji,jj,1)
172                  IF ( ln_pnd_H12 ) THEN   ! melt ponds
173                     pa_ip(ji,jj,jl) = ( pa_ip(ji,jj,jl) - zeps &
174                        &              + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt &
175                        &              ) * tmask(ji,jj,1)
176                  ENDIF
177              END DO
178            END DO
179            CALL lbc_lnk( pa_i(:,:,jl), 'T',  1. )
180            !
181         END DO
182
183      END DO
184      !
185   END SUBROUTINE ice_dyn_adv_umx
186
187   
188   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pfu, pfv )
189      !!----------------------------------------------------------------------
190      !!                  ***  ROUTINE adv_umx  ***
191      !!
192      !! **  Purpose :   Compute the now trend due to total advection of
193      !!       tracers and add it to the general trend of tracer equations
194      !!
195      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
196      !!       corrected flux (monotonic correction)
197      !!       note: - this advection scheme needs a leap-frog time scheme
198      !!
199      !! ** Action : - pt  the after advective tracer
200      !!----------------------------------------------------------------------
201      REAL(wp)                    , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0)
202      INTEGER                     , INTENT(in   )           ::   kn_umx         ! order of the scheme (1-5=UM or 20=CEN2)
203      INTEGER                     , INTENT(in   )           ::   jt             ! number of sub-iteration
204      INTEGER                     , INTENT(in   )           ::   kt             ! number of iteration
205      REAL(wp)                    , INTENT(in   )           ::   pdt            ! tracer time-step
206      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2
207      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u
208      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity
209      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pt             ! tracer field
210      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptc            ! tracer content field
211      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pfu, pfv       ! high order fluxes
212      !
213      INTEGER  ::   ji, jj           ! dummy loop indices 
214      REAL(wp) ::   ztra             ! local scalar
215!!clem      REAL(wp) ::   zeps = 1.e-02        ! local scalar
216      INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid)
217      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt
218      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups, zt_ups   ! only for nonosc
219      !!----------------------------------------------------------------------
220      !
221      ! add a constant value to avoid problems with zeros
222      DO jj = 1, jpj
223         DO ji = 1, jpi
224            zpt(ji,jj) = pt(ji,jj) !!clem + zeps * tmask(ji,jj,1)
225         END DO
226      END DO
227
228      !  upstream (_ups) advection with initial mass fluxes
229      ! ---------------------------------------------------
230      ! fluxes
231      DO jj = 1, jpjm1
232         DO ji = 1, fs_jpim1
233            zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj)
234            zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1)
235         END DO
236      END DO
237      ! guess after content field with upstream scheme
238      DO jj = 2, jpjm1
239         DO ji = fs_2, fs_jpim1
240            ztra          = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   &
241               &                + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1) ) * r1_e1e2t(ji,jj)
242            zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
243         END DO
244      END DO
245      CALL lbc_lnk( zt_ups, 'T', 1. )
246     
247      ! High order (_ho) fluxes
248      ! -----------------------
249      SELECT CASE( kn_umx )
250         !
251      CASE ( 20 )                          !== centered second order ==!
252         !
253         CALL cen2( kn_limiter, jt, kt, pdt, zpt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  &
254            &       zt_ups, zfu_ups, zfv_ups )
255         !
256      CASE ( 1:5 )                         !== 1st to 5th order ULTIMATE-MACHO scheme ==!
257         !
258         CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, zpt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho,  &
259            &        zt_ups, zfu_ups, zfv_ups )
260         !
261      END SELECT
262         
263      ! output high order fluxes
264      ! ------------------------
265      IF( PRESENT(pfu) ) THEN
266         DO jj = 1, jpjm1
267            DO ji = 1, fs_jpim1
268               pfu(ji,jj) =  zfu_ho(ji,jj) !!clem - zeps * puc(ji,jj)
269               pfv(ji,jj) =  zfv_ho(ji,jj) !!clem - zeps * pvc(ji,jj)
270            END DO
271         END DO
272         !!CALL lbc_lnk( pfu, 'U',  -1. ) ! clem: not needed I think
273         !!CALL lbc_lnk( pfv, 'V',  -1. )
274      ENDIF
275     
276      ! final trend with corrected fluxes
277      ! ------------------------------------
278      DO jj = 2, jpjm1
279         DO ji = fs_2, fs_jpim1 
280            ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )         &  !        Div(uaH) or        Div(ua)
281!!clem               &           + ( puc   (ji,jj) - puc   (ji-1,jj) + pvc   (ji,jj) - pvc   (ji,jj-1) ) * zeps  &  ! epsi * Div(ua)  or epsi * Div(u)
282               &         ) * r1_e1e2t(ji,jj) 
283            ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
284         END DO
285      END DO
286      CALL lbc_lnk( ptc, 'T',  1. )
287      !
288   END SUBROUTINE adv_umx
289
290   SUBROUTINE cen2( kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, &
291      &             pt_ups, pfu_ups, pfv_ups )
292      !!---------------------------------------------------------------------
293      !!                    ***  ROUTINE macho  ***
294      !!     
295      !! **  Purpose :   compute 
296      !!
297      !! **  Method  :   ... ???
298      !!                 TIM = transient interpolation Modeling
299      !!
300      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
301      !!----------------------------------------------------------------------
302      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter
303      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration
304      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration
305      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step
306      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt               ! tracer fields
307      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
308      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components
309      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc              ! tracer content at before time step
310      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes
311      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content
312      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes
313      !
314      INTEGER  ::   ji, jj    ! dummy loop indices
315      LOGICAL  ::   ll_xy = .TRUE.
316      REAL(wp), DIMENSION(jpi,jpj) ::   zzt
317      !!----------------------------------------------------------------------
318      !
319      IF( .NOT.ll_xy ) THEN   !-- no alternate directions --!
320         !
321         DO jj = 1, jpjm1
322            DO ji = 1, fs_jpim1
323               pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
324               pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
325            END DO
326         END DO
327         IF    ( kn_limiter == 1 ) THEN
328            CALL nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
329         ELSEIF( kn_limiter == 2 ) THEN
330            CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
331            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
332         ELSEIF( kn_limiter == 3 ) THEN
333            CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
334            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
335         ENDIF
336         !
337      ELSE                    !-- alternate directions --!
338         !
339         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
340            !
341            ! flux in x-direction
342            DO jj = 1, jpjm1
343               DO ji = 1, fs_jpim1
344                  pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
345               END DO
346            END DO
347            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
348            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
349
350            ! first guess of tracer content from u-flux
351            DO jj = 2, jpjm1
352               DO ji = fs_2, fs_jpim1   ! vector opt.
353                  zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
354               END DO
355            END DO
356            CALL lbc_lnk( zzt, 'T', 1. )
357
358            ! flux in y-direction
359            DO jj = 1, jpjm1
360               DO ji = 1, fs_jpim1
361                  pfv_ho(ji,jj) = 0.5 * pv(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) )
362               END DO
363            END DO
364            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
365            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
366
367         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
368            !
369            ! flux in y-direction
370            DO jj = 1, jpjm1
371               DO ji = 1, fs_jpim1
372                  pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
373               END DO
374            END DO
375            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
376            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
377            !
378            ! first guess of tracer content from v-flux
379            DO jj = 2, jpjm1
380               DO ji = fs_2, fs_jpim1   ! vector opt.
381                  zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
382               END DO
383            END DO
384            CALL lbc_lnk( zzt, 'T', 1. )
385            !
386            ! flux in x-direction
387            DO jj = 1, jpjm1
388               DO ji = 1, fs_jpim1
389                  pfu_ho(ji,jj) = 0.5 * pu(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) )
390               END DO
391            END DO
392            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
393            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
394
395         ENDIF
396         IF( kn_limiter == 1 )   CALL nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
397         
398      ENDIF
399   
400   END SUBROUTINE cen2
401
402   
403   SUBROUTINE macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, pt_u, pt_v, pfu_ho, pfv_ho, &
404      &              pt_ups, pfu_ups, pfv_ups )
405      !!---------------------------------------------------------------------
406      !!                    ***  ROUTINE macho  ***
407      !!     
408      !! **  Purpose :   compute 
409      !!
410      !! **  Method  :   ... ???
411      !!                 TIM = transient interpolation Modeling
412      !!
413      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
414      !!----------------------------------------------------------------------
415      REAL(wp)                    , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
416      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter
417      INTEGER                     , INTENT(in   ) ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2)
418      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration
419      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration
420      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step
421      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt               ! tracer fields
422      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
423      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components
424      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox     ! upstream velocity
425      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc              ! tracer content at before time step
426      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u, pt_v       ! tracer at u- and v-points
427      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes
428      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content
429      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes
430      !
431      INTEGER  ::   ji, jj    ! dummy loop indices
432      REAL(wp), DIMENSION(jpi,jpj) ::   zzt
433      !!----------------------------------------------------------------------
434      !
435      IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
436         !
437         !                                                        !--  ultimate interpolation of pt at u-point  --!
438         CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
439         !                                                        !--  limiter in x --!
440         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
441         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
442         !                                                        !--  advective form update in zzt  --!
443         DO jj = 2, jpjm1
444            DO ji = fs_2, fs_jpim1   ! vector opt.
445               zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  &
446                  &                   - pt   (ji,jj) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk
447               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
448            END DO
449         END DO
450         CALL lbc_lnk( zzt, 'T', 1. )
451         !                                                        !--  ultimate interpolation of pt at v-point  --!
452         CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho )
453         !                                                        !--  limiter in y --!
454         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
455         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
456         !
457      ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
458         !
459         !                                                        !--  ultimate interpolation of pt at v-point  --!
460         CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
461         !                                                        !--  limiter in y --!
462         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
463         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
464         !                                                        !--  advective form update in zzt  --!
465         DO jj = 2, jpjm1
466            DO ji = fs_2, fs_jpim1
467               zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  &
468                  &                    - pt  (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk
469               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
470            END DO
471         END DO
472         CALL lbc_lnk( zzt, 'T', 1. )
473         !                                                        !--  ultimate interpolation of pt at u-point  --!
474         CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho )
475         !                                                        !--  limiter in x --!
476         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
477         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
478         !
479      ENDIF
480      IF( kn_limiter == 1 )   CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
481      !
482   END SUBROUTINE macho
483
484
485   SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
486      !!---------------------------------------------------------------------
487      !!                    ***  ROUTINE ultimate_x  ***
488      !!     
489      !! **  Purpose :   compute 
490      !!
491      !! **  Method  :   ... ???
492      !!                 TIM = transient interpolation Modeling
493      !!
494      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
495      !!----------------------------------------------------------------------
496      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
497      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
498      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu        ! ice i-velocity component
499      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc       ! ice i-velocity * A component
500      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
501      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u      ! tracer at u-point
502      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho    ! high order flux
503      !
504      INTEGER  ::   ji, jj             ! dummy loop indices
505      REAL(wp) ::   zcu, zdx2, zdx4    !   -      -
506      REAL(wp), DIMENSION(jpi,jpj) ::   ztu1, ztu2, ztu3, ztu4
507      !!----------------------------------------------------------------------
508      !
509      !                                                     !--  Laplacian in i-direction  --!
510      DO jj = 2, jpjm1         ! First derivative (gradient)
511         DO ji = 1, fs_jpim1
512            ztu1(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
513         END DO
514         !                     ! Second derivative (Laplacian)
515         DO ji = fs_2, fs_jpim1
516            ztu2(ji,jj) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj)
517         END DO
518      END DO
519      CALL lbc_lnk( ztu2, 'T', 1. )
520      !
521      !                                                     !--  BiLaplacian in i-direction  --!
522      DO jj = 2, jpjm1         ! Third derivative
523         DO ji = 1, fs_jpim1
524            ztu3(ji,jj) = ( ztu2(ji+1,jj) - ztu2(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
525         END DO
526         !                     ! Fourth derivative
527         DO ji = fs_2, fs_jpim1
528            ztu4(ji,jj) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj)
529         END DO
530      END DO
531      CALL lbc_lnk( ztu4, 'T', 1. )
532      !
533      !
534      SELECT CASE (kn_umx )
535      !
536      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
537         !       
538         DO jj = 2, jpjm1
539            DO ji = 1, fs_jpim1   ! vector opt.
540               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                              pt(ji+1,jj) + pt(ji,jj)   &
541                  &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )
542            END DO
543         END DO
544         !
545      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
546         !
547         DO jj = 2, jpjm1
548            DO ji = 1, fs_jpim1   ! vector opt.
549               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
550               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                                   pt(ji+1,jj) + pt(ji,jj)   &
551                  &                                               -              zcu   * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
552            END DO
553         END DO
554         ! 
555      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
556         !
557         DO jj = 2, jpjm1
558            DO ji = 1, fs_jpim1   ! vector opt.
559               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
560               zdx2 = e1u(ji,jj) * e1u(ji,jj)
561!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
562               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                         pt  (ji+1,jj) + pt  (ji,jj)        &
563                  &                                               -              zcu   * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   &
564                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj) + ztu2(ji,jj)        &
565                  &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
566            END DO
567         END DO
568         !
569      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
570         !
571         DO jj = 2, jpjm1
572            DO ji = 1, fs_jpim1   ! vector opt.
573               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
574               zdx2 = e1u(ji,jj) * e1u(ji,jj)
575!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
576               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                   pt  (ji+1,jj) + pt  (ji,jj)        &
577                  &                                               -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   &
578                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj) + ztu2(ji,jj)        &
579                  &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
580            END DO
581         END DO
582         !
583      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
584         !
585         DO jj = 2, jpjm1
586            DO ji = 1, fs_jpim1   ! vector opt.
587               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
588               zdx2 = e1u(ji,jj) * e1u(ji,jj)
589!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
590               zdx4 = zdx2 * zdx2
591               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (               (                   pt  (ji+1,jj) + pt  (ji,jj)       &
592                  &                                                     -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) ) )   &
593                  &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj) + ztu2(ji,jj)       &
594                  &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) )   &
595                  &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj) + ztu4(ji,jj)       &
596                  &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj) - ztu4(ji,jj) ) ) )
597            END DO
598         END DO
599         !
600      END SELECT
601      !                                                     !-- High order flux in i-direction  --!
602      DO jj = 1, jpjm1
603         DO ji = 1, fs_jpim1   ! vector opt.
604            pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj)
605         END DO
606      END DO
607      !
608   END SUBROUTINE ultimate_x
609   
610 
611   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
612      !!---------------------------------------------------------------------
613      !!                    ***  ROUTINE ultimate_y  ***
614      !!     
615      !! **  Purpose :   compute 
616      !!
617      !! **  Method  :   ... ???
618      !!                 TIM = transient interpolation Modeling
619      !!
620      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
621      !!----------------------------------------------------------------------
622      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
623      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
624      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pv        ! ice j-velocity component
625      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvc       ! ice j-velocity*A component
626      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
627      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_v      ! tracer at v-point
628      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfv_ho    ! high order flux
629      !
630      INTEGER  ::   ji, jj       ! dummy loop indices
631      REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
632      REAL(wp), DIMENSION(jpi,jpj) ::   ztv1, ztv2, ztv3, ztv4
633      !!----------------------------------------------------------------------
634      !
635      !                                                     !--  Laplacian in j-direction  --!
636      DO jj = 1, jpjm1         ! First derivative (gradient)
637         DO ji = fs_2, fs_jpim1
638            ztv1(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
639         END DO
640      END DO
641      DO jj = 2, jpjm1         ! Second derivative (Laplacian)
642         DO ji = fs_2, fs_jpim1
643            ztv2(ji,jj) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj)
644         END DO
645      END DO
646      CALL lbc_lnk( ztv2, 'T', 1. )
647      !
648      !                                                     !--  BiLaplacian in j-direction  --!
649      DO jj = 1, jpjm1         ! First derivative
650         DO ji = fs_2, fs_jpim1
651            ztv3(ji,jj) = ( ztv2(ji,jj+1) - ztv2(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
652         END DO
653      END DO
654      DO jj = 2, jpjm1         ! Second derivative
655         DO ji = fs_2, fs_jpim1
656            ztv4(ji,jj) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj)
657         END DO
658      END DO
659      CALL lbc_lnk( ztv4, 'T', 1. )
660      !
661      !
662      SELECT CASE (kn_umx )
663      !
664      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
665         DO jj = 1, jpjm1
666            DO ji = fs_2, fs_jpim1
667               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  &
668                  &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
669            END DO
670         END DO
671         !
672      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
673         DO jj = 1, jpjm1
674            DO ji = fs_2, fs_jpim1
675               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
676               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (        ( pt(ji,jj+1) + pt(ji,jj) )  &
677                  &                                     - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) )
678            END DO
679         END DO
680         CALL lbc_lnk( pt_v, 'V',  1. )
681         !
682      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
683         DO jj = 1, jpjm1
684            DO ji = fs_2, fs_jpim1
685               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
686               zdy2 = e2v(ji,jj) * e2v(ji,jj)
687!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
688               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)       &
689                  &                                     -                        zcv   * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   &
690                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1) + ztv2(ji,jj)       &
691                  &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
692            END DO
693         END DO
694         !
695      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
696         DO jj = 1, jpjm1
697            DO ji = fs_2, fs_jpim1
698               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
699               zdy2 = e2v(ji,jj) * e2v(ji,jj)
700!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
701               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                           ( pt  (ji,jj+1) + pt  (ji,jj)       &
702                  &                                               -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   &
703                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                   ztv2(ji,jj+1) + ztv2(ji,jj)       &
704                  &                                               - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
705            END DO
706         END DO
707         !
708      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
709         DO jj = 1, jpjm1
710            DO ji = fs_2, fs_jpim1
711               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
712               zdy2 = e2v(ji,jj) * e2v(ji,jj)
713!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
714               zdy4 = zdy2 * zdy2
715               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)      &
716                  &                                                     -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )  &
717                  &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) *     (                   ztv2(ji,jj+1) + ztv2(ji,jj)      &
718                  &                                                     - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) )  &
719                  &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1) + ztv4(ji,jj)      &
720                  &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1) - ztv4(ji,jj) ) ) )
721            END DO
722         END DO
723         !
724      END SELECT
725      !                                                     !-- High order flux in j-direction  --!
726      DO jj = 1, jpjm1
727         DO ji = 1, fs_jpim1   ! vector opt.
728            pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj)
729         END DO
730      END DO
731      !
732   END SUBROUTINE ultimate_y
733     
734
735   SUBROUTINE nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
736      !!---------------------------------------------------------------------
737      !!                    ***  ROUTINE nonosc  ***
738      !!     
739      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
740      !!       scheme and the before field by a nonoscillatory algorithm
741      !!
742      !! **  Method  :   ... ???
743      !!       warning : ptc and pt_ups must be masked, but the boundaries
744      !!       conditions on the fluxes are not necessary zalezak (1979)
745      !!       drange (1995) multi-dimensional forward-in-time and upstream-
746      !!       in-space based differencing for fluid
747      !!----------------------------------------------------------------------
748      REAL(wp)                     , INTENT(in   ) ::   pdt              ! tracer time-step
749      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt_ups      ! before field & upstream guess of after field
750      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pfv_ups, pfu_ups ! upstream flux
751      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux
752      !
753      INTEGER  ::   ji, jj    ! dummy loop indices
754      REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt   ! local scalars
755      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo         !   -      -
756      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zdiv
757      !!----------------------------------------------------------------------
758      zbig = 1.e+40_wp
759      zsml = 1.e-15_wp
760
761      ! test on divergence
762      DO jj = 2, jpjm1
763         DO ji = fs_2, fs_jpim1   ! vector opt. 
764            zdiv(ji,jj) =  - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) + pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) 
765         END DO
766      END DO
767      CALL lbc_lnk( zdiv, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign)
768
769      ! antidiffusive flux : high order minus low order
770      ! --------------------------------------------------
771      DO jj = 1, jpjm1
772         DO ji = 1, fs_jpim1   ! vector opt.
773            pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_ups(ji,jj)
774            pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_ups(ji,jj)
775         END DO
776      END DO
777
778      ! Search local extrema
779      ! --------------------
780      ! max/min of ptc & pt_ups with large negative/positive value (-/+zbig) outside ice cover
781      DO jj = 1, jpj
782         DO ji = fs_2, fs_jpim1
783            IF( ptc(ji,jj) == 0._wp .AND. pt_ups(ji,jj) == 0._wp .AND. zdiv(ji,jj) == 0._wp ) THEN
784               zbup(ji,jj) = -zbig
785               zbdo(ji,jj) =  zbig
786            ELSE               
787               zbup(ji,jj) = MAX( ptc(ji,jj) , pt_ups(ji,jj) )
788               zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_ups(ji,jj) )
789            ENDIF
790         END DO
791      END DO
792
793      z1_dt = 1._wp / pdt
794      DO jj = 2, jpjm1
795         DO ji = fs_2, fs_jpim1   ! vector opt.
796            !
797            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
798            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) )
799            !
800            zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji  ,jj) ) &
801               & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj  ) )  ! positive/negative part of the flux
802            zneg = MAX( 0., pfu_ho(ji  ,jj) ) - MIN( 0., pfu_ho(ji-1,jj) ) &
803               & + MAX( 0., pfv_ho(ji,jj  ) ) - MIN( 0., pfv_ho(ji,jj-1) )
804            !
805            !                                  ! up & down beta terms
806!!clem            zbetup(ji,jj) = ( zup - pt_ups(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt
807!!clem            zbetdo(ji,jj) = ( pt_ups(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt
808            IF( zpos >= epsi20 ) THEN
809               zbetup(ji,jj) = ( zup - pt_ups(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt
810            ELSE
811               zbetup(ji,jj) = zbig
812            ENDIF
813            !
814            IF( zneg >= epsi20 ) THEN
815               zbetdo(ji,jj) = ( pt_ups(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
816            ELSE
817               zbetdo(ji,jj) = zbig
818            ENDIF
819            !
820         END DO
821      END DO
822      CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
823
824      ! monotonic flux in the y direction
825      ! ---------------------------------
826      DO jj = 1, jpjm1
827         DO ji = 1, fs_jpim1   ! vector opt.
828            zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) )
829            zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) )
830            zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj) )
831            !
832            pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_ups(ji,jj)
833         END DO
834      END DO
835
836      DO jj = 1, jpjm1
837         DO ji = 1, fs_jpim1   ! vector opt.
838            zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) )
839            zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) )
840            zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj) )
841            !
842            pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_ups(ji,jj)
843         END DO
844      END DO
845      !
846   END SUBROUTINE nonosc_2d
847
848   SUBROUTINE limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
849      !!---------------------------------------------------------------------
850      !!                    ***  ROUTINE limiter_x  ***
851      !!     
852      !! **  Purpose :   compute flux limiter
853      !!----------------------------------------------------------------------
854      REAL(wp)                     , INTENT(in   )           ::   pdt          ! tracer time-step
855      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pu           ! ice i-velocity => u*e2
856      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   puc          ! ice i-velocity *A => u*e2*a
857      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pt           ! ice tracer
858      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout)           ::   pfu_ho       ! high order flux
859      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ), OPTIONAL ::   pfu_ups      ! upstream flux
860      !
861      REAL(wp) ::   Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr
862      INTEGER  ::   ji, jj    ! dummy loop indices
863      REAL(wp), DIMENSION (jpi,jpj) ::   zslpx       ! tracer slopes
864      !!----------------------------------------------------------------------
865      !
866      DO jj = 2, jpjm1
867         DO ji = fs_2, fs_jpim1   ! vector opt.
868            zslpx(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * umask(ji,jj,1)
869         END DO
870      END DO
871      CALL lbc_lnk( zslpx, 'U', -1.)   ! lateral boundary cond.
872     
873      DO jj = 2, jpjm1
874         DO ji = fs_2, fs_jpim1   ! vector opt.
875            uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
876
877            Rjm = zslpx(ji-1,jj)
878            Rj  = zslpx(ji  ,jj)
879            Rjp = zslpx(ji+1,jj)
880
881            IF( PRESENT(pfu_ups) ) THEN
882
883               IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm
884               ELSE                        ;   Rr = Rjp
885               ENDIF
886
887               zh3 = pfu_ho(ji,jj) - pfu_ups(ji,jj)     
888               IF( Rj > 0. ) THEN
889                  zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(puc(ji,jj)),  &
890                     &        MIN( 2. * Rr * 0.5 * ABS(puc(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) )
891               ELSE
892                  zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(puc(ji,jj)),  &
893                     &        MIN(-2. * Rr * 0.5 * ABS(puc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) )
894               ENDIF
895               pfu_ho(ji,jj) = pfu_ups(ji,jj) + zlimiter
896               
897            ELSE
898               IF( Rj /= 0. ) THEN
899                  IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
900                  ELSE                        ;   Cr = Rjp / Rj
901                  ENDIF
902               ELSE
903                  Cr = 0.
904                  !IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20
905                  !ELSE                        ;   Cr = Rjp * 1.e20
906                  !ENDIF
907               ENDIF
908               
909               ! -- superbee --
910               zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
911               ! -- van albada 2 --
912               !!zpsi = 2.*Cr / (Cr*Cr+1.)
913
914               ! -- sweby (with beta=1) --
915               !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
916               ! -- van Leer --
917               !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
918               ! -- ospre --
919               !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
920               ! -- koren --
921               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
922               ! -- charm --
923               !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
924               !ELSE                 ;   zpsi = 0.
925               !ENDIF
926               ! -- van albada 1 --
927               !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
928               ! -- smart --
929               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
930               ! -- umist --
931               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
932
933               ! high order flux corrected by the limiter
934               pfu_ho(ji,jj) = pfu_ho(ji,jj) - ABS( puc(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5
935
936            ENDIF
937         END DO
938      END DO
939      CALL lbc_lnk( pfu_ho, 'U', -1.)   ! lateral boundary cond.
940      !
941   END SUBROUTINE limiter_x
942
943   SUBROUTINE limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
944      !!---------------------------------------------------------------------
945      !!                    ***  ROUTINE limiter_y  ***
946      !!     
947      !! **  Purpose :   compute flux limiter
948      !!----------------------------------------------------------------------
949      REAL(wp)                     , INTENT(in   )           ::   pdt          ! tracer time-step
950      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pv           ! ice i-velocity => u*e2
951      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pvc          ! ice i-velocity *A => u*e2*a
952      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pt           ! ice tracer
953      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout)           ::   pfv_ho       ! high order flux
954      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ), OPTIONAL ::   pfv_ups      ! upstream flux
955      !
956      REAL(wp) ::   Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr
957      INTEGER  ::   ji, jj    ! dummy loop indices
958      REAL(wp), DIMENSION (jpi,jpj) ::   zslpy       ! tracer slopes
959      !!----------------------------------------------------------------------
960      !
961      DO jj = 2, jpjm1
962         DO ji = fs_2, fs_jpim1   ! vector opt.
963            zslpy(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * vmask(ji,jj,1)
964         END DO
965      END DO
966      CALL lbc_lnk( zslpy, 'V', -1.)   ! lateral boundary cond.
967     
968      DO jj = 2, jpjm1
969         DO ji = fs_2, fs_jpim1   ! vector opt.
970            vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
971
972            Rjm = zslpy(ji,jj-1)
973            Rj  = zslpy(ji,jj  )
974            Rjp = zslpy(ji,jj+1)
975
976            IF( PRESENT(pfv_ups) ) THEN
977
978               IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm
979               ELSE                        ;   Rr = Rjp
980               ENDIF
981
982               zh3 = pfv_ho(ji,jj) - pfv_ups(ji,jj)     
983               IF( Rj > 0. ) THEN
984                  zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pvc(ji,jj)),  &
985                     &        MIN( 2. * Rr * 0.5 * ABS(pvc(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) )
986               ELSE
987                  zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pvc(ji,jj)),  &
988                     &        MIN(-2. * Rr * 0.5 * ABS(pvc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) )
989               ENDIF
990               pfv_ho(ji,jj) = pfv_ups(ji,jj) + zlimiter
991
992            ELSE
993
994               IF( Rj /= 0. ) THEN
995                  IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
996                  ELSE                        ;   Cr = Rjp / Rj
997                  ENDIF
998               ELSE
999                  Cr = 0.
1000                  !IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20
1001                  !ELSE                        ;   Cr = Rjp * 1.e20
1002                  !ENDIF
1003               ENDIF
1004
1005               ! -- superbee --
1006               zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1007               ! -- van albada 2 --
1008               !!zpsi = 2.*Cr / (Cr*Cr+1.)
1009
1010               ! -- sweby (with beta=1) --
1011               !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1012               ! -- van Leer --
1013               !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1014               ! -- ospre --
1015               !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
1016               ! -- koren --
1017               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
1018               ! -- charm --
1019               !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
1020               !ELSE                 ;   zpsi = 0.
1021               !ENDIF
1022               ! -- van albada 1 --
1023               !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
1024               ! -- smart --
1025               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
1026               ! -- umist --
1027               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1028
1029               ! high order flux corrected by the limiter
1030               pfv_ho(ji,jj) = pfv_ho(ji,jj) - ABS( pvc(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5
1031
1032            ENDIF
1033         END DO
1034      END DO
1035      CALL lbc_lnk( pfv_ho, 'V', -1.)   ! lateral boundary cond.
1036      !
1037   END SUBROUTINE limiter_y
1038
1039#else
1040   !!----------------------------------------------------------------------
1041   !!   Default option           Dummy module         NO SI3 sea-ice model
1042   !!----------------------------------------------------------------------
1043#endif
1044
1045   !!======================================================================
1046END MODULE icedyn_adv_umx
Note: See TracBrowser for help on using the repository browser.