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.
Changeset 10267 – NEMO

Changeset 10267


Ignore:
Timestamp:
2018-10-31T19:09:43+01:00 (6 years ago)
Author:
clem
Message:

improve ice thickness after advection (advection of H using u*a from advection of concentration)

Location:
NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn.F90

    r9604 r10267  
    7171      INTEGER, INTENT(in) ::   kt     ! ice time step 
    7272      !! 
    73       INTEGER ::   ji, jj, jl         ! dummy loop indices 
    74       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhmax 
     73      INTEGER  ::   ji, jj, jl        ! dummy loop indices 
     74      REAL(wp) ::   zcoefu, zcoefv 
     75      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhmax 
     76      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i 
    7577      !!-------------------------------------------------------------------- 
    7678      ! 
     
    118120 
    119121      CASE ( np_dynADV )           !==  pure advection ==!   (prescribed velocities) 
     122         ALLOCATE( zdivu_i(jpi,jpj) ) 
     123 
    120124         u_ice(:,:) = rn_uice * umask(:,:,1) 
    121125         v_ice(:,:) = rn_vice * vmask(:,:,1) 
    122          !!CALL RANDOM_NUMBER(u_ice(:,:)) 
    123          !!CALL RANDOM_NUMBER(v_ice(:,:)) 
     126         !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1) 
     127         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 
     128         ! --- monotonicity test from Lipscomb et al 2004 --- ! 
     129         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length 
     130         ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s  
     131         !DO jj = 1, jpj 
     132         !   DO ji = 1, jpi 
     133         !      zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. ) 
     134         !      zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. ) 
     135         !      u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1) 
     136         !      v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1) 
     137         !   END DO 
     138         !END DO 
     139         ! --- 
    124140         CALL ice_dyn_adv   ( kt )                            ! -- advection of ice 
     141 
     142         ! diagnostics: divergence at T points  
     143         DO jj = 2, jpjm1 
     144            DO ji = 2, jpim1 
     145               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     146                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     147            END DO 
     148         END DO 
     149         CALL lbc_lnk( zdivu_i, 'T', 1. ) 
     150         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) ) 
     151 
     152         DEALLOCATE( zdivu_i ) 
    125153 
    126154      END SELECT 
  • NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90

    r9929 r10267  
    1414   !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders 
    1515   !!   macho             : ??? 
    16    !!   nonosc_2d         : compute monotonic tracer fluxes by a non-oscillatory algorithm  
     16   !!   nonosc            : compute monotonic tracer fluxes by a non-oscillatory algorithm  
    1717   !!---------------------------------------------------------------------- 
    1818   USE phycst         ! physical constant 
     
    2020   USE sbc_oce , ONLY : nn_fsbc   ! update frequency of surface boundary condition 
    2121   USE ice            ! sea-ice variables 
     22   USE icevar         ! sea-ice: operations 
    2223   ! 
    2324   USE in_out_manager ! I/O manager 
     
    4344CONTAINS 
    4445 
    45    SUBROUTINE ice_dyn_adv_umx( k_order, kt, pu_ice, pv_ice,  & 
    46       &                    pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     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 ) 
    4748      !!---------------------------------------------------------------------- 
    4849      !!                  ***  ROUTINE ice_dyn_adv_umx  *** 
     
    5455      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    5556      !!---------------------------------------------------------------------- 
    56       INTEGER                     , INTENT(in   ) ::   k_order    ! order of the scheme (1-5 or 20) 
     57      INTEGER                     , INTENT(in   ) ::   kn_umx     ! order of the scheme (1-5=UM or 20=CEN2) 
    5758      INTEGER                     , INTENT(in   ) ::   kt         ! time step 
    5859      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity 
     
    7071      ! 
    7172      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
    72       INTEGER  ::   initad                  ! number of sub-timestep for the advection 
    73       REAL(wp) ::   zcfl , zusnit, zdt      !   -      - 
    74       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zudy, zvdx, zcu_box, zcv_box 
     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  
    7580      !!---------------------------------------------------------------------- 
    7681      ! 
    7782      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 
    78       ! 
    79       ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) ) 
    8083      ! 
    8184      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !         
     
    8487      IF( lk_mpp )   CALL mpp_max( zcfl ) 
    8588 
    86       IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp 
    87       ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp 
     89      IF( zcfl > 0.5 ) THEN   ;   icycle = 2  
     90      ELSE                    ;   icycle = 1  
    8891      ENDIF 
    8992 
    90       zdt = rdt_ice / REAL(initad) 
     93      zdt = rdt_ice / REAL(icycle) 
    9194 
    9295      ! --- transport --- ! 
     
    112115      !== advection ==! 
    113116      !---------------! 
    114       DO jt = 1, initad 
    115          CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) )             ! Open water area  
     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 
    116122         DO jl = 1, jpl 
    117             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) )         ! Ice area 
    118             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) )         ! Ice  volume 
    119             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, psv_i(:,:,jl) )        ! Salt content 
    120             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i(:,:,jl) )        ! Age content 
     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 
    121142            DO jk = 1, nlay_i 
    122                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) )   ! Ice  heat content 
    123             END DO 
    124             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume 
     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 
    125146            DO jk = 1, nlay_s 
    126                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,jk,jl) )   ! Snow heat content 
    127             END DO 
    128             IF ( ln_pnd_H12 ) THEN 
    129                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction 
    130                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume 
     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 
    131163            ENDIF 
    132          END DO 
    133       END DO 
    134       ! 
    135       DEALLOCATE( zudy, zvdx, zcu_box, zcv_box ) 
     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 
    136184      ! 
    137185   END SUBROUTINE ice_dyn_adv_umx 
    138186 
    139187    
    140    SUBROUTINE adv_umx( k_order, kt, pdt, puc, pvc, pubox, pvbox, ptc ) 
     188   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pfu, pfv ) 
    141189      !!---------------------------------------------------------------------- 
    142190      !!                  ***  ROUTINE adv_umx  *** 
     
    151199      !! ** Action : - pt  the after advective tracer 
    152200      !!---------------------------------------------------------------------- 
    153       INTEGER                     , INTENT(in   ) ::   k_order        ! order of the ULTIMATE scheme 
    154       INTEGER                     , INTENT(in   ) ::   kt             ! number of iteration 
    155       REAL(wp)                    , INTENT(in   ) ::   pdt            ! tracer time-step 
    156       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc  , pvc     ! 2 ice velocity components => u*e2 
    157       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox   ! upstream velocity 
    158       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptc            ! tracer content field 
     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 
    159212      ! 
    160213      INTEGER  ::   ji, jj           ! dummy loop indices   
    161214      REAL(wp) ::   ztra             ! local scalar 
    162       REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfu_ho, zt_u, zt_ups 
    163       REAL(wp), DIMENSION(jpi,jpj) ::   zfv_ups, zfv_ho, zt_v, ztrd 
    164       !!---------------------------------------------------------------------- 
    165       ! 
    166       !  upstream advection with initial mass fluxes & intermediate update 
    167       ! -------------------------------------------------------------------- 
    168       DO jj = 1, jpjm1         ! upstream tracer flux in the i and j direction 
    169          DO ji = 1, fs_jpim1   ! vector opt. 
    170             zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * ptc(ji,jj) + MIN( puc(ji,jj), 0._wp ) * ptc(ji+1,jj) 
    171             zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * ptc(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * ptc(ji,jj+1) 
    172          END DO 
    173       END DO 
    174        
    175       DO jj = 2, jpjm1            ! total intermediate advective trends 
    176          DO ji = fs_2, fs_jpim1   ! vector opt. 
    177             ztra = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   & 
    178                &       + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1)   ) * r1_e1e2t(ji,jj) 
    179             ! 
    180             ztrd(ji,jj) =                         ztra                         ! upstream trend [ -div(uh) or -div(uhT) ]   
    181             zt_ups (ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)      ! guess after content field with monotonic scheme 
    182          END DO 
    183       END DO 
    184       CALL lbc_lnk( zt_ups, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign) 
     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. ) 
    185246       
    186247      ! High order (_ho) fluxes  
    187248      ! ----------------------- 
    188       SELECT CASE( k_order ) 
    189       CASE ( 20 )                          ! centered second order 
    190          DO jj = 1, jpjm1 
    191             DO ji = 1, fs_jpim1   ! vector opt. 
    192                zfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( ptc(ji,jj) + ptc(ji+1,jj) ) 
    193                zfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj) + ptc(ji,jj+1) ) 
    194             END DO 
    195          END DO 
    196          ! 
    197       CASE ( 1:5 )                      ! 1st to 5th order ULTIMATE-MACHO scheme 
    198          CALL macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v ) 
    199          ! 
    200          DO jj = 1, jpjm1 
    201             DO ji = 1, fs_jpim1   ! vector opt. 
    202                zfu_ho(ji,jj) = puc(ji,jj) * zt_u(ji,jj) 
    203                zfv_ho(ji,jj) = pvc(ji,jj) * zt_v(ji,jj) 
    204             END DO 
    205          END DO 
     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 ) 
    206260         ! 
    207261      END SELECT 
    208262          
    209       ! antidiffusive flux : high order minus low order 
    210       ! -------------------------------------------------- 
    211       DO jj = 1, jpjm1 
    212          DO ji = 1, fs_jpim1   ! vector opt. 
    213             zfu_ho(ji,jj) = zfu_ho(ji,jj) - zfu_ups(ji,jj) 
    214             zfv_ho(ji,jj) = zfv_ho(ji,jj) - zfv_ups(ji,jj) 
    215          END DO 
    216       END DO 
    217        
    218       ! monotonicity algorithm 
    219       ! ------------------------- 
    220       CALL nonosc_2d( ptc, zfu_ho, zfv_ho, zt_ups, pdt ) 
     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 
    221275       
    222276      ! final trend with corrected fluxes 
    223277      ! ------------------------------------ 
    224278      DO jj = 2, jpjm1 
    225          DO ji = fs_2, fs_jpim1   ! vector opt.   
    226             ztra       = ztrd(ji,jj)  - (  zfu_ho(ji,jj) - zfu_ho(ji-1,jj  )   & 
    227                &                         + zfv_ho(ji,jj) - zfv_ho(ji  ,jj-1) ) * r1_e1e2t(ji,jj)   
     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)   
    228283            ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
    229284         END DO 
     
    233288   END SUBROUTINE adv_umx 
    234289 
    235  
    236    SUBROUTINE macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 
     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 ) 
    237486      !!--------------------------------------------------------------------- 
    238487      !!                    ***  ROUTINE ultimate_x  *** 
     
    245494      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    246495      !!---------------------------------------------------------------------- 
    247       INTEGER                     , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme 
    248       INTEGER                     , INTENT(in   ) ::   kt         ! number of iteration 
    249       REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step 
    250       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc        ! tracer fields 
    251       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc   ! 2 ice velocity components 
    252       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox   ! upstream velocity 
    253       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u, pt_v ! tracer at u- and v-points  
    254       ! 
    255       INTEGER  ::   ji, jj    ! dummy loop indices 
    256       REAL(wp) ::   zc_box    !   -      - 
    257       REAL(wp), DIMENSION(jpi,jpj) :: zzt 
    258       !!---------------------------------------------------------------------- 
    259       ! 
    260       IF( MOD( (kt - 1) / nn_fsbc , 2 ) == 0 ) THEN         !==  odd ice time step:  adv_x then adv_y  ==! 
    261          ! 
    262          !                                                           !--  ultimate interpolation of pt at u-point  --! 
    263          CALL ultimate_x( k_order, pdt, ptc, puc, pt_u ) 
    264          ! 
    265          !                                                           !--  advective form update in zzt  --! 
    266          DO jj = 2, jpjm1 
    267             DO ji = fs_2, fs_jpim1   ! vector opt. 
    268                zzt(ji,jj) = ptc(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  & 
    269                   &                    - ptc  (ji,jj) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj) 
    270                zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 
    271             END DO 
    272          END DO 
    273          CALL lbc_lnk( zzt, 'T', 1. ) 
    274          ! 
    275          !                                                           !--  ultimate interpolation of pt at v-point  --! 
    276          CALL ultimate_y( k_order, pdt, zzt, pvc, pt_v ) 
    277          ! 
    278       ELSE                                                  !==  even ice time step:  adv_y then adv_x  ==! 
    279          ! 
    280          !                                                           !--  ultimate interpolation of pt at v-point  --! 
    281          CALL ultimate_y( k_order, pdt, ptc, pvc, pt_v ) 
    282          ! 
    283          !                                                           !--  advective form update in zzt  --! 
    284          DO jj = 2, jpjm1 
    285             DO ji = fs_2, fs_jpim1 
    286                zzt(ji,jj) = ptc(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  & 
    287                   &                    - ptc  (ji,jj) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj) 
    288                zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 
    289             END DO 
    290          END DO 
    291          CALL lbc_lnk( zzt, 'T', 1. ) 
    292          ! 
    293          !                                                           !--  ultimate interpolation of pt at u-point  --! 
    294          CALL ultimate_x( k_order, pdt, zzt, puc, pt_u ) 
    295          !       
    296       ENDIF       
    297       ! 
    298    END SUBROUTINE macho 
    299  
    300  
    301    SUBROUTINE ultimate_x( k_order, pdt, pt, puc, pt_u ) 
    302       !!--------------------------------------------------------------------- 
    303       !!                    ***  ROUTINE ultimate_x  *** 
    304       !!      
    305       !! **  Purpose :   compute   
    306       !! 
    307       !! **  Method  :   ... ??? 
    308       !!                 TIM = transient interpolation Modeling  
    309       !! 
    310       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    311       !!---------------------------------------------------------------------- 
    312       INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index 
     496      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
    313497      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step 
    314       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc       ! ice i-velocity component 
     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 
    315500      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields 
    316501      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u      ! tracer at u-point  
    317       ! 
    318       INTEGER  ::   ji, jj       ! dummy loop indices 
     502      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho    ! high order flux  
     503      ! 
     504      INTEGER  ::   ji, jj             ! dummy loop indices 
    319505      REAL(wp) ::   zcu, zdx2, zdx4    !   -      - 
    320       REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu2, ztu3, ztu4 
     506      REAL(wp), DIMENSION(jpi,jpj) ::   ztu1, ztu2, ztu3, ztu4 
    321507      !!---------------------------------------------------------------------- 
    322508      ! 
     
    346532      ! 
    347533      ! 
    348       SELECT CASE (k_order ) 
     534      SELECT CASE (kn_umx ) 
    349535      ! 
    350536      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21) 
     
    352538         DO jj = 2, jpjm1 
    353539            DO ji = 1, fs_jpim1   ! vector opt. 
    354                pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj) + pt(ji,jj)   & 
    355                   &                                    - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
     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) ) ) 
    356542            END DO 
    357543         END DO 
     
    361547         DO jj = 2, jpjm1 
    362548            DO ji = 1, fs_jpim1   ! vector opt. 
    363                zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     549               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    364550               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                                   pt(ji+1,jj) + pt(ji,jj)   & 
    365551                  &                                               -              zcu   * ( pt(ji+1,jj) - pt(ji,jj) ) )  
     
    371557         DO jj = 2, jpjm1 
    372558            DO ji = 1, fs_jpim1   ! vector opt. 
    373                zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     559               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    374560               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    375561!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
     
    385571         DO jj = 2, jpjm1 
    386572            DO ji = 1, fs_jpim1   ! vector opt. 
    387                zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     573               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    388574               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    389575!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
     
    399585         DO jj = 2, jpjm1 
    400586            DO ji = 1, fs_jpim1   ! vector opt. 
    401                zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     587               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    402588               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    403589!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
     
    413599         ! 
    414600      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 
    415607      ! 
    416608   END SUBROUTINE ultimate_x 
    417609    
    418610  
    419    SUBROUTINE ultimate_y( k_order, pdt, pt, pvc, pt_v ) 
     611   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 
    420612      !!--------------------------------------------------------------------- 
    421613      !!                    ***  ROUTINE ultimate_y  *** 
     
    428620      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    429621      !!---------------------------------------------------------------------- 
    430       INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index 
     622      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
    431623      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step 
    432       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvc       ! ice j-velocity component 
     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 
    433626      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields 
    434627      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  
    435629      ! 
    436630      INTEGER  ::   ji, jj       ! dummy loop indices 
    437631      REAL(wp) ::   zcv, zdy2, zdy4    !   -      - 
    438       REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv2, ztv3, ztv4 
     632      REAL(wp), DIMENSION(jpi,jpj) ::   ztv1, ztv2, ztv3, ztv4 
    439633      !!---------------------------------------------------------------------- 
    440634      ! 
     
    466660      ! 
    467661      ! 
    468       SELECT CASE (k_order ) 
     662      SELECT CASE (kn_umx ) 
    469663      ! 
    470664      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21) 
    471665         DO jj = 1, jpjm1 
    472666            DO ji = fs_2, fs_jpim1 
    473                pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1) + pt(ji,jj) )  & 
    474                   &                                     - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 
     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) ) ) 
    475669            END DO 
    476670         END DO 
     
    479673         DO jj = 1, jpjm1 
    480674            DO ji = fs_2, fs_jpim1 
    481                zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     675               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    482676               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (        ( pt(ji,jj+1) + pt(ji,jj) )  & 
    483677                  &                                     - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) ) 
     
    489683         DO jj = 1, jpjm1 
    490684            DO ji = fs_2, fs_jpim1 
    491                zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     685               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    492686               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    493687!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
     
    502696         DO jj = 1, jpjm1 
    503697            DO ji = fs_2, fs_jpim1 
    504                zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     698               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    505699               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    506700!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
     
    515709         DO jj = 1, jpjm1 
    516710            DO ji = fs_2, fs_jpim1 
    517                zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     711               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    518712               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    519713!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
     
    529723         ! 
    530724      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 
    531731      ! 
    532732   END SUBROUTINE ultimate_y 
    533     
    534    
    535    SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, pdt ) 
     733      
     734 
     735   SUBROUTINE nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    536736      !!--------------------------------------------------------------------- 
    537737      !!                    ***  ROUTINE nonosc  *** 
     
    541741      !! 
    542742      !! **  Method  :   ... ??? 
    543       !!       warning : pbef and paft must be masked, but the boundaries 
     743      !!       warning : ptc and pt_ups must be masked, but the boundaries 
    544744      !!       conditions on the fluxes are not necessary zalezak (1979) 
    545745      !!       drange (1995) multi-dimensional forward-in-time and upstream- 
    546746      !!       in-space based differencing for fluid 
    547747      !!---------------------------------------------------------------------- 
    548       REAL(wp)                     , INTENT(in   ) ::   pdt          ! tracer time-step 
    549       REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pbef, paft   ! before & after field 
    550       REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   paa, pbb     ! monotonic fluxes in the 2 directions 
     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 
    551752      ! 
    552753      INTEGER  ::   ji, jj    ! dummy loop indices 
    553       INTEGER  ::   ikm1      ! local integer 
    554       REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zsml, z1_dt   ! local scalars 
    555       REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      - 
    556       REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zmsk, zdiv 
    557       !!---------------------------------------------------------------------- 
    558       ! 
     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      !!---------------------------------------------------------------------- 
    559758      zbig = 1.e+40_wp 
    560759      zsml = 1.e-15_wp 
     
    563762      DO jj = 2, jpjm1 
    564763         DO ji = fs_2, fs_jpim1   ! vector opt.   
    565             zdiv(ji,jj) =  - (  paa(ji,jj) - paa(ji-1,jj  )   & 
    566                &              + pbb(ji,jj) - pbb(ji  ,jj-1) )   
     764            zdiv(ji,jj) =  - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) + pfu_ho(ji,jj) - pfu_ho(ji-1,jj) )   
    567765         END DO 
    568766      END DO 
    569767      CALL lbc_lnk( zdiv, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign) 
    570768 
    571       ! Determine ice masks for before and after tracers  
    572       WHERE( pbef(:,:) == 0._wp .AND. paft(:,:) == 0._wp .AND. zdiv(:,:) == 0._wp )   ;   zmsk(:,:) = 0._wp 
    573       ELSEWHERE                                                                       ;   zmsk(:,:) = 1._wp * tmask(:,:,1) 
    574       END WHERE 
     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 
    575777 
    576778      ! Search local extrema 
    577779      ! -------------------- 
    578       ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 
    579 !      zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ),   & 
    580 !         &             paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) )  ) 
    581 !      zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ),   & 
    582 !         &             paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) )  ) 
    583       zbup(:,:) = MAX( pbef(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ),   & 
    584          &             paft(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) )  ) 
    585       zbdo(:,:) = MIN( pbef(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ),   & 
    586          &             paft(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) )  ) 
     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 
    587792 
    588793      z1_dt = 1._wp / pdt 
     
    590795         DO ji = fs_2, fs_jpim1   ! vector opt. 
    591796            ! 
    592             zup  = MAX(   zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ),   &        ! search max/min in neighbourhood 
    593                &                       zbup(ji  ,jj-1), zbup(ji  ,jj+1)    ) 
    594             zdo  = MIN(   zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ),   & 
    595                &                       zbdo(ji  ,jj-1), zbdo(ji  ,jj+1)    ) 
    596                ! 
    597             zpos = MAX( 0., paa(ji-1,jj  ) ) - MIN( 0., paa(ji  ,jj  ) )   &        ! positive/negative  part of the flux 
    598                & + MAX( 0., pbb(ji  ,jj-1) ) - MIN( 0., pbb(ji  ,jj  ) ) 
    599             zneg = MAX( 0., paa(ji  ,jj  ) ) - MIN( 0., paa(ji-1,jj  ) )   & 
    600                & + MAX( 0., pbb(ji  ,jj  ) ) - MIN( 0., pbb(ji  ,jj-1) ) 
    601                ! 
    602             zbt = e1e2t(ji,jj) * z1_dt                                   ! up & down beta terms 
    603             zbetup(ji,jj) = ( zup         - paft(ji,jj) ) / ( zpos + zsml ) * zbt 
    604             zbetdo(ji,jj) = ( paft(ji,jj) - zdo         ) / ( zneg + zsml ) * zbt 
     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            ! 
    605820         END DO 
    606821      END DO 
    607822      CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign) 
    608823 
    609       ! monotonic flux in the i & j direction (paa & pbb) 
    610       ! ------------------------------------- 
    611       DO jj = 2, jpjm1 
     824      ! monotonic flux in the y direction 
     825      ! --------------------------------- 
     826      DO jj = 1, jpjm1 
    612827         DO ji = 1, fs_jpim1   ! vector opt. 
    613828            zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) ) 
    614829            zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) ) 
    615             zcu = 0.5  + SIGN( 0.5 , paa(ji,jj) ) 
    616             ! 
    617             paa(ji,jj) = paa(ji,jj) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 
    618          END DO 
    619       END DO 
    620       ! 
     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 
    621836      DO jj = 1, jpjm1 
    622          DO ji = fs_2, fs_jpim1   ! vector opt. 
     837         DO ji = 1, fs_jpim1   ! vector opt. 
    623838            zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) ) 
    624839            zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) ) 
    625             zcv = 0.5  + SIGN( 0.5 , pbb(ji,jj) ) 
    626             ! 
    627             pbb(ji,jj) = pbb(ji,jj) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 
     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) 
    628843         END DO 
    629844      END DO 
    630845      ! 
    631846   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 
    6321038 
    6331039#else 
Note: See TracChangeset for help on using the changeset viewer.