Ignore:
Timestamp:
2018-12-17T12:25:09+01:00 (2 years ago)
Author:
clem
Message:

improve ice advection (toward an acceptable solution)

File:
1 edited

Legend:

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

    r10331 r10399  
    3535   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120 
    3636 
    37    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_a, z1_a_ups, zua_ups, zva_ups 
     37   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_ai, amaxu, amaxv 
    3838    
    39    ! rachid trick (only for nonosc) 
    40    LOGICAL :: ll_rachid = .TRUE. 
    41  
     39   LOGICAL ll_dens 
     40 
     41   ! advect H all the way (and get V=H*A at the end) 
     42   LOGICAL :: ll_thickness = .FALSE. 
     43    
     44   ! look for 9 points around in nonosc limiter   
     45   LOGICAL :: ll_9points = .FALSE.  ! false=better h? 
     46 
     47   ! use HgradU in nonosc limiter   
     48   LOGICAL :: ll_HgradU = .TRUE.   ! no effect? 
     49 
     50   ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme 
     51   LOGICAL :: ll_neg = .TRUE.      ! keep TRUE 
     52    
     53   ! limit the fluxes 
     54   LOGICAL :: ll_zeroup1 = .FALSE. ! false ok if Hbig otherwise needed for 2D sinon on a des valeurs de H trop fortes !! 
     55   LOGICAL :: ll_zeroup2 = .FALSE. ! false ok for 1D, 2D, 3D 
     56   LOGICAL :: ll_zeroup4 = .FALSE. ! false ok for 1D, 2D, 3D 
     57   LOGICAL :: ll_zeroup5 = .FALSE. ! false ok for 1D, 2D 
     58 
     59   ! fluxes that are limited are u*H, then (u*H)*(ua/u) is used for V (only for nonosc) 
     60   LOGICAL :: ll_clem   = .TRUE.   ! simpler than rachid and works 
     61 
     62   ! First advect H as H*=Hdiv(u), then use H* for H(n+1)=H(n)-div(uH*) 
     63   LOGICAL :: ll_gurvan = .FALSE.  ! must be false for 1D case !! 
     64 
     65   ! First guess as div(uH) (-true-) or Hdiv(u)+ugradH (-false-) 
     66   LOGICAL :: ll_1stguess_clem = .FALSE. ! better negative values but less good h 
     67 
     68   ! advect (or not) open water. If not, retrieve it from advection of A 
     69   LOGICAL :: ll_ADVopw = .FALSE.  ! 
     70    
    4271   ! alternate directions for upstream 
    43    ! clem: it gives results for Lipscomb test that are the same as "ll_upsxy=false" 
    44    ! clem: needs to be set to true in 2D when using prelimiter (otherwise "wavy solutions" are created) 
    4572   LOGICAL :: ll_upsxy = .TRUE. 
     73 
     74   ! alternate directions for high order 
     75   LOGICAL :: ll_hoxy = .TRUE. 
    4676    
    47    ! prelimiter 
    48    ! clem: use it to avoid overshoot in H 
    49    LOGICAL :: ll_prelimiter_zalesak = .TRUE.  ! from: Zalesak(1979) eq. 14 => better than Devore especially in 2D 
    50    LOGICAL :: ll_prelimiter_devore  = .FALSE.  ! from: Devore eq. 11 
     77   ! prelimiter: use it to avoid overshoot in H 
     78   LOGICAL :: ll_prelimiter_zalesak = .TRUE.  ! from: Zalesak(1979) eq. 14 => true is better for 1D but false is better in 3D (for h and negative values) => pb in x-y? 
     79   LOGICAL :: ll_prelimiter_devore  = .FALSE.  ! from: Devore eq. 11 (worth than zalesak) 
    5180 
    5281   ! iterate on the limiter (only nonosc) 
    53    ! clem: useless 
    5482   LOGICAL :: ll_limiter_it2 = .FALSE. 
    5583    
     
    94122      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
    95123      REAL(wp) ::   zcfl , zdt 
    96       REAL(wp) ::   zeps = 0.0_wp           ! shift in concentration to avoid division by 0 
    97       !                                     !    must be >= 0.01 and the best seems to be 0.1 
    98       REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx, zcu_box, zcv_box, zua_ho, zva_ho, za_ups 
    99       REAL(wp), DIMENSION(jpi,jpj) ::   zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, zh_ip  
     124      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx, zcu_box, zcv_box, zua_ho, zva_ho 
     125      REAL(wp), DIMENSION(jpi,jpj) ::   zhvar 
     126      REAL(wp), DIMENSION(jpi,jpj) ::   zai_b, zai_a, z1_ai_b 
    100127      !!---------------------------------------------------------------------- 
    101128      ! 
    102129      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 
     130      ! 
    103131      ! 
    104132      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !         
     
    132160      END DO 
    133161 
    134       IF(.NOT. ALLOCATED(z1_a)    )   ALLOCATE(z1_a    (jpi,jpj)) 
    135       IF(.NOT. ALLOCATED(z1_a_ups))   ALLOCATE(z1_a_ups(jpi,jpj)) 
    136       IF(.NOT. ALLOCATED(zua_ups) )   ALLOCATE(zua_ups (jpi,jpj)) 
    137       IF(.NOT. ALLOCATED(zva_ups) )   ALLOCATE(zva_ups (jpi,jpj)) 
     162      IF(.NOT. ALLOCATED(z1_ai))       ALLOCATE(z1_ai(jpi,jpj)) 
     163      IF(.NOT. ALLOCATED(amaxu))       ALLOCATE(amaxu (jpi,jpj)) 
     164      IF(.NOT. ALLOCATED(amaxv))       ALLOCATE(amaxv (jpi,jpj)) 
     165 
    138166      !---------------! 
    139167      !== advection ==! 
     
    141169      DO jt = 1, icycle 
    142170 
    143          zamsk = 1._wp ; zua_ups(:,:) = zudy(:,:) ; zva_ups(:,:) = zvdx(:,:)  
    144          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:), zua_ups, zva_ups, za_ups, zua_ho, zva_ho )   ! Open water area  
    145  
     171         IF( ll_ADVopw ) THEN 
     172            ll_dens=.FALSE. 
     173            zamsk = 1._wp 
     174            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) )   ! Open water area  
     175         ELSE 
     176            zai_b(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
     177         ENDIF 
     178          
    146179         DO jl = 1, jpl 
    147             ! to avoid a problem with the limiter nonosc when A gets close to 0 
    148             pa_i(:,:,jl) = pa_i(:,:,jl) + zeps * tmask(:,:,1) 
    149             ! 
    150             WHERE( pa_i(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_i(:,:,jl) 
    151             ELSEWHERE                        ;   z1_a(:,:) = 0. 
     180            ! 
     181            WHERE( pa_i(:,:,jl) >= epsi20 )   ;   z1_ai_b(:,:) = 1._wp / pa_i(:,:,jl) 
     182            ELSEWHERE                         ;   z1_ai_b(:,:) = 0. 
    152183            END WHERE 
    153184            ! 
    154             zamsk = 1._wp ; zua_ups(:,:) = zudy(:,:) ; zva_ups(:,:) = zvdx(:,:)  
    155             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zua_ups, zva_ups, za_ups, zua_ho, zva_ho ) ! Ice area 
    156  
    157             ! 1/A_ups 
    158 !!            IF( .NOT. ll_rachid ) THEN 
    159                WHERE( za_ups(:,:) > epsi20 )   ;   z1_a_ups(:,:) = 1._wp / za_ups(:,:) 
    160                ELSEWHERE                       ;   z1_a_ups(:,:) = 0. 
    161                END WHERE 
    162 !!            ELSE 
    163 !!               WHERE( pa_i(:,:,jl) > epsi20 )   ;   z1_a_ups(:,:) = 1._wp / pa_i(:,:,jl) 
    164 !!               ELSEWHERE                        ;   z1_a_ups(:,:) = 0. 
    165 !!               END WHERE               
    166 !!            ENDIF 
    167 !! 
    168 !!            IF( ll_rachid )   zua_ups = zua_ho ; zva_ups = zva_ho 
    169             zamsk = 0._wp ; zh_i(:,:) = pv_i (:,:,jl) * z1_a(:,:) 
    170             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zh_i(:,:), pv_i (:,:,jl), zua_ups, zva_ups )              ! Ice volume 
    171             zamsk = 0._wp ; zh_s(:,:) = pv_s (:,:,jl) * z1_a(:,:) 
    172             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zh_s(:,:), pv_s (:,:,jl), zua_ups, zva_ups )              ! Snw volume 
    173             zamsk = 0._wp ; zs_i(:,:) = psv_i(:,:,jl) * z1_a(:,:) 
    174             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zs_i(:,:), psv_i(:,:,jl), zua_ups, zva_ups )              ! Salt content 
    175             zamsk = 0._wp ; zo_i(:,:) = poa_i(:,:,jl) * z1_a(:,:) 
    176             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zo_i(:,:), poa_i(:,:,jl), zua_ups, zva_ups )              ! Age content 
     185            IF( ll_zeroup2 ) THEN 
     186               DO jj = 2, jpjm1 
     187                  DO ji = fs_2, fs_jpim1 
     188                     amaxu(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji,jj-1,jl), pa_i(ji,jj+1,jl), & 
     189                        &                              pa_i(ji+1,jj,jl), pa_i(ji+1,jj-1,jl), pa_i(ji+1,jj+1,jl) ) 
     190                     amaxv(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji-1,jj,jl), pa_i(ji+1,jj,jl), & 
     191                        &                              pa_i(ji,jj+1,jl), pa_i(ji-1,jj+1,jl), pa_i(ji+1,jj+1,jl) ) 
     192                 END DO 
     193               END DO 
     194               CALL lbc_lnk_multi(amaxu, 'T', 1., amaxv, 'T', 1.) 
     195            ENDIF 
     196            ! 
     197            zamsk = 1._wp 
     198            ll_dens=.TRUE. 
     199            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zua_ho, zva_ho ) ! Ice area 
     200            ll_dens=.FALSE. 
     201 
     202            WHERE( pa_i(:,:,jl) >= epsi20 )   ;   z1_ai(:,:) = 1._wp / pa_i(:,:,jl) 
     203            ELSEWHERE                         ;   z1_ai(:,:) = 0. 
     204            END WHERE               
     205 
     206            IF( ll_thickness ) THEN 
     207               zua_ho(:,:) = zudy(:,:) 
     208               zva_ho(:,:) = zvdx(:,:) 
     209            ENDIF 
     210             
     211            zamsk = 0._wp ; zhvar(:,:) = pv_i (:,:,jl) * z1_ai_b(:,:) 
     212            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_i (:,:,jl) )              ! Ice volume 
     213            IF( ll_thickness )   pv_i(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl) 
     214             
     215            zamsk = 0._wp ; zhvar(:,:) = pv_s (:,:,jl) * z1_ai_b(:,:) 
     216            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_s (:,:,jl) )              ! Snw volume 
     217            IF( ll_thickness )   pv_s(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl) 
     218 
     219            zamsk = 0._wp ; zhvar(:,:) = psv_i(:,:,jl) * z1_ai_b(:,:) 
     220            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), psv_i(:,:,jl) )              ! Salt content 
     221 
     222            zamsk = 0._wp ; zhvar(:,:) = poa_i(:,:,jl) * z1_ai_b(:,:) 
     223            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), poa_i(:,:,jl) )              ! Age content 
     224 
    177225            DO jk = 1, nlay_i 
    178                zamsk = 0._wp ; ze_i(:,:) = pe_i(:,:,jk,jl) * z1_a(:,:) 
    179                CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, ze_i(:,:), pe_i(:,:,jk,jl), zua_ups, zva_ups )           ! Ice heat content 
    180             END DO 
     226               zamsk = 0._wp ; zhvar(:,:) = pe_i(:,:,jk,jl) * z1_ai_b(:,:) 
     227               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar(:,:), pe_i(:,:,jk,jl) )           ! Ice heat content 
     228            END DO 
     229 
    181230            DO jk = 1, nlay_s 
    182                zamsk = 0._wp ; ze_s(:,:) = pe_s(:,:,jk,jl) * z1_a(:,:) 
    183                CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, ze_s(:,:), pe_s(:,:,jk,jl), zua_ups, zva_ups )           ! Snw heat content 
    184             END DO 
    185             ! 
    186             IF ( ln_pnd_H12 ) THEN   ! melt ponds (must be the last ones to be advected because of z1_a...) 
    187                ! to avoid a problem with the limiter nonosc when A gets close to 0 
    188                pa_ip(:,:,jl) = pa_ip(:,:,jl) + zeps * tmask(:,:,1) 
     231               zamsk = 0._wp ; zhvar(:,:) = pe_s(:,:,jk,jl) * z1_ai_b(:,:) 
     232               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar(:,:), pe_s(:,:,jk,jl) )           ! Snw heat content 
     233            END DO 
     234            ! 
     235            IF ( ln_pnd_H12 ) THEN   ! melt ponds (must be the last ones to be advected because of z1_ai_b...) 
    189236               ! 
    190                WHERE( pa_ip(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_ip(:,:,jl) 
    191                ELSEWHERE                         ;   z1_a(:,:) = 0. 
     237               WHERE( pa_ip(:,:,jl) >= epsi20 )   ;   z1_ai_b(:,:) = 1._wp / pa_ip(:,:,jl) 
     238               ELSEWHERE                          ;   z1_ai_b(:,:) = 0. 
    192239               END WHERE 
    193240               ! 
    194                zamsk = 1._wp ; zua_ups(:,:) = zudy(:,:) ; zva_ups(:,:) = zvdx(:,:)  
    195                CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zua_ups, zva_ups, za_ups, zua_ho, zva_ho ) ! mp fraction 
    196  
    197                ! 1/A_ups 
    198                WHERE( za_ups(:,:) > epsi20 )   ;   z1_a_ups(:,:) = 1._wp / za_ups(:,:) 
    199                ELSEWHERE                       ;   z1_a_ups(:,:) = 0. 
    200                END WHERE 
    201  
    202 !!            IF( ll_rachid )   zua_ups = zua_ho ; zva_ups = zva_ho 
    203                zamsk = 0._wp ; zh_ip(:,:) = pv_ip(:,:,jl) * z1_a(:,:) 
    204                CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zh_ip(:,:), pv_ip(:,:,jl), zua_ups, zva_ups )              ! mp volume 
    205             ENDIF 
    206             ! 
    207             ! to avoid a problem with the limiter nonosc when A gets close to 0 
     241               zamsk = 1._wp 
     242               ll_dens=.TRUE. 
     243               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zua_ho, zva_ho ) ! mp fractio!n 
     244               ll_dens=.FALSE. 
     245 
     246               WHERE( pa_ip(:,:,jl) >= epsi20 )   ;   z1_ai(:,:) = 1._wp / pa_ip(:,:,jl) 
     247               ELSEWHERE                          ;   z1_ai(:,:) = 0. 
     248               END WHERE               
     249 
     250               zamsk = 0._wp ; zhvar(:,:) = pv_ip(:,:,jl) * z1_ai_b(:,:) 
     251               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_ip(:,:,jl) )              ! mp volume 
     252            ENDIF 
     253            ! 
     254            ! 
     255         END DO 
     256 
     257         IF( .NOT. ll_ADVopw ) THEN 
     258            zai_a(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
    208259            DO jj = 2, jpjm1 
    209260               DO ji = fs_2, fs_jpim1 
    210                   !pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps ) * tmask(ji,jj,1) 
    211                   pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps & 
    212                      &             + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt & 
    213                      &             ) * tmask(ji,jj,1) 
    214                   IF ( ln_pnd_H12 ) THEN   ! melt ponds 
    215                      pa_ip(ji,jj,jl) = ( pa_ip(ji,jj,jl) - zeps & 
    216                         &              + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt & 
    217                         &              ) * tmask(ji,jj,1) 
    218                   ENDIF 
    219               END DO 
    220             END DO 
    221             CALL lbc_lnk_multi( pa_i(:,:,jl), 'T',  1., pa_ip(:,:,jl), 'T',  1. ) 
    222             ! 
    223          END DO 
    224  
     261                  pato_i(ji,jj) = pato_i(ji,jj) - ( zai_a(ji,jj) - zai_b(ji,jj) ) & 
     262                     &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
     263               END DO 
     264            END DO 
     265            CALL lbc_lnk( pato_i(:,:), 'T',  1. ) 
     266         ENDIF 
     267          
    225268      END DO 
    226269      ! 
     
    228271 
    229272    
    230    SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ups, pva_ups, pa_ups, pua_ho, pva_ho ) 
     273   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho ) 
    231274      !!---------------------------------------------------------------------- 
    232275      !!                  ***  ROUTINE adv_umx  *** 
     
    249292      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u 
    250293      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity 
    251       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pt             ! tracer field 
     294      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   pt             ! tracer field 
    252295      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptc            ! tracer content field 
    253       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   pua_ups, pva_ups ! upstream u*a fluxes or u 
    254       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pa_ups         ! concentration advected upstream 
    255296      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes 
    256297      ! 
     
    264305      !  upstream (_ups) advection with initial mass fluxes 
    265306      ! --------------------------------------------------- 
    266       IF( ll_rachid )    zfu_ups=0.; zfv_ups=0. 
     307      IF( ll_clem )    zfu_ups=0.; zfv_ups=0. 
     308 
     309      IF( ll_gurvan .AND. pamsk==0. ) THEN 
     310         DO jj = 2, jpjm1 
     311            DO ji = fs_2, fs_jpim1 
     312               pt(ji,jj) = ( pt (ji,jj) + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) & 
     313                  &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     314            END DO 
     315         END DO 
     316         CALL lbc_lnk( pt, 'T', 1. ) 
     317      ENDIF 
     318 
     319       
    267320      IF( .NOT. ll_upsxy ) THEN 
    268321 
     
    270323         DO jj = 1, jpjm1 
    271324            DO ji = 1, fs_jpim1 
    272                zfu_ups(ji,jj) = MAX( pua_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pua_ups(ji,jj), 0._wp ) * pt(ji+1,jj) 
    273                zfv_ups(ji,jj) = MAX( pva_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pva_ups(ji,jj), 0._wp ) * pt(ji,jj+1) 
     325               IF( ll_clem ) THEN 
     326                  zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj) 
     327                  zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1) 
     328               ELSE 
     329                  zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj) 
     330                  zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1) 
     331               ENDIF 
    274332            END DO 
    275333         END DO 
     
    277335      ELSE 
    278336         ! 1 if advection of A 
    279          ! z1_a_ups already defined IF advection of other tracers 
    280          IF( pamsk == 1. )   z1_a_ups(:,:) = 1._wp 
     337         ! z1_ai already defined IF advection of other tracers 
     338         IF( pamsk == 1. )   z1_ai(:,:) = 1._wp 
    281339         ! 
    282340         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     
    284342            DO jj = 1, jpjm1 
    285343               DO ji = 1, fs_jpim1 
    286                   zfu_ups(ji,jj) = MAX( pua_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pua_ups(ji,jj), 0._wp ) * pt(ji+1,jj) 
     344                  IF( ll_clem ) THEN 
     345                     zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj) 
     346                  ELSE 
     347                     zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj) 
     348                  ENDIF 
    287349               END DO 
    288350            END DO 
     351             
    289352            ! first guess of tracer content from u-flux 
    290353            DO jj = 2, jpjm1 
    291354               DO ji = fs_2, fs_jpim1   ! vector opt. 
    292                   zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) & 
    293                      &         * tmask(ji,jj,1) * z1_a_ups(ji,jj) 
     355                  IF( ll_clem ) THEN 
     356                     IF( ll_gurvan ) THEN 
     357                        zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     358                     ELSE 
     359                        zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     360                           &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
     361                           &         * tmask(ji,jj,1) 
     362                     ENDIF 
     363                  ELSE 
     364                     zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) & 
     365                        &         * tmask(ji,jj,1) 
     366                  ENDIF 
     367                  IF( ji==26 .AND. jj==86) THEN 
     368                     WRITE(numout,*) '************************' 
     369                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 
     370                  ENDIF 
    294371               END DO 
    295372            END DO 
     
    299376            DO jj = 1, jpjm1 
    300377               DO ji = 1, fs_jpim1 
    301                   zfv_ups(ji,jj) = MAX( pva_ups(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pva_ups(ji,jj), 0._wp ) * zpt(ji,jj+1) 
     378                  IF( ll_clem ) THEN 
     379                     zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1) 
     380                  ELSE 
     381                     zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1) 
     382                  ENDIF 
    302383               END DO 
    303384            END DO 
    304             ! 
     385 
     386         ! 
    305387         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
    306388            ! flux in y-direction 
    307389            DO jj = 1, jpjm1 
    308390               DO ji = 1, fs_jpim1 
    309                   zfv_ups(ji,jj) = MAX( pva_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pva_ups(ji,jj), 0._wp ) * pt(ji,jj+1) 
     391                  IF( ll_clem ) THEN 
     392                     zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1) 
     393                  ELSE 
     394                     zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1) 
     395                  ENDIF 
    310396               END DO 
    311397            END DO 
     398 
    312399            ! first guess of tracer content from v-flux 
    313400            DO jj = 2, jpjm1 
    314401               DO ji = fs_2, fs_jpim1   ! vector opt. 
    315                   zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) & 
    316                      &         * tmask(ji,jj,1) * z1_a_ups(ji,jj)  
    317                END DO 
     402                  IF( ll_clem ) THEN 
     403                     IF( ll_gurvan ) THEN 
     404                        zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     405                     ELSE 
     406                        zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 
     407                        &                        + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
     408                        &         * tmask(ji,jj,1) 
     409                     ENDIF 
     410                  ELSE 
     411                     zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) & 
     412                        &         * tmask(ji,jj,1) 
     413                  ENDIF 
     414                  IF( ji==26 .AND. jj==86) THEN 
     415                     WRITE(numout,*) '************************' 
     416                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 
     417                  ENDIF 
     418                END DO 
    318419            END DO 
    319420            CALL lbc_lnk( zpt, 'T', 1. ) 
     
    322423            DO jj = 1, jpjm1 
    323424               DO ji = 1, fs_jpim1 
    324                   zfu_ups(ji,jj) = MAX( pua_ups(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pua_ups(ji,jj), 0._wp ) * zpt(ji+1,jj) 
     425                  IF( ll_clem ) THEN 
     426                     zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj) 
     427                  ELSE 
     428                     zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj) 
     429                  ENDIF 
    325430               END DO 
    326431            END DO 
     
    330435      ENDIF 
    331436 
    332       ! Rachid trick 
    333       ! ------------ 
    334       IF( ll_rachid .AND. kn_limiter /= 1 )  & 
    335          & CALL ctl_stop( 'STOP', 'icedyn_adv_umx: ll_rachid incompatible with limiters other than nonosc' ) 
    336       IF( ll_rachid ) THEN 
    337          call lbc_lnk (zfu_ups,'U',-1.) 
    338          call lbc_lnk (zfv_ups,'V',-1.) 
    339       IF( pamsk == 0. ) THEN 
    340          WHERE( ABS( puc(:,:) ) > 0._wp )   ;   zfu_ups(:,:) = zfu_ups(:,:) / puc(:,:) * pu(:,:) 
    341          ELSEWHERE                          ;   zfu_ups(:,:) = 0._wp 
    342          END WHERE 
    343  
    344          WHERE( ABS( pvc(:,:) ) > 0._wp )   ;   zfv_ups(:,:) = zfv_ups(:,:) / pvc(:,:) * pv(:,:) 
    345          ELSEWHERE                          ;   zfv_ups(:,:) = 0._wp 
    346          END WHERE          
    347       ENDIF 
     437      IF( ll_clem .AND. kn_limiter /= 1 )  & 
     438         & CALL ctl_stop( 'STOP', 'icedyn_adv_umx: ll_clem incompatible with limiters other than nonosc' ) 
     439 
     440      IF( ll_zeroup2 ) THEN 
     441         DO jj = 1, jpjm1 
     442            DO ji = 1, fs_jpim1   ! vector opt. 
     443               IF( amaxu(ji,jj) == 0._wp )   zfu_ups(ji,jj) = 0._wp 
     444               IF( amaxv(ji,jj) == 0._wp )   zfv_ups(ji,jj) = 0._wp 
     445            END DO 
     446         END DO 
    348447      ENDIF 
    349448 
     
    353452            ztra          = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   & 
    354453               &                + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1) ) * r1_e1e2t(ji,jj) 
    355             IF( ll_rachid ) THEN 
    356                zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
     454            IF( ll_clem ) THEN 
     455               IF( ll_gurvan ) THEN 
     456                  zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
     457               ELSE 
     458                  zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     459                     &                                      + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     460               ENDIF 
    357461            ELSE 
    358462               zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
     463            ENDIF 
     464            IF( ji==26 .AND. jj==86) THEN 
     465               WRITE(numout,*) '**************************' 
     466               WRITE(numout,*) 'zt upstream',zt_ups(ji,jj) 
    359467            ENDIF 
    360468         END DO 
     
    377485         ! 
    378486      END SELECT 
    379           
     487 
     488      IF( ll_thickness ) THEN 
     489         ! final trend with corrected fluxes 
     490         ! ------------------------------------ 
     491         DO jj = 2, jpjm1 
     492            DO ji = fs_2, fs_jpim1 
     493               IF( ll_gurvan ) THEN 
     494                  ztra       = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj)  
     495               ELSE 
     496                  ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  &  
     497                     &           + ( pt(ji,jj) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) & 
     498                     &           + ( pt(ji,jj) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj) 
     499               ENDIF 
     500               pt(ji,jj) = ( pt(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
     501 
     502               IF( pt(ji,jj) < -epsi20 ) THEN 
     503                  WRITE(numout,*) 'T<0 ',pt(ji,jj) 
     504               ENDIF 
     505                
     506               IF( pt(ji,jj) < 0._wp .AND. pt(ji,jj) >= -epsi20 )   pt(ji,jj) = 0._wp 
     507                
     508               IF( ji==26 .AND. jj==86) THEN 
     509                  WRITE(numout,*) 'zt high order',pt(ji,jj) 
     510               ENDIF 
     511            END DO 
     512         END DO 
     513         CALL lbc_lnk( pt, 'T',  1. ) 
     514      ENDIF 
     515    
     516      ! Rachid trick 
     517      ! ------------ 
     518     IF( ll_clem ) THEN 
     519         IF( pamsk == 0. ) THEN 
     520            DO jj = 1, jpjm1 
     521               DO ji = 1, fs_jpim1 
     522                  IF( ABS( puc(ji,jj) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN 
     523                     zfu_ho (ji,jj) = zfu_ho (ji,jj) * puc(ji,jj) / pu(ji,jj) 
     524                     zfu_ups(ji,jj) = zfu_ups(ji,jj) * puc(ji,jj) / pu(ji,jj) 
     525                  ELSE 
     526                     zfu_ho (ji,jj) = 0._wp 
     527                     zfu_ups(ji,jj) = 0._wp 
     528                  ENDIF 
     529                  ! 
     530                  IF( ABS( pvc(ji,jj) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN 
     531                     zfv_ho (ji,jj) = zfv_ho (ji,jj) * pvc(ji,jj) / pv(ji,jj) 
     532                     zfv_ups(ji,jj) = zfv_ups(ji,jj) * pvc(ji,jj) / pv(ji,jj) 
     533                  ELSE 
     534                     zfv_ho (ji,jj) = 0._wp   
     535                     zfv_ups(ji,jj) = 0._wp   
     536                  ENDIF 
     537               ENDDO 
     538            ENDDO 
     539         ENDIF 
     540      ENDIF 
     541 
     542      IF( ll_zeroup5 ) THEN 
     543         DO jj = 2, jpjm1 
     544            DO ji = 2, fs_jpim1   ! vector opt. 
     545               zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     546                  &                      - ( zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  ) * tmask(ji,jj,1) 
     547               IF( zpt(ji,jj) < 0. ) THEN 
     548                  zfu_ho(ji,jj)   = zfu_ups(ji,jj) 
     549                  zfu_ho(ji-1,jj) = zfu_ups(ji-1,jj) 
     550                  zfv_ho(ji,jj)   = zfv_ups(ji,jj) 
     551                  zfv_ho(ji,jj-1) = zfv_ups(ji,jj-1) 
     552               ENDIF 
     553            END DO 
     554         END DO 
     555         CALL lbc_lnk_multi( zfu_ho, 'U',  -1., zfv_ho, 'V',  -1. ) 
     556      ENDIF 
     557 
    380558      ! output upstream trend of concentration and high order fluxes 
    381559      ! ------------------------------------------------------------ 
    382       IF( pamsk == 1. ) THEN 
    383          ! upstream trend of concentration 
    384          pa_ups(:,:) = zt_ups(:,:) 
    385          ! upstream and high order u*a 
     560      IF( ll_dens ) THEN 
     561         ! high order u*a 
    386562         DO jj = 1, jpjm1 
    387563            DO ji = 1, fs_jpim1 
    388                pua_ups(ji,jj) = zfu_ups(ji,jj) 
    389                pva_ups(ji,jj) = zfv_ups(ji,jj) 
    390564               pua_ho (ji,jj) = zfu_ho (ji,jj) 
    391565               pva_ho (ji,jj) = zfv_ho (ji,jj) 
     
    395569         !!CALL lbc_lnk( pva_ho, 'V',  -1. ) 
    396570      ENDIF 
     571 
     572 
     573      IF( .NOT.ll_thickness ) THEN 
     574         ! final trend with corrected fluxes 
     575         ! ------------------------------------ 
     576         DO jj = 2, jpjm1 
     577            DO ji = fs_2, fs_jpim1  
     578               ztra = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  &  !        Div(uaH) or        Div(ua) 
     579                  &     * r1_e1e2t(ji,jj) * pdt   
     580 
     581               !!IF( ptc(ji,jj)+ztra < 0._wp ) THEN 
     582               !!   ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) )  &  !        Div(uaH) or        Div(ua) 
     583               !!      &     * r1_e1e2t(ji,jj) * pdt                  
     584               !!ENDIF 
     585               !!IF( ptc(ji,jj)+ztra < 0._wp ) THEN 
     586               !!   WRITE(numout,*) 'Tc<0 ',ptc(ji,jj)+ztra 
     587               !!   ztra = 0._wp 
     588               !!ENDIF 
     589 
     590               ptc(ji,jj) = ( ptc(ji,jj) + ztra ) * tmask(ji,jj,1) 
     591                              
     592               IF( ji==26 .AND. jj==86) THEN 
     593                  WRITE(numout,*) 'ztc high order',ptc(ji,jj) 
     594               ENDIF 
     595                
     596            END DO 
     597         END DO 
     598         CALL lbc_lnk( ptc, 'T',  1. ) 
     599      ENDIF 
    397600       
    398       ! final trend with corrected fluxes 
    399       ! ------------------------------------ 
    400       DO jj = 2, jpjm1 
    401          DO ji = fs_2, fs_jpim1  
    402             ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  &  !        Div(uaH) or        Div(ua) 
    403                &         ) * r1_e1e2t(ji,jj)   
    404             ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
    405          END DO 
    406       END DO 
    407       CALL lbc_lnk( ptc, 'T',  1. ) 
    408601      ! 
    409602   END SUBROUTINE adv_umx 
     
    448641         END DO 
    449642         IF    ( kn_limiter == 1 ) THEN 
    450             CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     643            IF( ll_clem ) THEN 
     644               CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     645            ELSE 
     646               CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     647            ENDIF 
    451648         ELSEIF( kn_limiter == 2 ) THEN 
    452649            CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     
    459656      ELSE                    !-- alternate directions --! 
    460657         ! 
     658         IF( pamsk == 1. )   z1_ai(:,:) = 1._wp 
     659         ! 
    461660         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
    462661            ! 
     
    473672            DO jj = 2, jpjm1 
    474673               DO ji = fs_2, fs_jpim1   ! vector opt. 
    475                   zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)  
     674                  IF( ll_clem ) THEN 
     675                     zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) * z1_ai(ji,jj) ) * tmask(ji,jj,1) 
     676                  ELSE                      
     677                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * z1_ai(ji,jj) 
     678                  ENDIF 
    476679               END DO 
    477680            END DO 
     
    481684            DO jj = 1, jpjm1 
    482685               DO ji = 1, fs_jpim1 
    483                   pfv_ho(ji,jj) = 0.5 * pv(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) ) 
     686                  pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) ) 
    484687               END DO 
    485688            END DO 
     
    501704            DO jj = 2, jpjm1 
    502705               DO ji = fs_2, fs_jpim1   ! vector opt. 
    503                   zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)  
     706                  IF( ll_clem ) THEN 
     707                     zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) * z1_ai(ji,jj) ) * tmask(ji,jj,1) 
     708                  ELSE 
     709                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * z1_ai(ji,jj) 
     710                  ENDIF 
    504711               END DO 
    505712            END DO 
     
    509716            DO jj = 1, jpjm1 
    510717               DO ji = 1, fs_jpim1 
    511                   pfu_ho(ji,jj) = 0.5 * pu(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) ) 
     718                  pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) ) 
    512719               END DO 
    513720            END DO 
     
    516723 
    517724         ENDIF 
    518          IF( kn_limiter == 1 )   CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     725         IF( ll_clem ) THEN 
     726            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 ) 
     727         ELSE 
     728            IF( kn_limiter == 1 )   CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     729         ENDIF 
    519730          
    520731      ENDIF 
     
    564775         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
    565776         !                                                        !--  advective form update in zzt  --! 
    566          DO jj = 2, jpjm1 
    567             DO ji = fs_2, fs_jpim1   ! vector opt. 
    568                zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  & 
    569                   &                   - pt   (ji,jj) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk 
    570                zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 
    571             END DO 
    572          END DO 
    573          CALL lbc_lnk( zzt, 'T', 1. ) 
     777 
     778         IF( ll_1stguess_clem ) THEN 
     779 
     780            ! first guess of tracer content from u-flux 
     781            DO jj = 2, jpjm1 
     782               DO ji = fs_2, fs_jpim1   ! vector opt. 
     783                  IF( ll_clem ) THEN 
     784                     IF( ll_gurvan ) THEN 
     785                        zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     786                     ELSE 
     787                        zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     788                           &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     789                     ENDIF 
     790                  ELSE 
     791                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     792                  ENDIF 
     793               END DO 
     794            END DO 
     795            CALL lbc_lnk( zzt, 'T', 1. ) 
     796 
     797         ELSE 
     798 
     799            DO jj = 2, jpjm1 
     800               DO ji = fs_2, fs_jpim1   ! vector opt. 
     801                  IF( ll_gurvan ) THEN 
     802                     zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  & 
     803                        &                   - pt   (ji,jj) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) 
     804                  ELSE 
     805                     zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  & 
     806                        &                   - pt   (ji,jj) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk 
     807                  ENDIF 
     808                  zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 
     809               END DO 
     810            END DO 
     811            CALL lbc_lnk( zzt, 'T', 1. ) 
     812         ENDIF 
     813         ! 
    574814         !                                                        !--  ultimate interpolation of pt at v-point  --! 
    575          CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho ) 
     815         IF( ll_hoxy ) THEN 
     816            CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho ) 
     817         ELSE 
     818            CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 
     819         ENDIF 
    576820         !                                                        !--  limiter in y --! 
    577821         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
    578822         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     823         !          
    579824         ! 
    580825      ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     
    586831         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
    587832         !                                                        !--  advective form update in zzt  --! 
    588          DO jj = 2, jpjm1 
    589             DO ji = fs_2, fs_jpim1 
    590                zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  & 
    591                   &                    - pt  (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk 
    592                zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 
    593             END DO 
    594          END DO 
    595          CALL lbc_lnk( zzt, 'T', 1. ) 
     833         IF( ll_1stguess_clem ) THEN 
     834             
     835            ! first guess of tracer content from v-flux  
     836            DO jj = 2, jpjm1 
     837               DO ji = fs_2, fs_jpim1   ! vector opt. 
     838                  IF( ll_clem ) THEN 
     839                     IF( ll_gurvan ) THEN 
     840                        zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     841                     ELSE 
     842                        zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 
     843                           &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     844                     ENDIF 
     845                  ELSE 
     846                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) & 
     847                        &         * tmask(ji,jj,1) 
     848                  ENDIF 
     849                END DO 
     850            END DO 
     851            CALL lbc_lnk( zzt, 'T', 1. ) 
     852             
     853         ELSE 
     854             
     855            DO jj = 2, jpjm1 
     856               DO ji = fs_2, fs_jpim1 
     857                  IF( ll_gurvan ) THEN 
     858                     zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  & 
     859                        &                   - pt   (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) 
     860                  ELSE 
     861                     zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  & 
     862                        &                   - pt   (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk 
     863                  ENDIF 
     864                  zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 
     865               END DO 
     866            END DO 
     867            CALL lbc_lnk( zzt, 'T', 1. ) 
     868         ENDIF 
     869         ! 
    596870         !                                                        !--  ultimate interpolation of pt at u-point  --! 
    597          CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho ) 
     871         IF( ll_hoxy ) THEN 
     872            CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho ) 
     873         ELSE 
     874            CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 
     875         ENDIF 
    598876         !                                                        !--  limiter in x --! 
    599877         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
    600878         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
    601879         ! 
     880         ! 
    602881      ENDIF 
     882 
     883      
    603884      IF( kn_limiter == 1 ) THEN 
    604885         IF( .NOT. ll_limiter_it2 ) THEN 
    605             CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     886            IF( ll_clem ) THEN 
     887               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     888            ELSE 
     889               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     890            ENDIF 
    606891         ELSE 
    607892            zzfu_ho(:,:) = pfu_ho(:,:) 
    608893            zzfv_ho(:,:) = pfv_ho(:,:) 
    609894            ! 1st iteration of nonosc (limit the flux with the upstream solution) 
    610             CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 
     895            IF( ll_clem ) THEN 
     896               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 
     897            ELSE 
     898               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 
     899            ENDIF 
    611900            ! guess after content field with high order 
    612901            DO jj = 2, jpjm1 
     
    618907            CALL lbc_lnk( zzt, 'T', 1. ) 
    619908            ! 2nd iteration of nonosc (limit the flux with the limited high order solution) 
    620             CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 
     909            IF( ll_clem ) THEN 
     910               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 
     911            ELSE 
     912               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 
     913            ENDIF 
    621914         ENDIF 
    622915      ENDIF 
     
    678971      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21) 
    679972         !         
    680          DO jj = 2, jpjm1 
     973         DO jj = 1, jpjm1 
    681974            DO ji = 1, fs_jpim1   ! vector opt. 
    682975               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                              pt(ji+1,jj) + pt(ji,jj)   & 
     
    687980      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23) 
    688981         ! 
    689          DO jj = 2, jpjm1 
     982         DO jj = 1, jpjm1 
    690983            DO ji = 1, fs_jpim1   ! vector opt. 
    691984               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     
    697990      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24) 
    698991         ! 
    699          DO jj = 2, jpjm1 
     992         DO jj = 1, jpjm1 
    700993            DO ji = 1, fs_jpim1   ! vector opt. 
    701994               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     
    7111004      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27) 
    7121005         ! 
    713          DO jj = 2, jpjm1 
     1006         DO jj = 1, jpjm1 
    7141007            DO ji = 1, fs_jpim1   ! vector opt. 
    7151008               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     
    7251018      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29) 
    7261019         ! 
    727          DO jj = 2, jpjm1 
     1020         DO jj = 1, jpjm1 
    7281021            DO ji = 1, fs_jpim1   ! vector opt. 
    7291022               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     
    7421035      END SELECT 
    7431036      !                                                     !-- High order flux in i-direction  --! 
     1037      IF( ll_neg ) THEN 
     1038         DO jj = 1, jpjm1 
     1039            DO ji = 1, fs_jpim1 
     1040               IF( pt_u(ji,jj) < 0._wp ) THEN 
     1041                  pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                              pt(ji+1,jj) + pt(ji,jj)   & 
     1042                     &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
     1043               ENDIF 
     1044            END DO 
     1045         END DO 
     1046      ENDIF 
     1047 
    7441048      DO jj = 1, jpjm1 
    7451049         DO ji = 1, fs_jpim1   ! vector opt. 
    746             pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj) 
     1050            IF( ll_clem ) THEN 
     1051               pfu_ho(ji,jj) = pu(ji,jj) * pt_u(ji,jj) 
     1052            ELSE 
     1053               pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj) 
     1054            ENDIF 
    7471055         END DO 
    7481056      END DO 
     
    8061114      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21) 
    8071115         DO jj = 1, jpjm1 
    808             DO ji = fs_2, fs_jpim1 
     1116            DO ji = 1, fs_jpim1 
    8091117               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  & 
    8101118                  &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 
     
    8141122      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23) 
    8151123         DO jj = 1, jpjm1 
    816             DO ji = fs_2, fs_jpim1 
     1124            DO ji = 1, fs_jpim1 
    8171125               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    8181126               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (        ( pt(ji,jj+1) + pt(ji,jj) )  & 
     
    8241132      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24) 
    8251133         DO jj = 1, jpjm1 
    826             DO ji = fs_2, fs_jpim1 
     1134            DO ji = 1, fs_jpim1 
    8271135               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    8281136               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     
    8371145      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27) 
    8381146         DO jj = 1, jpjm1 
    839             DO ji = fs_2, fs_jpim1 
     1147            DO ji = 1, fs_jpim1 
    8401148               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    8411149               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     
    8501158      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29) 
    8511159         DO jj = 1, jpjm1 
    852             DO ji = fs_2, fs_jpim1 
     1160            DO ji = 1, fs_jpim1 
    8531161               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    8541162               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     
    8661174      END SELECT 
    8671175      !                                                     !-- High order flux in j-direction  --! 
     1176      IF( ll_neg ) THEN 
     1177         DO jj = 1, jpjm1 
     1178            DO ji = 1, fs_jpim1 
     1179               IF( pt_v(ji,jj) < 0._wp ) THEN 
     1180                  pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  & 
     1181                     &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 
     1182               ENDIF 
     1183            END DO 
     1184         END DO 
     1185      ENDIF 
     1186 
    8681187      DO jj = 1, jpjm1 
    8691188         DO ji = 1, fs_jpim1   ! vector opt. 
    870             pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj) 
     1189            IF( ll_clem ) THEN 
     1190               pfv_ho(ji,jj) = pv(ji,jj) * pt_v(ji,jj) 
     1191            ELSE 
     1192               pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj) 
     1193            ENDIF 
    8711194         END DO 
    8721195      END DO 
     
    8751198      
    8761199 
    877    SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho ) 
     1200   SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho ) 
    8781201      !!--------------------------------------------------------------------- 
    8791202      !!                    ***  ROUTINE nonosc  *** 
     
    8831206      !! 
    8841207      !! **  Method  :   ... ??? 
    885       !!       warning : ptc and pt_low must be masked, but the boundaries 
     1208      !!       warning : pt and pt_low must be masked, but the boundaries 
    8861209      !!       conditions on the fluxes are not necessary zalezak (1979) 
    8871210      !!       drange (1995) multi-dimensional forward-in-time and upstream- 
     
    8941217      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pv               ! ice j-velocity => v*e1 
    8951218      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pvc              ! ice j-velocity *A => v*e1*a 
    896       REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt_low      ! before field & upstream guess of after field 
     1219      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt, pt_low  ! before field & upstream guess of after field 
    8971220      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_low, pfu_low ! upstream flux 
    8981221      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux 
    8991222      ! 
    9001223      INTEGER  ::   ji, jj    ! dummy loop indices 
    901       REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt   ! local scalars 
    902       REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign        !   -      - 
    903       REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zti_low, ztj_low 
     1224      REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt, zpos2, zneg2   ! local scalars 
     1225      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign, zcoef        !   -      - 
     1226      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zti_low, ztj_low, zzt 
    9041227      !!---------------------------------------------------------------------- 
    9051228      zbig = 1.e+40_wp 
    9061229      zsml = epsi20 
    9071230 
    908       ! Rachid trick (upstream already done in macho) 
    909       ! ------------ 
    910       IF( ll_rachid ) THEN 
    911       IF( pamsk == 0. ) THEN 
    912          WHERE( ABS( puc(:,:) ) > 0._wp )   ;   pfu_ho(:,:) = pfu_ho(:,:) / puc(:,:) * pu(:,:) 
    913          ELSEWHERE                          ;   pfu_ho(:,:) = 0._wp 
    914          END WHERE 
    915  
    916          WHERE( ABS( pvc(:,:) ) > 0._wp )   ;   pfv_ho(:,:) = pfv_ho(:,:) / pvc(:,:) * pv(:,:) 
    917          ELSEWHERE                          ;   pfv_ho(:,:) = 0._wp 
    918          END WHERE          
     1231      IF( ll_zeroup2 ) THEN 
     1232         DO jj = 1, jpjm1 
     1233            DO ji = 1, fs_jpim1   ! vector opt. 
     1234               IF( amaxu(ji,jj) == 0._wp )   pfu_ho(ji,jj) = 0._wp 
     1235               IF( amaxv(ji,jj) == 0._wp )   pfv_ho(ji,jj) = 0._wp 
     1236            END DO 
     1237         END DO 
    9191238      ENDIF 
     1239       
     1240      IF( ll_zeroup4 ) THEN 
     1241         DO jj = 1, jpjm1 
     1242            DO ji = 1, fs_jpim1   ! vector opt. 
     1243               IF( pfu_low(ji,jj) == 0._wp )   pfu_ho(ji,jj) = 0._wp 
     1244               IF( pfv_low(ji,jj) == 0._wp )   pfv_ho(ji,jj) = 0._wp 
     1245            END DO 
     1246         END DO 
    9201247      ENDIF 
     1248 
     1249 
     1250      IF( ll_zeroup1 ) THEN 
     1251         DO jj = 2, jpjm1 
     1252            DO ji = fs_2, fs_jpim1 
     1253               IF( ll_gurvan ) THEN 
     1254                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1255                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1256               ELSE 
     1257                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1258                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  & 
     1259                     &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1260                     &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     1261               ENDIF 
     1262               IF( zzt(ji,jj) < 0._wp ) THEN 
     1263                  pfu_ho(ji,jj)   = pfu_low(ji,jj) 
     1264                  pfv_ho(ji,jj)   = pfv_low(ji,jj) 
     1265                  WRITE(numout,*) '*** 1 negative high order zzt ***',ji,jj,zzt(ji,jj) 
     1266               ENDIF 
     1267               IF( ji==26 .AND. jj==86) THEN 
     1268                  WRITE(numout,*) 'zzt high order',zzt(ji,jj) 
     1269                  WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1270                  WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1271                  WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt 
     1272                  WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt 
     1273               ENDIF 
     1274               IF( ll_gurvan ) THEN 
     1275                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1276                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1277               ELSE 
     1278                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1279                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  & 
     1280                     &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1281                     &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     1282               ENDIF 
     1283               IF( zzt(ji,jj) < 0._wp ) THEN 
     1284                  pfu_ho(ji-1,jj) = pfu_low(ji-1,jj) 
     1285                  pfv_ho(ji,jj-1) = pfv_low(ji,jj-1) 
     1286                  WRITE(numout,*) '*** 2 negative high order zzt ***',ji,jj,zzt(ji,jj) 
     1287               ENDIF 
     1288               IF( ll_gurvan ) THEN 
     1289                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1290                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1291               ELSE 
     1292                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1293                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  & 
     1294                     &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1295                     &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     1296               ENDIF 
     1297               IF( zzt(ji,jj) < 0._wp ) THEN 
     1298                  WRITE(numout,*) '*** 3 negative high order zzt ***',ji,jj,zzt(ji,jj) 
     1299               ENDIF 
     1300            END DO 
     1301         END DO 
     1302         CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) 
     1303      ENDIF 
     1304 
    9211305       
    9221306      ! antidiffusive flux : high order minus low order 
     
    9261310            pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_low(ji,jj) 
    9271311            pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_low(ji,jj) 
    928          END DO 
     1312          END DO 
    9291313      END DO 
    9301314 
     
    10141398      ! Search local extrema 
    10151399      ! -------------------- 
    1016       ! max/min of ptc & pt_low with large negative/positive value (-/+zbig) outside ice cover 
     1400      ! max/min of pt & pt_low with large negative/positive value (-/+zbig) outside ice cover 
    10171401      DO jj = 1, jpj 
    10181402         DO ji = 1, jpi 
    1019 !!clem            IF    ( ptc(ji,jj) == 0._wp .AND. pt_low(ji,jj) == 0._wp ) THEN 
    1020             IF    ( ptc(ji,jj) <= epsi20 .AND. pt_low(ji,jj) <= epsi20 ) THEN 
     1403            IF    ( pt(ji,jj) <= 0._wp .AND. pt_low(ji,jj) <= 0._wp ) THEN 
    10211404               zbup(ji,jj) = -zbig 
    10221405               zbdo(ji,jj) =  zbig 
    1023 !!clem            ELSEIF( ptc(ji,jj) == 0._wp .AND. pt_low(ji,jj) /= 0._wp ) THEN 
    1024             ELSEIF( ptc(ji,jj) <= epsi20 .AND. pt_low(ji,jj) > epsi20 ) THEN 
     1406            ELSEIF( pt(ji,jj) <= 0._wp .AND. pt_low(ji,jj) > 0._wp ) THEN 
    10251407               zbup(ji,jj) = pt_low(ji,jj) 
    10261408               zbdo(ji,jj) = pt_low(ji,jj) 
    1027 !!clem            ELSEIF( ptc(ji,jj) /= 0._wp .AND. pt_low(ji,jj) == 0._wp ) THEN 
    1028             ELSEIF( ptc(ji,jj) > epsi20 .AND. pt_low(ji,jj) <= epsi20 ) THEN 
    1029                zbup(ji,jj) = ptc(ji,jj) 
    1030                zbdo(ji,jj) = ptc(ji,jj) 
     1409            ELSEIF( pt(ji,jj) > 0._wp .AND. pt_low(ji,jj) <= 0._wp ) THEN 
     1410               zbup(ji,jj) = pt(ji,jj) 
     1411               zbdo(ji,jj) = pt(ji,jj) 
    10311412            ELSE 
    1032                zbup(ji,jj) = MAX( ptc(ji,jj) , pt_low(ji,jj) ) 
    1033                zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_low(ji,jj) ) 
    1034             ENDIF 
    1035          END DO 
    1036       END DO 
    1037  
     1413               zbup(ji,jj) = MAX( pt(ji,jj) , pt_low(ji,jj) ) 
     1414               zbdo(ji,jj) = MIN( pt(ji,jj) , pt_low(ji,jj) ) 
     1415            ENDIF 
     1416         END DO 
     1417      END DO 
     1418 
     1419  
    10381420      z1_dt = 1._wp / pdt 
    10391421      DO jj = 2, jpjm1 
    10401422         DO ji = fs_2, fs_jpim1   ! vector opt. 
    10411423            ! 
    1042 !            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 
    1043 !            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1) ) 
    1044             ! 
     1424            IF( .NOT. ll_9points ) THEN 
     1425            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 
     1426            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1) ) 
     1427            ! 
     1428            ELSE 
    10451429            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 
    10461430               &                     zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1)  ) 
    10471431            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1), & 
    10481432               &                     zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1)  ) 
     1433            ENDIF 
    10491434            ! 
    10501435            zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji  ,jj) ) &  ! positive/negative part of the flux 
     
    10531438               & + MAX( 0., pfv_ho(ji,jj  ) ) - MIN( 0., pfv_ho(ji,jj-1) ) 
    10541439            ! 
     1440            IF( ll_HgradU .AND. .NOT.ll_gurvan ) THEN 
     1441               zneg2 = (   pt(ji,jj) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) ) * ( 1. - pamsk ) 
     1442               zpos2 = ( - pt(ji,jj) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) - pt(ji,jj) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) ) * ( 1. - pamsk ) 
     1443            ELSE 
     1444               zneg2 = 0. ; zpos2 = 0. 
     1445            ENDIF 
     1446            ! 
    10551447            !                                  ! up & down beta terms 
    10561448!            zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt 
    10571449!            zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt 
    1058             IF( zpos >= epsi10 ) THEN ; zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt 
    1059             ELSE                      ; zbetup(ji,jj) = 0. 
    1060             ENDIF 
    1061             ! 
    1062             IF( zneg >= epsi10 ) THEN ; zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 
    1063             ELSE                      ; zbetdo(ji,jj) = 0. 
    1064             ENDIF 
    1065             ! 
    1066             IF( zbetdo(ji,jj) < 0._wp ) zbetdo(ji,jj)=0. 
    1067             IF( zbetup(ji,jj) < 0._wp ) zbetup(ji,jj)=0. 
    1068             ! 
     1450 
     1451            IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj) = MAX( 0._wp, zup - pt_low(ji,jj) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt 
     1452            ELSE                         ; zbetup(ji,jj) = 0. ! zbig 
     1453            ENDIF 
     1454            ! 
     1455            IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj) = MAX( 0._wp, pt_low(ji,jj) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt 
     1456            ELSE                         ; zbetdo(ji,jj) = 0. ! zbig 
     1457            ENDIF 
     1458            ! 
     1459            ! if all the points are outside ice cover 
     1460            IF( zup == -zbig )   zbetup(ji,jj) = 0. ! zbig 
     1461            IF( zdo ==  zbig )   zbetdo(ji,jj) = 0. ! zbig             
     1462            ! 
     1463 
     1464            IF( ji==26 .AND. jj==86) THEN 
     1465               WRITE(numout,*) '-----------------' 
     1466               WRITE(numout,*) 'zpos',zpos,zpos2 
     1467               WRITE(numout,*) 'zneg',zneg,zneg2 
     1468               WRITE(numout,*) 'puc/pu',ABS(puc(ji,jj))/MAX(epsi20, ABS(pu(ji,jj))) 
     1469               WRITE(numout,*) 'pvc/pv',ABS(pvc(ji,jj))/MAX(epsi20, ABS(pv(ji,jj))) 
     1470               WRITE(numout,*) 'pucm1/pu',ABS(puc(ji-1,jj))/MAX(epsi20, ABS(pu(ji-1,jj))) 
     1471               WRITE(numout,*) 'pvcm1/pv',ABS(pvc(ji,jj-1))/MAX(epsi20, ABS(pv(ji,jj-1))) 
     1472               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)+pfu_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1473               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)+pfv_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1474               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)+pfu_low(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt 
     1475               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)+pfv_low(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt 
     1476               WRITE(numout,*) 'pfu_low',pfu_low(ji,jj) * r1_e1e2t(ji,jj) * pdt 
     1477               WRITE(numout,*) 'pfv_low',pfv_low(ji,jj) * r1_e1e2t(ji,jj) * pdt 
     1478               WRITE(numout,*) 'pfu_lowm1',pfu_low(ji-1,jj) * r1_e1e2t(ji,jj) * pdt 
     1479               WRITE(numout,*) 'pfv_lowm1',pfv_low(ji,jj-1) * r1_e1e2t(ji,jj) * pdt 
     1480                
     1481               WRITE(numout,*) 'pt',pt(ji,jj) 
     1482               WRITE(numout,*) 'ptim1',pt(ji-1,jj) 
     1483               WRITE(numout,*) 'ptjm1',pt(ji,jj-1) 
     1484               WRITE(numout,*) 'pt_low',pt_low(ji,jj) 
     1485               WRITE(numout,*) 'zbetup',zbetup(ji,jj) 
     1486               WRITE(numout,*) 'zbetdo',zbetdo(ji,jj) 
     1487               WRITE(numout,*) 'zup',zup 
     1488               WRITE(numout,*) 'zdo',zdo 
     1489            ENDIF 
    10691490            ! 
    10701491         END DO 
     
    10721493      CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign) 
    10731494 
    1074       ! Rachid trick 
    1075       ! ------------ 
    1076       IF( ll_rachid ) THEN 
    1077       IF( pamsk == 0. ) THEN 
    1078          WHERE( ABS( pu(:,:) ) > 0._wp ) 
    1079             pfu_ho (:,:) = pfu_ho (:,:) * puc(:,:) / pu(:,:) 
    1080             pfu_low(:,:) = pfu_low(:,:) * puc(:,:) / pu(:,:) 
    1081          ELSEWHERE 
    1082             pfu_ho (:,:) = 0._wp 
    1083             pfu_low(:,:) = 0._wp 
    1084          END WHERE 
    1085  
    1086          WHERE( ABS( pv(:,:) ) > 0._wp ) 
    1087             pfv_ho (:,:) = pfv_ho (:,:) * pvc(:,:) / pv(:,:) 
    1088             pfv_low(:,:) = pfv_low(:,:) * pvc(:,:) / pv(:,:) 
    1089          ELSEWHERE 
    1090             pfv_ho (:,:) = 0._wp 
    1091             pfv_low(:,:) = 0._wp 
    1092          END WHERE          
    1093       ENDIF 
    1094       ENDIF 
    10951495       
    10961496      ! monotonic flux in the y direction 
     
    11021502            zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj) ) 
    11031503            ! 
    1104             pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_low(ji,jj) 
     1504            zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu ) 
     1505             
     1506            pfu_ho(ji,jj) = pfu_ho(ji,jj) * zcoef + pfu_low(ji,jj) 
     1507 
     1508            IF( ji==26 .AND. jj==86) THEN 
     1509               WRITE(numout,*) 'coefU',zcoef 
     1510               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1511               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt 
     1512            ENDIF 
     1513 
    11051514         END DO 
    11061515      END DO 
     
    11121521            zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj) ) 
    11131522            ! 
    1114             pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_low(ji,jj) 
    1115          END DO 
    1116       END DO 
     1523            zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv ) 
     1524             
     1525            pfv_ho(ji,jj) = pfv_ho(ji,jj) * zcoef + pfv_low(ji,jj) 
     1526 
     1527            IF( ji==26 .AND. jj==86) THEN 
     1528               WRITE(numout,*) 'coefV',zcoef 
     1529               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1530               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt 
     1531            ENDIF 
     1532         END DO 
     1533      END DO 
     1534 
     1535      ! clem test 
     1536      DO jj = 2, jpjm1 
     1537         DO ji = 2, fs_jpim1   ! vector opt. 
     1538            IF( ll_gurvan ) THEN 
     1539               zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1540                  &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1541            ELSE 
     1542               zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1543                  &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  & 
     1544                  &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1545                  &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     1546            ENDIF 
     1547            IF( zzt(ji,jj) < -epsi20 ) THEN 
     1548               WRITE(numout,*) 'T<0 nonosc',zzt(ji,jj) 
     1549            ENDIF 
     1550         END DO 
     1551      END DO 
     1552       
     1553      ! 
    11171554      ! 
    11181555   END SUBROUTINE nonosc_2d 
Note: See TracChangeset for help on using the changeset viewer.