Changeset 10446


Ignore:
Timestamp:
2018-12-21T16:07:19+01:00 (21 months ago)
Author:
clem
Message:

clean the routine icedyn_adv_umx

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icedyn_adv_umx.F90

    r10439 r10446  
    3535   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120 
    3636    
     37   ! limiter: 1=nonosc, 2=superbee, 3=h3(rachid) 
     38   INTEGER ::   kn_limiter = 1 
     39 
    3740   ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme 
    38    LOGICAL :: ll_neg = .TRUE. 
     41   LOGICAL ::   ll_neg = .TRUE. 
    3942    
    4043   ! alternate directions for upstream 
    41    LOGICAL :: ll_upsxy = .TRUE. 
     44   LOGICAL ::   ll_upsxy = .TRUE. 
    4245 
    4346   ! alternate directions for high order 
    44    LOGICAL :: ll_hoxy = .TRUE. 
     47   LOGICAL ::   ll_hoxy = .TRUE. 
    4548    
    4649   ! prelimiter: use it to avoid overshoot in H 
    47    LOGICAL :: ll_prelimiter_zalesak = .TRUE.  ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D 
     50   LOGICAL ::   ll_prelimiter_zalesak = .TRUE.  ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D 
    4851 
    4952 
     
    213216      !!  
    214217      !! **  Purpose :   Compute the now trend due to total advection of  
    215       !!       tracers and add it to the general trend of tracer equations 
     218      !!                 tracers and add it to the general trend of tracer equations 
    216219      !! 
    217       !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with 
    218       !!       corrected flux (monotonic correction) 
    219       !!       note: - this advection scheme needs a leap-frog time scheme 
     220      !! **  Method  :   - calculate upstream fluxes and upstream solution for tracer H 
     221      !!                 - calculate tracer H at u and v points (Ultimate) 
     222      !!                 - calculate the high order fluxes using alterning directions (Macho?) 
     223      !!                 - apply a limiter on the fluxes (nonosc) 
     224      !!                 - convert this tracer flux to a tracer content flux (uH -> uV) 
     225      !!                 - calculate the high order solution for tracer content V 
    220226      !! 
    221       !! ** Action : - pt  the after advective tracer 
     227      !! ** Action : solve 2 equations => a) da/dt = -div(ua) 
     228      !!                                  b) dV/dt = -div(uV) using dH/dt = -u.grad(H) 
     229      !!             in eq. b), - fluxes uH are evaluated (with UMx) and limited (with nonosc). This step is necessary to get a good H. 
     230      !!                        - then we convert this flux to a "volume" flux this way => uH*ua/u 
     231      !!                             where ua is the flux from eq. a) 
     232      !!                        - at last we estimate dV/dt = -div(uH*ua/u) 
     233      !! 
     234      !! ** Note : - this method can lead to small negative V (since we only limit H) => corrected in icedyn_adv.F90 conserving mass etc. 
     235      !!           - negative tracers at u-v points can also occur from the Ultimate scheme (usually at the ice edge) and the solution for now 
     236      !!             is to apply an upstream scheme when it occurs. A better solution would be to degrade the order of 
     237      !!             the scheme automatically by applying a mask of the ice cover inside Ultimate (not done). 
     238      !!           - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good. 
     239      !!             Large values of H can appear for very small ice concentration, and when it does it messes the things up since we 
     240      !!             work on H (and not V). It probably comes from the prelimiter of zalesak which is coded for 1D and not 2D. 
     241      !!             Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 
     242      !!             concentration is small). 
     243      !! To-do: expand the prelimiter from zalesak to make it work in 2D 
    222244      !!---------------------------------------------------------------------- 
    223245      REAL(wp)                        , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0) 
     
    235257      INTEGER  ::   ji, jj, jl       ! dummy loop indices   
    236258      REAL(wp) ::   ztra             ! local scalar 
    237       INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid) 
    238       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt 
     259      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ho , zfv_ho , zpt 
    239260      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ups, zfv_ups, zt_ups 
    240261      !!---------------------------------------------------------------------- 
    241262      ! 
    242       !  upstream (_ups) advection with initial mass fluxes 
    243       ! --------------------------------------------------- 
    244       ! 
    245       IF( .NOT. ll_upsxy ) THEN         !** no alternate directions **! 
    246          DO jl = 1, jpl 
    247             DO jj = 1, jpjm1 
    248                DO ji = 1, fs_jpim1 
    249                   zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
    250                   zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
    251                END DO 
    252             END DO 
    253          END DO 
    254  
    255       ELSE                              !** alternate directions **! 
    256          ! 
    257          IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
    258             ! 
    259             DO jl = 1, jpl              !-- flux in x-direction 
    260                DO jj = 1, jpjm1 
    261                   DO ji = 1, fs_jpim1 
    262                      zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
    263                   END DO 
    264                END DO 
    265             END DO 
    266             ! 
    267             DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    268                DO jj = 2, jpjm1 
    269                   DO ji = fs_2, fs_jpim1 
    270                      zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
    271                         &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
    272                         &            ) * tmask(ji,jj,1) 
    273                   END DO 
    274                END DO 
    275             END DO 
    276             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    277             ! 
    278             DO jl = 1, jpl              !-- flux in y-direction 
    279                DO jj = 1, jpjm1 
    280                   DO ji = 1, fs_jpim1 
    281                      zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 
    282                   END DO 
    283                END DO 
    284             END DO 
    285             ! 
    286          ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
    287             ! 
    288             DO jl = 1, jpl              !-- flux in y-direction 
    289                DO jj = 1, jpjm1 
    290                   DO ji = 1, fs_jpim1 
    291                      zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
    292                   END DO 
    293                END DO 
    294             END DO 
    295             ! 
    296             DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    297                DO jj = 2, jpjm1 
    298                   DO ji = fs_2, fs_jpim1 
    299                      zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 
    300                         &            + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
    301                         &            * tmask(ji,jj,1) 
    302                   END DO 
    303                END DO 
    304             END DO 
    305             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    306             ! 
    307             DO jl = 1, jpl              !-- flux in x-direction 
    308                DO jj = 1, jpjm1 
    309                   DO ji = 1, fs_jpim1 
    310                      zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 
    311                   END DO 
    312                END DO 
    313             END DO 
    314             ! 
    315          ENDIF 
    316           
    317       ENDIF 
    318       ! 
    319       DO jl = 1, jpl                    !-- after tracer with upstream scheme 
    320          DO jj = 2, jpjm1 
    321             DO ji = fs_2, fs_jpim1 
    322                ztra          = - (   zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj  ,jl)   & 
    323                   &                + zfv_ups(ji,jj,jl) - zfv_ups(ji  ,jj-1,jl) ) * r1_e1e2t(ji,jj) 
    324                zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra + ( pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) )   & 
    325                   &                                            +   pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) ) & 
    326                   &                                              * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
    327             END DO 
    328          END DO 
    329       END DO 
    330       CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. ) 
    331  
     263      ! Upstream (_ups) fluxes  
     264      ! ----------------------- 
     265      CALL upstream( pamsk, jt, kt, pdt, pt, pu, pv, zt_ups, zfu_ups, zfv_ups ) 
     266       
    332267      ! High order (_ho) fluxes  
    333268      ! ----------------------- 
     
    336271      CASE ( 20 )                          !== centered second order ==! 
    337272         ! 
    338          CALL cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  & 
    339             &       zt_ups, zfu_ups, zfv_ups ) 
     273         CALL cen2( pamsk, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho ) 
    340274         ! 
    341275      CASE ( 1:5 )                         !== 1st to 5th order ULTIMATE-MACHO scheme ==! 
    342276         ! 
    343          CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho,  & 
    344             &        zt_ups, zfu_ups, zfv_ups ) 
     277         CALL macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho ) 
    345278         ! 
    346279      END SELECT 
     
    372305         END DO 
    373306      ENDIF 
    374       ! 
    375       ! in case of advection of A: output u*a(ho) 
    376       ! ----------------------------------------- 
     307      !                                   --ho 
     308      ! in case of advection of A: output u*a 
     309      ! ------------------------------------- 
    377310      IF( PRESENT( pua_ho ) ) THEN 
    378311         DO jl = 1, jpl 
     
    391324         DO jj = 2, jpjm1 
    392325            DO ji = fs_2, fs_jpim1  
    393                ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) * pdt  
    394                 
    395                ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra ) * tmask(ji,jj,1)                
     326               ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )  
     327               ! 
     328               ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)                
    396329            END DO 
    397330         END DO 
     
    402335 
    403336 
    404    SUBROUTINE cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, & 
    405       &             pt_ups, pfu_ups, pfv_ups ) 
     337   SUBROUTINE upstream( pamsk, jt, kt, pdt, pt, pu, pv, pt_ups, pfu_ups, pfv_ups ) 
    406338      !!--------------------------------------------------------------------- 
    407       !!                    ***  ROUTINE macho  *** 
     339      !!                    ***  ROUTINE upstream  *** 
    408340      !!      
    409       !! **  Purpose :   compute   
    410       !! 
    411       !! **  Method  :   ... ??? 
    412       !!                 TIM = transient interpolation Modeling  
    413       !! 
    414       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
     341      !! **  Purpose :   compute the upstream fluxes and upstream guess of tracer 
    415342      !!---------------------------------------------------------------------- 
    416343      REAL(wp)                        , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
    417       INTEGER                         , INTENT(in   ) ::   kn_limiter       ! limiter 
     344      INTEGER                         , INTENT(in   ) ::   jt               ! number of sub-iteration 
     345      INTEGER                         , INTENT(in   ) ::   kt               ! number of iteration 
     346      REAL(wp)                        , INTENT(in   ) ::   pdt              ! tracer time-step 
     347      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt               ! tracer fields 
     348      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components 
     349      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_ups           ! upstream guess of tracer  
     350      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ups, pfv_ups ! upstream fluxes  
     351      ! 
     352      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     353      REAL(wp) ::   ztra          ! local scalar 
     354      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zpt 
     355      !!---------------------------------------------------------------------- 
     356 
     357      IF( .NOT. ll_upsxy ) THEN         !** no alternate directions **! 
     358         ! 
     359         DO jl = 1, jpl 
     360            DO jj = 1, jpjm1 
     361               DO ji = 1, fs_jpim1 
     362                  pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
     363                  pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
     364               END DO 
     365            END DO 
     366         END DO 
     367         ! 
     368      ELSE                              !** alternate directions **! 
     369         ! 
     370         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     371            ! 
     372            DO jl = 1, jpl              !-- flux in x-direction 
     373               DO jj = 1, jpjm1 
     374                  DO ji = 1, fs_jpim1 
     375                     pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
     376                  END DO 
     377               END DO 
     378            END DO 
     379            ! 
     380            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
     381               DO jj = 2, jpjm1 
     382                  DO ji = fs_2, fs_jpim1 
     383                     ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) )              & 
     384                        &   + ( pu     (ji,jj   ) - pu     (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     385                     ! 
     386                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     387                  END DO 
     388               END DO 
     389            END DO 
     390            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
     391            ! 
     392            DO jl = 1, jpl              !-- flux in y-direction 
     393               DO jj = 1, jpjm1 
     394                  DO ji = 1, fs_jpim1 
     395                     pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 
     396                  END DO 
     397               END DO 
     398            END DO 
     399            ! 
     400         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     401            ! 
     402            DO jl = 1, jpl              !-- flux in y-direction 
     403               DO jj = 1, jpjm1 
     404                  DO ji = 1, fs_jpim1 
     405                     pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
     406                  END DO 
     407               END DO 
     408            END DO 
     409            ! 
     410            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
     411               DO jj = 2, jpjm1 
     412                  DO ji = fs_2, fs_jpim1 
     413                     ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) )  & 
     414                        &   + ( pv     (ji,jj   ) - pv     (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     415                     ! 
     416                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     417                  END DO 
     418               END DO 
     419            END DO 
     420            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
     421            ! 
     422            DO jl = 1, jpl              !-- flux in x-direction 
     423               DO jj = 1, jpjm1 
     424                  DO ji = 1, fs_jpim1 
     425                     pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 
     426                  END DO 
     427               END DO 
     428            END DO 
     429            ! 
     430         ENDIF 
     431          
     432      ENDIF 
     433      ! 
     434      DO jl = 1, jpl                    !-- after tracer with upstream scheme 
     435         DO jj = 2, jpjm1 
     436            DO ji = fs_2, fs_jpim1 
     437               ztra = - (   pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj  ,jl)   & 
     438                  &       + pfv_ups(ji,jj,jl) - pfv_ups(ji  ,jj-1,jl) ) & 
     439                  &   + (   pu     (ji,jj   ) - pu     (ji-1,jj     )   & 
     440                  &       + pv     (ji,jj   ) - pv     (ji  ,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     441               ! 
     442               pt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     443            END DO 
     444         END DO 
     445      END DO 
     446      CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1. ) 
     447 
     448   END SUBROUTINE upstream 
     449 
     450    
     451   SUBROUTINE cen2( pamsk, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     452      !!--------------------------------------------------------------------- 
     453      !!                    ***  ROUTINE cen2  *** 
     454      !!      
     455      !! **  Purpose :   compute the high order fluxes using a centered 
     456      !!                 second order scheme  
     457      !!---------------------------------------------------------------------- 
     458      REAL(wp)                        , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
    418459      INTEGER                         , INTENT(in   ) ::   jt               ! number of sub-iteration 
    419460      INTEGER                         , INTENT(in   ) ::   kt               ! number of iteration 
     
    423464      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components 
    424465      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   ptc              ! tracer content at before time step  
     466      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt_ups           ! upstream guess of tracer  
     467      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes  
    425468      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes  
    426       REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content  
    427       REAL(wp), DIMENSION(:,:,:)      , INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes  
    428469      ! 
    429470      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    430       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zzt 
     471      REAL(wp) ::   ztra          ! local scalar 
     472      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zpt 
    431473      !!---------------------------------------------------------------------- 
    432474      ! 
     
    442484         END DO 
    443485         IF    ( kn_limiter == 1 ) THEN 
    444             CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    445          ELSEIF( kn_limiter == 2 ) THEN 
    446             CALL limiter_x( pdt, pu, pt, pfu_ho ) 
    447             CALL limiter_y( pdt, pv, pt, pfv_ho ) 
    448          ELSEIF( kn_limiter == 3 ) THEN 
    449             CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 
    450             CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 
     486            CALL nonosc( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     487         ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN 
     488            CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     489            CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    451490         ENDIF 
    452491         ! 
     
    462501               END DO 
    463502            END DO 
    464             IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, pt, pfu_ho ) 
    465             IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 
     503            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    466504 
    467505            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    468506               DO jj = 2, jpjm1 
    469507                  DO ji = fs_2, fs_jpim1 
    470                      zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)        & 
    471                         &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
    472                         &         * tmask(ji,jj,1) 
     508                     ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) )              & 
     509                        &   + ( pu    (ji,jj   ) - pu    (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     510                     ! 
     511                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
    473512                  END DO 
    474513               END DO 
    475514            END DO 
    476             CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     515            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    477516 
    478517            DO jl = 1, jpl              !-- flux in y-direction 
    479518               DO jj = 1, jpjm1 
    480519                  DO ji = 1, fs_jpim1 
    481                      pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) ) 
     520                     pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 
    482521                  END DO 
    483522               END DO 
    484523            END DO 
    485             IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pt, pfv_ho ) 
    486             IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 
     524            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    487525 
    488526         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     
    495533               END DO 
    496534            END DO 
    497             IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pt, pfv_ho ) 
    498             IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 
     535            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    499536            ! 
    500537            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    501538               DO jj = 2, jpjm1 
    502539                  DO ji = fs_2, fs_jpim1 
    503                      zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 
    504                         &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
    505                         &         * tmask(ji,jj,1) 
     540                     ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) )  & 
     541                        &   + ( pv    (ji,jj   ) - pv    (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     542                     ! 
     543                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
    506544                  END DO 
    507545               END DO 
    508546            END DO 
    509             CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     547            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    510548            ! 
    511549            DO jl = 1, jpl              !-- flux in x-direction 
    512550               DO jj = 1, jpjm1 
    513551                  DO ji = 1, fs_jpim1 
    514                      pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) ) 
     552                     pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 
    515553                  END DO 
    516554               END DO 
    517555            END DO 
    518             IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, pt, pfu_ho ) 
    519             IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 
     556            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    520557 
    521558         ENDIF 
    522          IF( kn_limiter == 1 )   CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     559         IF( kn_limiter == 1 )   CALL nonosc( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    523560          
    524561      ENDIF 
     
    527564 
    528565    
    529    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, & 
    530       &              pt_ups, pfu_ups, pfv_ups ) 
     566   SUBROUTINE macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    531567      !!--------------------------------------------------------------------- 
    532568      !!                    ***  ROUTINE macho  *** 
    533569      !!      
    534       !! **  Purpose :   compute  
     570      !! **  Purpose :   compute the high order fluxes using Ultimate-Macho scheme  
    535571      !! 
    536       !! **  Method  :   ... ??? 
    537       !!                 TIM = transient interpolation Modeling  
     572      !! **  Method  :   ... 
    538573      !! 
    539574      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    540575      !!---------------------------------------------------------------------- 
    541576      REAL(wp)                        , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
    542       INTEGER                         , INTENT(in   ) ::   kn_limiter       ! limiter 
    543577      INTEGER                         , INTENT(in   ) ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2) 
    544578      INTEGER                         , INTENT(in   ) ::   jt               ! number of sub-iteration 
     
    550584      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pubox, pvbox     ! upstream velocity 
    551585      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   ptc              ! tracer content at before time step  
    552       REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_u, pt_v       ! tracer at u- and v-points  
     586      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt_ups           ! upstream guess of tracer  
     587      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes  
    553588      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes  
    554       REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content  
    555       REAL(wp), DIMENSION(:,:,:)      , INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes  
    556589      ! 
    557590      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    558       REAL(wp) ::   ztra 
    559       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zzt, zzfu_ho, zzfv_ho 
     591      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zt_u, zt_v, zpt 
    560592      !!---------------------------------------------------------------------- 
    561593      ! 
     
    563595         ! 
    564596         !                                                        !--  ultimate interpolation of pt at u-point  --! 
    565          CALL ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 
     597         CALL ultimate_x( kn_umx, pdt, pt, pu, zt_u, pfu_ho ) 
    566598         !                                                        !--  limiter in x --! 
    567          IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, pt, pfu_ho ) 
    568          IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 
    569          !                                                        !--  advective form update in zzt  --! 
     599         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     600         !                                                        !--  advective form update in zpt  --! 
    570601         DO jl = 1, jpl 
    571602            DO jj = 2, jpjm1 
    572603               DO ji = fs_2, fs_jpim1 
    573                   zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj   ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj)  & 
    574                      &                         - pt   (ji,jj,jl) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk 
    575                   zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 
    576                END DO 
    577             END DO 
    578          END DO 
    579          CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     604                  zpt(ji,jj,jl) = ( pt(ji,jj,jl) - (  pubox(ji,jj   ) * ( zt_u(ji,jj,jl) - zt_u(ji-1,jj,jl) ) * r1_e1t  (ji,jj) & 
     605                     &                              + pt   (ji,jj,jl) * ( pu  (ji,jj   ) - pu  (ji-1,jj   ) ) * r1_e1e2t(ji,jj) & 
     606                     &                                                                                        * pamsk           & 
     607                     &                             ) * pdt ) * tmask(ji,jj,1)  
     608               END DO 
     609            END DO 
     610         END DO 
     611         CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    580612         ! 
    581613         !                                                        !--  ultimate interpolation of pt at v-point  --! 
    582614         IF( ll_hoxy ) THEN 
    583             CALL ultimate_y( kn_umx, pdt, zzt, pv, pt_v, pfv_ho ) 
     615            CALL ultimate_y( kn_umx, pdt, zpt, pv, zt_v, pfv_ho ) 
    584616         ELSE 
    585             CALL ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 
     617            CALL ultimate_y( kn_umx, pdt, pt , pv, zt_v, pfv_ho ) 
    586618         ENDIF 
    587619         !                                                        !--  limiter in y --! 
    588          IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pt, pfv_ho ) 
    589          IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 
     620         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    590621         !          
    591622         ! 
     
    593624         ! 
    594625         !                                                        !--  ultimate interpolation of pt at v-point  --! 
    595          CALL ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 
     626         CALL ultimate_y( kn_umx, pdt, pt, pv, zt_v, pfv_ho ) 
    596627         !                                                        !--  limiter in y --! 
    597          IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pt, pfv_ho ) 
    598          IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 
    599          !                                                        !--  advective form update in zzt  --! 
     628         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     629         !                                                        !--  advective form update in zpt  --! 
    600630         DO jl = 1, jpl 
    601631            DO jj = 2, jpjm1 
    602632               DO ji = fs_2, fs_jpim1 
    603                   zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj   ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj)  & 
    604                      &                         - pt   (ji,jj,jl) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk 
    605                   zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 
    606                END DO 
    607             END DO 
    608          END DO 
    609          CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     633                  zpt(ji,jj,jl) = ( pt(ji,jj,jl) - (  pvbox(ji,jj   ) * ( zt_v(ji,jj,jl) - zt_v(ji,jj-1,jl) ) * r1_e2t  (ji,jj) & 
     634                     &                              + pt   (ji,jj,jl) * ( pv  (ji,jj   ) - pv  (ji,jj-1   ) ) * r1_e1e2t(ji,jj) & 
     635                     &                                                                                        * pamsk           & 
     636                     &                             ) * pdt ) * tmask(ji,jj,1)  
     637               END DO 
     638            END DO 
     639         END DO 
     640         CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    610641         ! 
    611642         !                                                        !--  ultimate interpolation of pt at u-point  --! 
    612643         IF( ll_hoxy ) THEN 
    613             CALL ultimate_x( kn_umx, pdt, zzt, pu, pt_u, pfu_ho ) 
     644            CALL ultimate_x( kn_umx, pdt, zpt, pu, zt_u, pfu_ho ) 
    614645         ELSE 
    615             CALL ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 
     646            CALL ultimate_x( kn_umx, pdt, pt , pu, zt_u, pfu_ho ) 
    616647         ENDIF 
    617648         !                                                        !--  limiter in x --! 
    618          IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, pt, pfu_ho ) 
    619          IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 
    620          ! 
     649         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    621650         ! 
    622651      ENDIF 
    623652 
    624       IF( kn_limiter == 1 ) THEN 
    625          CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    626       ENDIF 
     653      IF( kn_limiter == 1 )   CALL nonosc( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    627654      ! 
    628655   END SUBROUTINE macho 
     
    633660      !!                    ***  ROUTINE ultimate_x  *** 
    634661      !!      
    635       !! **  Purpose :   compute  
     662      !! **  Purpose :   compute tracer at u-points  
    636663      !! 
    637       !! **  Method  :   ... ??? 
    638       !!                 TIM = transient interpolation Modeling  
     664      !! **  Method  :   ... 
    639665      !! 
    640666      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
     
    714740                  zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    715741!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    716                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        & 
    717                      &                                               -              zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )  & 
    718                      &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        & 
    719                      &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )  ) 
     742                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     & 
     743                     &                                               -              zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) & 
     744                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     & 
     745                     &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 
    720746               END DO 
    721747            END DO 
     
    730756                  zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    731757!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    732                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        & 
    733                      &                                               -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )  & 
    734                      &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        & 
    735                      &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )  ) 
     758                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     & 
     759                     &                                               -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) & 
     760                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     & 
     761                     &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 
    736762               END DO 
    737763            END DO 
     
    747773!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    748774                  zdx4 = zdx2 * zdx2 
    749                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (            (                   pt  (ji+1,jj,jl) + pt  (ji,jj,jl)       & 
    750                      &                                                     -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) )   & 
    751                      &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)       & 
    752                      &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) )   & 
    753                      &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)       & 
     775                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (            (                   pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     & 
     776                     &                                                     -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) & 
     777                     &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     & 
     778                     &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) & 
     779                     &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)     & 
    754780                     &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 
    755781               END DO 
     
    790816      !!                    ***  ROUTINE ultimate_y  *** 
    791817      !!      
    792       !! **  Purpose :   compute  
     818      !! **  Purpose :   compute tracer at v-points  
    793819      !! 
    794       !! **  Method  :   ... ??? 
    795       !!                 TIM = transient interpolation Modeling  
     820      !! **  Method  :   ... 
    796821      !! 
    797822      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
     
    871896                  zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    872897!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    873                   pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)       & 
    874                      &                                     -                        zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) )   & 
    875                      &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)       & 
     898                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     & 
     899                     &                                     -                        zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) & 
     900                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     & 
    876901                     &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 
    877902               END DO 
     
    886911                  zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    887912!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    888                   pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                        ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)       & 
    889                      &                                               -          zcv * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) )   & 
    890                      &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                   ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)       & 
     913                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                        ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     & 
     914                     &                                               -          zcv * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) & 
     915                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                   ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     & 
    891916                     &                                               - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 
    892917               END DO 
     
    902927!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    903928                  zdy4 = zdy2 * zdy2 
    904                   pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)      & 
    905                      &                                                     -          zcv * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) )  & 
    906                      &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) *     (                   ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)      & 
    907                      &                                                     - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) )  & 
    908                      &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl)      & 
     929                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     & 
     930                     &                                                     -          zcv * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) & 
     931                     &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) *     (                   ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     & 
     932                     &                                                     - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) & 
     933                     &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl)     & 
    909934                     &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) ) 
    910935               END DO 
     
    941966      
    942967 
    943    SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     968   SUBROUTINE nonosc( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    944969      !!--------------------------------------------------------------------- 
    945970      !!                    ***  ROUTINE nonosc  *** 
    946971      !!      
    947        !! **  Purpose :   compute monotonic tracer fluxes from the upstream  
     972      !! **  Purpose :   compute monotonic tracer fluxes from the upstream  
    948973      !!       scheme and the before field by a nonoscillatory algorithm  
    949974      !! 
    950       !! **  Method  :   ... ??? 
    951       !!       warning : pt and pt_ups must be masked, but the boundaries 
    952       !!       conditions on the fluxes are not necessary zalezak (1979) 
    953       !!       drange (1995) multi-dimensional forward-in-time and upstream- 
    954       !!       in-space based differencing for fluid 
     975      !! **  Method  :   ... 
    955976      !!---------------------------------------------------------------------- 
    956977      REAL(wp)                   , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
    957978      REAL(wp)                   , INTENT(in   ) ::   pdt              ! tracer time-step 
    958979      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pu               ! ice i-velocity => u*e2 
    959       REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   puc              ! ice i-velocity *A => u*e2*a 
    960980      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pv               ! ice j-velocity => v*e1 
    961       REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pvc              ! ice j-velocity *A => v*e1*a 
    962       REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   ptc, pt, pt_ups  ! before field & upstream guess of after field 
    963       REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ups, pfu_ups ! upstream flux 
     981      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pt, pt_ups       ! before field & upstream guess of after field 
     982      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pfv_ups, pfu_ups ! upstream flux 
    964983      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux 
    965984      ! 
     
    10041023         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. ) 
    10051024 
    1006          !! this does not work ?? 
    1007          !!            DO jj = 2, jpjm1 
    1008          !!               DO ji = fs_2, fs_jpim1 
    1009          !!                  IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji+1,jj  ) - pt_ups (ji  ,jj) ) .AND.     & 
    1010          !!               &      SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji  ,jj+1) - pt_ups (ji  ,jj) )           & 
    1011          !!               &    ) THEN 
    1012          !!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_ups(ji+1,jj ) - zti_ups(ji  ,jj) ) .AND.   & 
    1013          !!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_ups(ji,jj+1 ) - ztj_ups(ji  ,jj) )         & 
    1014          !!               &       ) THEN 
    1015          !!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0. 
    1016          !!                     ENDIF 
    1017          !!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji  ,jj) - pt_ups (ji-1,jj  ) ) .AND.  & 
    1018          !!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji  ,jj) - pt_ups (ji  ,jj-1) )        & 
    1019          !!               &       )  THEN 
    1020          !!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0. 
    1021          !!                     ENDIF 
    1022          !!                  ENDIF 
    1023          !!                END DO 
    1024          !!            END DO             
    1025  
    10261025         DO jl = 1, jpl 
    10271026            DO jj = 2, jpjm1 
     
    10301029                     & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0. ) THEN 
    10311030                     ! 
    1032                      IF(  pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj,jl) - zti_ups(ji,jj,jl) ) <= 0 .AND.  & 
    1033                         & pfv_ho(ji,jj,jl) * ( ztj_ups(ji,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0) THEN 
     1031                     IF(  pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj,jl) - zti_ups(ji,jj,jl) ) <= 0. .AND.  & 
     1032                        & pfv_ho(ji,jj,jl) * ( ztj_ups(ji,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0. ) THEN 
    10341033                        pfu_ho(ji,jj,jl)=0. 
    10351034                        pfv_ho(ji,jj,jl)=0. 
    10361035                     ENDIF 
    10371036                     ! 
    1038                      IF(  pfu_ho(ji,jj,jl) * ( pt_ups(ji  ,jj,jl) - pt_ups(ji-1,jj,jl) ) <= 0 .AND.  & 
    1039                         & pfv_ho(ji,jj,jl) * ( pt_ups(ji  ,jj,jl) - pt_ups(ji,jj-1,jl) ) <= 0) THEN 
     1037                     IF(  pfu_ho(ji,jj,jl) * ( pt_ups(ji  ,jj,jl) - pt_ups(ji-1,jj,jl) ) <= 0. .AND.  & 
     1038                        & pfv_ho(ji,jj,jl) * ( pt_ups(ji  ,jj,jl) - pt_ups(ji,jj-1,jl) ) <= 0. ) THEN 
    10401039                        pfu_ho(ji,jj,jl)=0. 
    10411040                        pfv_ho(ji,jj,jl)=0. 
     
    10491048 
    10501049      ENDIF 
    1051  
    10521050 
    10531051      ! Search local extrema 
     
    10861084                  & + MAX( 0., pfv_ho(ji,jj  ,jl) ) - MIN( 0., pfv_ho(ji,jj-1,jl) ) 
    10871085               ! 
    1088                zpos = zpos - ( pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) & 
     1086               zpos = zpos - (pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1)) & 
    10891087                  &          ) * ( 1. - pamsk ) 
    1090                zneg = zneg + ( pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) & 
     1088               zneg = zneg + (pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1)) & 
    10911089                  &          ) * ( 1. - pamsk ) 
    10921090               ! 
     
    11541152 
    11551153      END DO 
    1156  
    1157       ! 
    1158       ! 
    1159    END SUBROUTINE nonosc_2d 
    1160  
    1161    SUBROUTINE limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 
     1154      ! 
     1155   END SUBROUTINE nonosc 
     1156 
     1157    
     1158   SUBROUTINE limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    11621159      !!--------------------------------------------------------------------- 
    11631160      !!                    ***  ROUTINE limiter_x  *** 
     
    11651162      !! **  Purpose :   compute flux limiter  
    11661163      !!---------------------------------------------------------------------- 
    1167       REAL(wp)                   , INTENT(in   )          ::   pdt          ! tracer time-step 
    1168       REAL(wp), DIMENSION (:,:  ), INTENT(in   )          ::   pu           ! ice i-velocity => u*e2 
    1169       REAL(wp), DIMENSION (:,:,:), INTENT(in   )          ::   pt           ! ice tracer 
    1170       REAL(wp), DIMENSION (:,:,:), INTENT(inout)           ::   pfu_ho       ! high order flux 
    1171       REAL(wp), DIMENSION (:,:,:), INTENT(in   ), OPTIONAL ::   pfu_ups      ! upstream flux 
     1164      REAL(wp)                  , INTENT(in   ) ::   pdt          ! tracer time-step 
     1165      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu           ! ice i-velocity => u*e2 
     1166      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt           ! ice tracer 
     1167      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pfu_ups      ! upstream flux 
     1168      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pfu_ho       ! high order flux 
    11721169      ! 
    11731170      REAL(wp) ::   Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr 
     
    11941191               Rjp = zslpx(ji+1,jj,jl) 
    11951192 
    1196                IF( PRESENT(pfu_ups) ) THEN 
     1193               IF( kn_limiter == 3 ) THEN 
    11971194 
    11981195                  IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
     
    12101207                  pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 
    12111208 
    1212                ELSE 
     1209               ELSEIF( kn_limiter == 2 ) THEN 
    12131210                  IF( Rj /= 0. ) THEN 
    12141211                     IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj 
     
    12231220                  ! -- van albada 2 -- 
    12241221                  !!zpsi = 2.*Cr / (Cr*Cr+1.) 
    1225  
    12261222                  ! -- sweby (with beta=1) -- 
    12271223                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 
     
    12541250   END SUBROUTINE limiter_x 
    12551251 
    1256    SUBROUTINE limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 
     1252    
     1253   SUBROUTINE limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    12571254      !!--------------------------------------------------------------------- 
    12581255      !!                    ***  ROUTINE limiter_y  *** 
     
    12601257      !! **  Purpose :   compute flux limiter  
    12611258      !!---------------------------------------------------------------------- 
    1262       REAL(wp)                   , INTENT(in   )           ::   pdt          ! tracer time-step 
    1263       REAL(wp), DIMENSION (:,:  ), INTENT(in   )           ::   pv           ! ice i-velocity => u*e2 
    1264       REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pt           ! ice tracer 
    1265       REAL(wp), DIMENSION (:,:,:), INTENT(inout)           ::   pfv_ho       ! high order flux 
    1266       REAL(wp), DIMENSION (:,:,:), INTENT(in   ), OPTIONAL ::   pfv_ups      ! upstream flux 
     1259      REAL(wp)                   , INTENT(in   ) ::   pdt          ! tracer time-step 
     1260      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pv           ! ice i-velocity => u*e2 
     1261      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pt           ! ice tracer 
     1262      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pfv_ups      ! upstream flux 
     1263      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ho       ! high order flux 
    12671264      ! 
    12681265      REAL(wp) ::   Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr 
     
    12891286               Rjp = zslpy(ji,jj+1,jl) 
    12901287 
    1291                IF( PRESENT(pfv_ups) ) THEN 
     1288               IF( kn_limiter == 3 ) THEN 
    12921289 
    12931290                  IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
     
    13051302                  pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 
    13061303 
    1307                ELSE 
     1304               ELSEIF( kn_limiter == 2 ) THEN 
    13081305 
    13091306                  IF( Rj /= 0. ) THEN 
     
    13191316                  ! -- van albada 2 -- 
    13201317                  !!zpsi = 2.*Cr / (Cr*Cr+1.) 
    1321  
    13221318                  ! -- sweby (with beta=1) -- 
    13231319                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 
Note: See TracChangeset for help on using the changeset viewer.