New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 10399 – NEMO

Changeset 10399


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

improve ice advection (toward an acceptable solution)

Location:
NEMO/branches/2018/dev_r9947_SI3_advection
Files:
3 deleted
9 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2018/dev_r9947_SI3_advection/cfgs/SHARED/namelist_ice_ref

    r10312 r10399  
    4949&namdyn         !   Ice dynamics 
    5050!------------------------------------------------------------------------------ 
    51    ln_dynFULL       = .true.          !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction) 
    52    ln_dynRHGADV     = .false.         !  dyn.: no ridge/raft & no corrections  (rheology + advection) 
    53    ln_dynADV        = .false.         !  dyn.: only advection w prescribed vel.(rn_uvice + advection) 
     51   ln_dynALL        = .true.          !  dyn.: full ice dynamics                  (rheology + advection + ridging/rafting + correction) 
     52   ln_dynRHGADV     = .false.         !  dyn.: no ridge/raft & no corrections     (rheology + advection) 
     53   ln_dynADV1D      = .false.         !  dyn.: only advection 1D                  (Schar & Smolarkiewicz 1996 test case) 
     54   ln_dynADV2D      = .false.         !  dyn.: only advection 2D w prescribed vel.(rn_uvice + advection) 
    5455      rn_uice       =   0.5           !        prescribed ice u-velocity 
    5556      rn_vice       =   0.5           !        prescribed ice v-velocity 

  • NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/ice.F90

    r10312 r10399  
    145145   INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling 
    146146   REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp)  
     147   ! 
     148   !                                     !!** ice-advection namelist (namdyn_adv) ** 
     149   LOGICAL , PUBLIC ::   ln_adv_Pra       !: Prather        advection scheme 
     150   LOGICAL , PUBLIC ::   ln_adv_UMx       !: Ultimate-Macho advection scheme 
    147151   ! 
    148152   !                                     !!** ice-surface forcing namelist (namforcing) ** 

  • NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn.F90

    r10312 r10399  
    2121   USE icecor         ! sea-ice: corrections 
    2222   USE icevar         ! sea-ice: operations 
     23   USE icectl         ! sea-ice: control prints 
    2324   ! 
    2425   USE in_out_manager ! I/O manager 
    … …  
    3738   INTEGER ::              nice_dyn   ! choice of the type of dynamics 
    3839   !                                        ! associated indices: 
    39    INTEGER, PARAMETER ::   np_dynFULL    = 1   ! full ice dynamics               (rheology + advection + ridging/rafting + correction) 
     40   INTEGER, PARAMETER ::   np_dynALL     = 1   ! full ice dynamics               (rheology + advection + ridging/rafting + correction) 
    4041   INTEGER, PARAMETER ::   np_dynRHGADV  = 2   ! pure dynamics                   (rheology + advection)  
    41    INTEGER, PARAMETER ::   np_dynADV     = 3   ! only advection w prescribed vel.(rn_uvice + advection) 
     42   INTEGER, PARAMETER ::   np_dynADV1D   = 3   ! only advection 1D - test case from Schar & Smolarkiewicz 1996 
     43   INTEGER, PARAMETER ::   np_dynADV2D   = 4   ! only advection 2D w prescribed vel.(rn_uvice + advection) 
    4244   ! 
    4345   ! ** namelist (namdyn) ** 
    44    LOGICAL  ::   ln_dynFULL       ! full ice dynamics               (rheology + advection + ridging/rafting + correction) 
    45    LOGICAL  ::   ln_dynRHGADV     ! no ridge/raft & no corrections  (rheology + advection) 
    46    LOGICAL  ::   ln_dynADV        ! only advection w prescribed vel.(rn_uvice + advection) 
    47    REAL(wp) ::   rn_uice          !    prescribed u-vel (case np_dynADV) 
    48    REAL(wp) ::   rn_vice          !    prescribed v-vel (case np_dynADV) 
     46   LOGICAL  ::   ln_dynALL        ! full ice dynamics                      (rheology + advection + ridging/rafting + correction) 
     47   LOGICAL  ::   ln_dynRHGADV     ! no ridge/raft & no corrections         (rheology + advection) 
     48   LOGICAL  ::   ln_dynADV1D      ! only advection in 1D w ice convergence (test case from Schar & Smolarkiewicz 1996) 
     49   LOGICAL  ::   ln_dynADV2D      ! only advection in 2D w prescribed vel. (rn_uvice + advection) 
     50   REAL(wp) ::   rn_uice          !    prescribed u-vel (case np_dynADV1D & np_dynADV2D) 
     51   REAL(wp) ::   rn_vice          !    prescribed v-vel (case np_dynADV1D & np_dynADV2D) 
    4952    
    5053   !! * Substitutions 
    … …  
    7376      INTEGER  ::   ji, jj, jl        ! dummy loop indices 
    7477      REAL(wp) ::   zcoefu, zcoefv 
    75       REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhmax 
     78      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max 
    7679      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i 
    7780      !!-------------------------------------------------------------------- 
    7881      ! 
    79       IF( ln_timing )   CALL timing_start('icedyn') 
     82      ! controls 
     83      IF( ln_timing    )   CALL timing_start('icedyn')                                                             ! timing 
     84      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    8085      ! 
    8186      IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    8893         tau_icebfr(:,:) = 0._wp 
    8994         DO jl = 1, jpl 
    90             WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
    91          END DO 
    92       ENDIF 
    93       ! 
    94       zhmax(:,:,:) = h_i_b(:,:,:)      !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig) 
     95            WHERE( h_i_b(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
     96         END DO 
     97      ENDIF 
     98      ! 
     99      !                                !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig) 
    95100      DO jl = 1, jpl 
    96101         DO jj = 2, jpjm1 
    97             DO ji = 2, jpim1 
    98 !!gm use of MAXVAL here is very probably less efficient than expending the 9 values 
    99                zhmax(ji,jj,jl) = MAX( epsi20, MAXVAL( h_i_b(ji-1:ji+1,jj-1:jj+1,jl) ) ) 
     102            DO ji = fs_2, fs_jpim1 
     103               zhi_max(ji,jj,jl) = MAX( epsi20, h_i_b(ji,jj,jl), h_i_b(ji+1,jj  ,jl), h_i_b(ji  ,jj+1,jl), & 
     104                  &                                              h_i_b(ji-1,jj  ,jl), h_i_b(ji  ,jj-1,jl), & 
     105                  &                                              h_i_b(ji+1,jj+1,jl), h_i_b(ji-1,jj-1,jl), & 
     106                  &                                              h_i_b(ji+1,jj-1,jl), h_i_b(ji-1,jj+1,jl) ) 
     107               zhs_max(ji,jj,jl) = MAX( epsi20, h_s_b(ji,jj,jl), h_s_b(ji+1,jj  ,jl), h_s_b(ji  ,jj+1,jl), & 
     108                  &                                              h_s_b(ji-1,jj  ,jl), h_s_b(ji  ,jj-1,jl), & 
     109                  &                                              h_s_b(ji+1,jj+1,jl), h_s_b(ji-1,jj-1,jl), & 
     110                  &                                              h_s_b(ji+1,jj-1,jl), h_s_b(ji-1,jj+1,jl) ) 
    100111            END DO 
    101112         END DO 
    102113      END DO 
    103       CALL lbc_lnk( zhmax(:,:,:), 'T', 1. ) 
     114      CALL lbc_lnk_multi( zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1. ) 
    104115      ! 
    105116      ! 
    106117      SELECT CASE( nice_dyn )           !-- Set which dynamics is running 
    107118 
    108       CASE ( np_dynFULL )          !==  all dynamical processes  ==! 
    109          CALL ice_dyn_rhg   ( kt )                            ! -- rheology   
    110          CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhmax )   ! -- advection of ice + correction on ice thickness 
    111          CALL ice_dyn_rdgrft( kt )                            ! -- ridging/rafting  
    112          CALL ice_cor       ( kt , 1 )                        ! -- Corrections 
     119      CASE ( np_dynALL )           !==  all dynamical processes  ==! 
     120         CALL ice_dyn_rhg   ( kt )                                       ! -- rheology   
     121         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max )   ! -- advection of ice + correction on ice thickness 
     122         CALL ice_dyn_rdgrft( kt )                                       ! -- ridging/rafting  
     123         CALL ice_cor       ( kt , 1 )                                   ! -- Corrections 
    113124 
    114125      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==! 
    115          CALL ice_dyn_rhg   ( kt )                            ! -- rheology   
    116          CALL ice_dyn_adv   ( kt )                            ! -- advection of ice 
    117          CALL Hpiling                                         ! -- simple pile-up (replaces ridging/rafting) 
    118  
    119       CASE ( np_dynADV )           !==  pure advection ==!   (prescribed velocities) 
     126         CALL ice_dyn_rhg   ( kt )                                       ! -- rheology   
     127         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max )   ! -- advection of ice + correction on ice thickness 
     128         CALL Hpiling                                                    ! -- simple pile-up (replaces ridging/rafting) 
     129 
     130      CASE ( np_dynADV1D )         !==  pure advection ==!   (1D) 
    120131         ALLOCATE( zdivu_i(jpi,jpj) ) 
    121  
    122          u_ice(:,:) = rn_uice * umask(:,:,1) 
    123          v_ice(:,:) = rn_vice * vmask(:,:,1) 
    124          !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1) 
    125          !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 
    126          ! --- monotonicity test from Lipscomb et al 2004 --- ! 
     132         ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- ! 
    127133         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length 
    128134         ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s  
    129          !DO jj = 1, jpj 
    130          !   DO ji = 1, jpi 
    131          !      zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. ) 
    132          !      zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. ) 
    133          !      u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1) 
    134          !      v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1) 
    135          !   END DO 
    136          !END DO 
     135         DO jj = 1, jpj 
     136            DO ji = 1, jpi 
     137               zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. ) 
     138               zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. ) 
     139               u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1) 
     140               v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1) 
     141            END DO 
     142         END DO 
    137143         ! --- 
    138          CALL ice_dyn_adv   ( kt )                            ! -- advection of ice 
     144         CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice + correction on ice thickness 
    139145 
    140146         ! diagnostics: divergence at T points  
    … …  
    150156         DEALLOCATE( zdivu_i ) 
    151157 
     158      CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities) 
     159         ALLOCATE( zdivu_i(jpi,jpj) ) 
     160         u_ice(:,:) = rn_uice * umask(:,:,1) 
     161         v_ice(:,:) = rn_vice * vmask(:,:,1) 
     162         !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1) 
     163         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 
     164         ! --- 
     165         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max )   ! -- advection of ice + correction on ice thickness 
     166 
     167         ! diagnostics: divergence at T points  
     168         DO jj = 2, jpjm1 
     169            DO ji = 2, jpim1 
     170               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     171                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     172            END DO 
     173         END DO 
     174         CALL lbc_lnk( zdivu_i, 'T', 1. ) 
     175         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) ) 
     176 
     177         DEALLOCATE( zdivu_i ) 
     178 
    152179      END SELECT 
    153       ! 
    154       IF( ln_timing )   CALL timing_stop('icedyn') 
     180       ! 
     181      ! controls 
     182      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     183      IF( ln_timing    )   CALL timing_stop ('icedyn')                                                             ! timing 
    155184      ! 
    156185   END SUBROUTINE ice_dyn 
    157186 
    158187 
    159    SUBROUTINE Hbig( phmax ) 
     188   SUBROUTINE Hbig( phi_max, phs_max ) 
    160189      !!------------------------------------------------------------------- 
    161190      !!                  ***  ROUTINE Hbig  *** 
    162191      !! 
    163192      !! ** Purpose : Thickness correction in case advection scheme creates 
    164       !!              abnormally tick ice 
    165       !! 
    166       !! ** Method  : 1- check whether ice thickness resulting from advection is 
    167       !!                 larger than the surrounding 9-points before advection 
    168       !!                 and reduce it if a) divergence or b) convergence & at_i>0.8 
    169       !!              2- bound ice thickness with hi_max (99m) 
     193      !!              abnormally tick ice or snow 
     194      !! 
     195      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points 
     196      !!                 (before advection) and reduce it by adapting ice concentration 
     197      !!              2- check whether snow thickness is larger than the surrounding 9-points 
     198      !!                 (before advection) and reduce it by sending the excess in the ocean 
     199      !!              3- check whether snow load deplets the snow-ice interface below sea level$ 
     200      !!                 and reduce it by sending the excess in the ocean 
     201      !!              4- correct pond fraction to avoid a_ip > a_i 
    170202      !! 
    171203      !! ** input   : Max thickness of the surrounding 9-points 
    172204      !!------------------------------------------------------------------- 
    173       REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phmax   ! max ice thick from surrounding 9-pts 
     205      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max   ! max ice thick from surrounding 9-pts 
    174206      ! 
    175207      INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    176       REAL(wp) ::   zh, zdv 
     208      REAL(wp) ::   zhi, zhs, zvs_excess, zfra 
    177209      !!------------------------------------------------------------------- 
    178210      ! 
    … …  
    182214         DO jj = 1, jpj 
    183215            DO ji = 1, jpi 
    184                IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to hmax 
     216               IF ( v_i(ji,jj,jl) > 0._wp ) THEN 
    185217                  ! 
    186                   zh  = v_i (ji,jj,jl) / a_i(ji,jj,jl) 
    187                   zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl)   
    188                   ! 
    189                   IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. & 
    190                      & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN 
    191                      a_i (ji,jj,jl) = v_i(ji,jj,jl) / MIN( phmax(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
     218                  !                               ! -- check h_i -- ! 
     219                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
     220                  zhi = v_i (ji,jj,jl) / a_i(ji,jj,jl) 
     221!!clem                  zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl)   
     222!!clem                  IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. & 
     223!!clem                     & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN 
     224                  IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN 
     225                     a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
    192226                  ENDIF 
    193227                  ! 
     228                  !                               ! -- check h_s -- ! 
     229                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
     230                  zhs = v_s (ji,jj,jl) / a_i(ji,jj,jl) 
     231                  IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN 
     232                     zfra = a_i(ji,jj,jl) * phs_max(ji,jj,jl) / MAX( v_s(ji,jj,jl), epsi20 ) 
     233                     ! 
     234                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice 
     235                     hfx_res(ji,jj) = hfx_res(ji,jj) - e_s(ji,jj,1,jl) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
     236                     ! 
     237                     e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zfra 
     238                     v_s(ji,jj,jl)   = a_i(ji,jj,jl) * phs_max(ji,jj,jl) 
     239                  ENDIF            
     240                  ! 
     241                  !                               ! -- check snow load -- ! 
     242                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 
     243                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 
     244                  !    this imposed mini can artificially make the snow thickness very high (if concentration decreases drastically) 
     245                  zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
     246                  IF( zvs_excess > 0._wp ) THEN 
     247                     zfra = zvs_excess / MAX( v_s(ji,jj,jl), epsi20 ) 
     248                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice 
     249                     hfx_res(ji,jj) = hfx_res(ji,jj) - e_s(ji,jj,1,jl) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
     250                     ! 
     251                     e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zfra 
     252                     v_s(ji,jj,jl)   = v_s(ji,jj,jl) - zvs_excess 
     253                  ENDIF 
     254                   
    194255               ENDIF 
    195256            END DO 
    … …  
    241302      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
    242303      !! 
    243       NAMELIST/namdyn/ ln_dynFULL, ln_dynRHGADV, ln_dynADV, rn_uice, rn_vice,  & 
    244          &             rn_ishlat ,                                             & 
     304      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  & 
     305         &             rn_ishlat ,                                                           & 
    245306         &             ln_landfast_L16, ln_landfast_home, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile 
    246307      !!------------------------------------------------------------------- 
    … …  
    259320         WRITE(numout,*) '~~~~~~~~~~~~' 
    260321         WRITE(numout,*) '   Namelist namdyn:' 
    261          WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynFULL      = ', ln_dynFULL 
     322         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynALL       = ', ln_dynALL 
    262323         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV    = ', ln_dynRHGADV 
    263          WRITE(numout,*) '      Advection only         (rn_uvice + adv)                 ln_dynADV       = ', ln_dynADV 
     324         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)     ln_dynADV1D     = ', ln_dynADV1D 
     325         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                 ln_dynADV2D     = ', ln_dynADV2D 
    264326         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)    = (', rn_uice,',', rn_vice,')' 
    265327         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat       = ', rn_ishlat 
    … …  
    275337      ioptio = 0  
    276338      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction) 
    277       IF( ln_dynFULL   ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynFULL      ;   ENDIF 
     339      IF( ln_dynALL    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynALL       ;   ENDIF 
    278340      !      !--- dynamics without ridging/rafting and corr   (rheology + advection) 
    279341      IF( ln_dynRHGADV ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV    ;   ENDIF 
    280       !      !--- advection only with prescribed ice velocities (from namelist) 
    281       IF( ln_dynADV    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV       ;   ENDIF 
     342      !      !--- advection 1D only - test case from Schar & Smolarkiewicz 1996 
     343      IF( ln_dynADV1D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV1D     ;   ENDIF 
     344      !      !--- advection 2D only with prescribed ice velocities (from namelist) 
     345      IF( ln_dynADV2D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV2D     ;   ENDIF 
    282346      ! 
    283347      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' ) 

  • NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv.F90

    r9943 r10399  
    4040   ! 
    4141   ! ** namelist (namdyn_adv) ** 
    42    LOGICAL ::   ln_adv_Pra   ! Prather        advection scheme 
    43    LOGICAL ::   ln_adv_UMx   ! Ultimate-Macho advection scheme 
    44    INTEGER ::   nn_UMx       ! order of the UMx advection scheme    
     42   INTEGER         ::   nn_UMx       ! order of the UMx advection scheme    
    4543   ! 
    4644   !! * Substitution 
    … …  
    8987      CASE( np_advUMx )                ! ULTIMATE-MACHO scheme ! 
    9088         !                             !-----------------------! 
    91          CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice,  & 
    92             &                  ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     89         CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    9390      !                                !-----------------------! 
    9491      CASE( np_advPRA )                ! PRATHER scheme        ! 
    9592         !                             !-----------------------! 
    96          CALL ice_dyn_adv_pra( kt, u_ice, v_ice,  & 
    97             &                  ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     93         CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    9894      END SELECT 
    9995 
    … …  
    10197      ! Debug the advection schemes 
    10298      !---------------------------- 
    103       ! clem: At least one advection scheme above is not strictly positive => UM from 3d to 5th order 
    104       !       In Prather, I am not sure if the fields are bounded by 0 or not (it seems not) 
    105       !       In UM3-5  , advected fields are not bounded and negative values can appear. 
     99      ! clem: At least one advection scheme above is not strictly positive => UMx 
     100      !       In Prather, I am not sure if the fields are bounded by 0 or not (it seems yes) 
     101      !       In UMx    , advected fields are not perfectly bounded and negative values can appear. 
    106102      !                   These values are usually very small but in some occasions they can also be non-negligible 
    107103      !                   Therefore one needs to bound the advected fields by 0 (though this is not a clean fix) 
    … …  
    118114      ! 
    119115      ! ==> conservation is ensured by calling this routine below, 
    120       !     however the global ice volume is then changed by advection (but errors are very small)  
     116      !     however the global ice volume is then changed by advection (but errors are small)  
    121117      CALL ice_var_zapneg( ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    122118 

  • 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 

  • NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_rdgrft.F90

    r10312 r10399  
    192192            ! divergence given by the advection scheme 
    193193            !   (which may not be equal to divu as computed from the velocity field) 
    194             zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 
     194            IF( ln_adv_Pra ) THEN 
     195               zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 
     196            ELSEIF( ln_adv_UMx ) THEN 
     197               zdivu_adv(ji) = zdivu(ji) 
     198            ENDIF 
    195199            ! 
    196200            IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )   ! make sure the closing rate is large enough 

  • NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/iceistate.F90

    r9935 r10399  
    174174         ! then we check whether the distribution fullfills 
    175175         ! volume and area conservation, positivity and ice categories bounds 
    176          zh_i_ini(:,:,:) = 0._wp  
    177          za_i_ini(:,:,:) = 0._wp 
    178          ! 
    179          DO jj = 1, jpj 
    180             DO ji = 1, jpi 
    181                ! 
    182                IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN 
    183  
    184                   ! find which category (jl0) the input ice thickness falls into 
    185                   jl0 = jpl 
    186                   DO jl = 1, jpl 
    187                      IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 
    188                         jl0 = jl 
    189                         CYCLE 
    190                      ENDIF 
    191                   END DO 
     176 
     177         IF( jpl == 1 ) THEN 
     178            ! 
     179            zh_i_ini(:,:,1) = zht_i_ini(:,:) 
     180            za_i_ini(:,:,1) = zat_i_ini(:,:)             
     181            ! 
     182         ELSE 
     183            zh_i_ini(:,:,:) = 0._wp  
     184            za_i_ini(:,:,:) = 0._wp 
     185            ! 
     186            DO jj = 1, jpj 
     187               DO ji = 1, jpi 
    192188                  ! 
    193                   itest(:) = 0 
    194                   i_fill   = jpl + 1                                            !------------------------------------ 
    195                   DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories 
    196                      !                                                          !------------------------------------ 
    197                      i_fill = i_fill - 1 
     189                  IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN 
     190 
     191                     ! find which category (jl0) the input ice thickness falls into 
     192                     jl0 = jpl 
     193                     DO jl = 1, jpl 
     194                        IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 
     195                           jl0 = jl 
     196                           CYCLE 
     197                        ENDIF 
     198                     END DO 
    198199                     ! 
    199                      zh_i_ini(ji,jj,:) = 0._wp  
    200                      za_i_ini(ji,jj,:) = 0._wp 
    201200                     itest(:) = 0 
    202                      ! 
    203                      IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1 
    204                         zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj) 
    205                         za_i_ini(ji,jj,1) = zat_i_ini(ji,jj) 
    206                      ELSE                         !-- case ice is thicker: fill categories >1 
    207                         ! thickness 
    208                         DO jl = 1, i_fill-1 
    209                            zh_i_ini(ji,jj,jl) = hi_mean(jl) 
    210                         END DO 
     201                     i_fill   = jpl + 1                                            !------------------------------------ 
     202                     DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories 
     203                        !                                                          !------------------------------------ 
     204                        i_fill = i_fill - 1 
    211205                        ! 
    212                         ! concentration 
    213                         za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 
    214                         DO jl = 1, i_fill - 1 
    215                            IF( jl /= jl0 )THEN 
    216                               zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) ) 
    217                               za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 
     206                        zh_i_ini(ji,jj,:) = 0._wp  
     207                        za_i_ini(ji,jj,:) = 0._wp 
     208                        itest(:) = 0 
     209                        ! 
     210                        IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1 
     211                           zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj) 
     212                           za_i_ini(ji,jj,1) = zat_i_ini(ji,jj) 
     213                        ELSE                         !-- case ice is thicker: fill categories >1 
     214                           ! thickness 
     215                           DO jl = 1, i_fill-1 
     216                              zh_i_ini(ji,jj,jl) = hi_mean(jl) 
     217                           END DO 
     218                           ! 
     219                           ! concentration 
     220                           za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 
     221                           DO jl = 1, i_fill - 1 
     222                              IF( jl /= jl0 )THEN 
     223                                 zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) ) 
     224                                 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 
     225                              ENDIF 
     226                           END DO 
     227 
     228                           ! last category 
     229                           za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) ) 
     230                           zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) ) 
     231                           zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 )  
     232 
     233                           ! correction if concentration of upper cat is greater than lower cat 
     234                           !   (it should be a gaussian around jl0 but sometimes it is not) 
     235                           IF ( jl0 /= jpl ) THEN 
     236                              DO jl = jpl, jl0+1, -1 
     237                                 IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN 
     238                                    zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl) 
     239                                    zh_i_ini(ji,jj,jl    ) = 0._wp 
     240                                    za_i_ini(ji,jj,jl    ) = 0._wp 
     241                                    za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  & 
     242                                       &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 ) 
     243                                 END IF 
     244                              ENDDO 
    218245                           ENDIF 
    219                         END DO 
    220  
    221                         ! last category 
    222                         za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) ) 
    223                         zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) ) 
    224                         zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 )  
    225  
    226                         ! correction if concentration of upper cat is greater than lower cat 
    227                         !   (it should be a gaussian around jl0 but sometimes it is not) 
    228                         IF ( jl0 /= jpl ) THEN 
    229                            DO jl = jpl, jl0+1, -1 
    230                               IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN 
    231                                  zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl) 
    232                                  zh_i_ini(ji,jj,jl    ) = 0._wp 
    233                                  za_i_ini(ji,jj,jl    ) = 0._wp 
    234                                  za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  & 
    235                                     &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 ) 
    236                               END IF 
    237                            ENDDO 
     246                           ! 
    238247                        ENDIF 
    239248                        ! 
     249                        ! Compatibility tests 
     250                        zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )           ! Test 1: area conservation 
     251                        IF ( zconv < epsi06 ) itest(1) = 1 
     252                        ! 
     253                        zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &         ! Test 2: volume conservation 
     254                           &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) ) 
     255                        IF ( zconv < epsi06 ) itest(2) = 1 
     256                        ! 
     257                        IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1           ! Test 3: thickness of the last category is in-bounds ? 
     258                        ! 
     259                        itest(4) = 1 
     260                        DO jl = 1, i_fill 
     261                           IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0                        ! Test 4: positivity of ice concentrations 
     262                        END DO 
     263                        !                                                          !---------------------------- 
     264                     END DO                                                        ! end iteration on categories 
     265                     !                                                             !---------------------------- 
     266                     IF( lwp .AND. SUM(itest) /= 4 ) THEN  
     267                        WRITE(numout,*) 
     268                        WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! ' 
     269                        WRITE(numout,*) ' !!!! Something is wrong in the SI3 initialization procedure ' 
     270                        WRITE(numout,*) 
     271                        WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:) 
     272                        WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 
     273                        WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 
    240274                     ENDIF 
    241275                     ! 
    242                      ! Compatibility tests 
    243                      zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )           ! Test 1: area conservation 
    244                      IF ( zconv < epsi06 ) itest(1) = 1 
    245                      ! 
    246                      zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &         ! Test 2: volume conservation 
    247                         &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) ) 
    248                      IF ( zconv < epsi06 ) itest(2) = 1 
    249                      ! 
    250                      IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1           ! Test 3: thickness of the last category is in-bounds ? 
    251                      ! 
    252                      itest(4) = 1 
    253                      DO jl = 1, i_fill 
    254                         IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0                        ! Test 4: positivity of ice concentrations 
    255                      END DO 
    256                      !                                                          !---------------------------- 
    257                   END DO                                                        ! end iteration on categories 
    258                   !                                                             !---------------------------- 
    259                   IF( lwp .AND. SUM(itest) /= 4 ) THEN  
    260                      WRITE(numout,*) 
    261                      WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! ' 
    262                      WRITE(numout,*) ' !!!! Something is wrong in the SI3 initialization procedure ' 
    263                      WRITE(numout,*) 
    264                      WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:) 
    265                      WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 
    266                      WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 
    267276                  ENDIF 
    268277                  ! 
    269                ENDIF 
    270                ! 
    271             END DO    
    272          END DO    
    273  
     278               END DO 
     279            END DO 
     280         ENDIF 
     281          
    274282         !--------------------------------------------------------------------- 
    275283         ! 4) Fill in sea ice arrays 

  • NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/EXPREF/namelist_ice_cfg

    r10277 r10399  
    3434&namdyn         !   Ice dynamics 
    3535!------------------------------------------------------------------------------ 
    36    ln_dynFULL       = .false.         !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction) 
     36   ln_dynALL        = .false.         !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction) 
    3737   ln_dynRHGADV     = .false.         !  dyn.: no ridge/raft & no corrections  (rheology + advection) 
    38    ln_dynADV        = .true.          !  dyn.: only advection w prescribed vel.(rn_uvice + advection) 
     38   ln_dynADV1D      = .true.          !  dyn.: only advection 1D                  (Schar & Smolarkiewicz 1996 test case) 
     39   ln_dynADV2D      = .false.         !  dyn.: only advection 2D w prescribed vel.(rn_uvice + advection) 
    3940      rn_uice       =   1.            !        prescribed ice u-velocity 
    4041      rn_vice       =   0.            !        prescribed ice v-velocity 

  • NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEDYN/EXPREF/namelist_ice_cfg

    r9801 r10399  
    3333&namdyn         !   Ice dynamics 
    3434!------------------------------------------------------------------------------ 
    35    ln_dynFULL       = .false.         !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction) 
     35   ln_dynALL        = .false.         !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction) 
    3636   ln_dynRHGADV     = .true.          !  dyn.: no ridge/raft & no corrections  (rheology + advection) 
    37    ln_dynADV        = .false.         !  dyn.: only advection w prescribed vel.(rn_uvice + advection) 
     37   ln_dynADV1D      = .false.         !  dyn.: only advection 1D                  (Schar & Smolarkiewicz 1996 test case) 
     38   ln_dynADV2D      = .false.         !  dyn.: only advection 2D w prescribed vel.(rn_uvice + advection) 
    3839      rn_uice       =   0.5           !        prescribed ice u-velocity 
    3940      rn_vice       =   0.            !        prescribed ice v-velocity 
Note: See TracChangeset for help on using the changeset viewer.