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 10413 for NEMO – NEMO

Changeset 10413 for NEMO


Ignore:
Timestamp:
2018-12-18T18:59:59+01:00 (5 years ago)
Author:
clem
Message:

merge dev_r9947_SI3_advection with the trunk. All sette tests passed. There is probably a conservation issue with the new advection scheme that I should solve but the routine icedyn_adv_umx.F90 is going to change anyway in the next couple of weeks. At worst, the old routine can be plugged without harm

Location:
NEMO/trunk
Files:
27 added
18 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/cfgs/SHARED/field_def_nemo-ice.xml

    r9943 r10413  
    2929          <field id="icemask15"    long_name="Ice mask (0 if ice conc. lower than 15%, 1 otherwise)"   standard_name="sea_ice_mask15"                            unit="" /> 
    3030     <field id="icepres"      long_name="Fraction of time steps with sea ice"                     standard_name="sea_ice_time_fraction"                     unit="" /> 
    31       
     31          <field id="fasticepres"  long_name="Fraction of time steps with landfast ice"                standard_name="fast_ice_time_fraction"                    unit="" /> 
     32  
    3233     <!-- general fields --> 
    3334          <field id="icemass"      long_name="Sea-ice mass per area"                                   standard_name="sea_ice_amount"                            unit="kg/m2"/> 

  • NEMO/trunk/cfgs/SHARED/namelist_ice_ref

    r10075 r10413  
    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) 
    54       rn_uice       =   0.00001       !        prescribed ice u-velocity 
    55       rn_vice       =   0.            !        prescribed ice v-velocity 
    56    rn_ishlat        =   2.            !  free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2) 
    57    ln_landfast      = .false.         !  landfast ice parameterization (T or F)                            
    58       rn_gamma      =   0.15          !     fraction of ocean depth that ice must reach to initiate landfast 
    59                                       !        recommended range: [0.1 ; 0.25] 
    60       rn_icebfr     =  10.            !     maximum bottom stress per unit area of contact [N/m2]                  
    61                                       !        a very large value ensures ice velocity=0 even with a small contact area 
    62                                       !        recommended range: ?? (should be greater than atm-ice stress => >0.1 N/m2) 
    63       rn_lfrelax    =   1.e-5         !     relaxation time scale to reach static friction [s-1] 
     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) 
     55      rn_uice       =   0.5           !        prescribed ice u-velocity 
     56      rn_vice       =   0.5           !        prescribed ice v-velocity 
     57   rn_ishlat        =   2.            !  lbc : free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2) 
     58   ln_landfast_L16  = .false.         !  landfast: parameterization from Lemieux 2016 
     59   ln_landfast_home = .false.         !  landfast: parameterization from "home made" 
     60      rn_depfra     =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast 
     61                                      !          recommended range: [0.1 ; 0.25] - L16=0.125 - home=0.15 
     62      rn_icebfr     =  15.            !        ln_landfast_L16:  maximum bottom stress per unit volume [N/m3] 
     63                                      !        ln_landfast_home: maximum bottom stress per unit area of contact [N/m2] 
     64                                      !          recommended range: ?? L16=15 - home=10 
     65      rn_lfrelax    =   1.e-5         !        relaxation time scale to reach static friction [s-1] 
     66      rn_tensile    =   0.2           !        isotropic tensile strength 
    6467/ 
    6568!------------------------------------------------------------------------------ 
    … …  
    205208   sn_tmi = 'Ice_initialization'    , -12 ,'tmi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    206209   sn_smi = 'Ice_initialization'    , -12 ,'smi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    207    cn_dir ='./' 
     210   cn_dir='./' 
    208211/ 
    209212!------------------------------------------------------------------------------ 

  • NEMO/trunk/cfgs/SPITZ12/EXPREF/namelist_ice_cfg

    r9902 r10413  
    3535&namdyn         !   Ice dynamics 
    3636!------------------------------------------------------------------------------ 
    37    ln_landfast      = .false.         !  landfast ice parameterization (T or F) 
    3837/ 
    3938!------------------------------------------------------------------------------ 

  • NEMO/trunk/src/ICE/ice.F90

    r10068 r10413  
    129129   !                                     !!** ice-dynamics namelist (namdyn) ** 
    130130   REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice 
    131    LOGICAL , PUBLIC ::   ln_landfast      !: landfast ice parameterization (T or F)  
    132    REAL(wp), PUBLIC ::   rn_gamma         !:    fraction of ocean depth that ice must reach to initiate landfast ice 
    133    REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (landfast ice)  
    134    REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction (landfast ice)  
    135    ! 
    136    !                                     !!** ice-rheology namelist (namrhg) ** 
     131   LOGICAL , PUBLIC ::   ln_landfast_L16  !: landfast ice parameterizationfrom lemieux2016  
     132   LOGICAL , PUBLIC ::   ln_landfast_home !: landfast ice parameterizationfrom home made  
     133   REAL(wp), PUBLIC ::   rn_depfra        !:    fraction of ocean depth that ice must reach to initiate landfast ice 
     134   REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (lemieux2016) or per unit volume (home)  
     135   REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction 
     136   REAL(wp), PUBLIC ::   rn_tensile       !:    isotropic tensile strength 
     137   ! 
     138   !                                     !!** ice-ridging/rafting namelist (namdyn_rdgrft) ** 
     139   REAL(wp), PUBLIC ::   rn_crhg          !: determines changes in ice strength (also used for landfast param) 
     140   ! 
     141   !                                     !!** ice-rheology namelist (namdyn_rhg) ** 
    137142   LOGICAL , PUBLIC ::   ln_aEVP          !: using adaptive EVP (T or F)  
    138143   REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9 
    … …  
    140145   INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling 
    141146   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 
    142151   ! 
    143152   !                                     !!** ice-surface forcing namelist (namforcing) ** 

  • NEMO/trunk/src/ICE/icedyn.F90

    r10069 r10413  
    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 
    … …  
    5356   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    5457   !! $Id$ 
    55    !! Software governed by the CeCILL license (see ./LICENSE) 
     58   !! Software governed by the CeCILL licence     (./LICENSE) 
    5659   !!---------------------------------------------------------------------- 
    5760CONTAINS 
    … …  
    7174      INTEGER, INTENT(in) ::   kt     ! ice time step 
    7275      !! 
    73       INTEGER ::   ji, jj, jl         ! dummy loop indices 
    74       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhmax 
     76      INTEGER  ::   ji, jj, jl        ! dummy loop indices 
     77      REAL(wp) ::   zcoefu, zcoefv 
     78      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max 
     79      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i 
    7580      !!-------------------------------------------------------------------- 
    7681      ! 
    77       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 
    7885      ! 
    7986      IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    8289         WRITE(numout,*)'~~~~~~~' 
    8390      ENDIF 
    84  
    8591      !                       
    86       IF( ln_landfast ) THEN            !-- Landfast ice parameterization: define max bottom friction 
     92      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization 
    8793         tau_icebfr(:,:) = 0._wp 
    8894         DO jl = 1, jpl 
    89             WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_gamma )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
    90          END DO 
    91          IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr )   
    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 
     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) ) 
     111            END DO 
     112         END DO 
     113      END DO 
     114      CALL lbc_lnk_multi( zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1. ) 
     115      ! 
     116      ! 
     117      SELECT CASE( nice_dyn )           !-- Set which dynamics is running 
     118 
     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 
     124 
     125      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==! 
     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) 
     131         ALLOCATE( zdivu_i(jpi,jpj) ) 
     132         ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- ! 
     133         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length 
     134         ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s  
     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 
     143         ! --- 
     144         CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice + correction on ice thickness 
     145 
     146         ! diagnostics: divergence at T points  
     147         DO jj = 2, jpjm1 
    97148            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) ) ) 
    100             END DO 
    101          END DO 
    102       END DO 
    103       CALL lbc_lnk( zhmax(:,:,:), 'T', 1. ) 
    104       ! 
    105       ! 
    106       SELECT CASE( nice_dyn )           !-- Set which dynamics is running 
    107  
    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 
    113  
    114       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) 
     149               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     150                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     151            END DO 
     152         END DO 
     153         CALL lbc_lnk( zdivu_i, 'T', 1. ) 
     154         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) ) 
     155 
     156         DEALLOCATE( zdivu_i ) 
     157 
     158      CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities) 
     159         ALLOCATE( zdivu_i(jpi,jpj) ) 
    120160         u_ice(:,:) = rn_uice * umask(:,:,1) 
    121161         v_ice(:,:) = rn_vice * vmask(:,:,1) 
    122          !!CALL RANDOM_NUMBER(u_ice(:,:)) 
    123          !!CALL RANDOM_NUMBER(v_ice(:,:)) 
    124          CALL ice_dyn_adv   ( kt )                            ! -- advection of ice 
     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 ) 
    125178 
    126179      END SELECT 
    127       ! 
    128       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 
    129184      ! 
    130185   END SUBROUTINE ice_dyn 
    131186 
    132187 
    133    SUBROUTINE Hbig( phmax ) 
     188   SUBROUTINE Hbig( phi_max, phs_max ) 
    134189      !!------------------------------------------------------------------- 
    135190      !!                  ***  ROUTINE Hbig  *** 
    136191      !! 
    137192      !! ** Purpose : Thickness correction in case advection scheme creates 
    138       !!              abnormally tick ice 
    139       !! 
    140       !! ** Method  : 1- check whether ice thickness resulting from advection is 
    141       !!                 larger than the surrounding 9-points before advection 
    142       !!                 and reduce it if a) divergence or b) convergence & at_i>0.8 
    143       !!              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 
    144202      !! 
    145203      !! ** input   : Max thickness of the surrounding 9-points 
    146204      !!------------------------------------------------------------------- 
    147       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 
    148206      ! 
    149207      INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    150       REAL(wp) ::   zh, zdv 
     208      REAL(wp) ::   zhi, zhs, zvs_excess, zfra 
    151209      !!------------------------------------------------------------------- 
    152210      ! 
    … …  
    156214         DO jj = 1, jpj 
    157215            DO ji = 1, jpi 
    158                IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to hmax 
     216               IF ( v_i(ji,jj,jl) > 0._wp ) THEN 
    159217                  ! 
    160                   zh  = v_i (ji,jj,jl) / a_i(ji,jj,jl) 
    161                   zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl)   
    162                   ! 
    163                   IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. & 
    164                      & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN 
    165                      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) 
    166226                  ENDIF 
    167227                  ! 
     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                   
    168255               ENDIF 
    169256            END DO 
    … …  
    215302      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
    216303      !! 
    217       NAMELIST/namdyn/ ln_dynFULL, ln_dynRHGADV, ln_dynADV, rn_uice, rn_vice,  & 
    218          &             rn_ishlat  ,                                            & 
    219          &             ln_landfast, rn_gamma , rn_icebfr, rn_lfrelax 
     304      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  & 
     305         &             rn_ishlat ,                                                           & 
     306         &             ln_landfast_L16, ln_landfast_home, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile 
    220307      !!------------------------------------------------------------------- 
    221308      ! 
    … …  
    233320         WRITE(numout,*) '~~~~~~~~~~~~' 
    234321         WRITE(numout,*) '   Namelist namdyn:' 
    235          WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynFULL   = ', ln_dynFULL 
    236          WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV = ', ln_dynRHGADV 
    237          WRITE(numout,*) '      Advection only         (rn_uvice + adv)                 ln_dynADV    = ', ln_dynADV 
    238          WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')' 
    239          WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat    = ', rn_ishlat 
    240          WRITE(numout,*) '      Landfast: param (T or F)                                ln_landfast  = ', ln_landfast 
    241          WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_gamma     = ', rn_gamma 
    242          WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr    = ', rn_icebfr 
    243          WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax   = ', rn_lfrelax 
     322         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynALL       = ', ln_dynALL 
     323         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV    = ', ln_dynRHGADV 
     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 
     326         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)    = (', rn_uice,',', rn_vice,')' 
     327         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat       = ', rn_ishlat 
     328         WRITE(numout,*) '      Landfast: param from Lemieux 2016                       ln_landfast_L16 = ', ln_landfast_L16 
     329         WRITE(numout,*) '      Landfast: param from home made                          ln_landfast_home= ', ln_landfast_home 
     330         WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_depfra       = ', rn_depfra 
     331         WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr       = ', rn_icebfr 
     332         WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax      = ', rn_lfrelax 
     333         WRITE(numout,*) '         isotropic tensile strength                           rn_tensile      = ', rn_tensile 
    244334         WRITE(numout,*) 
    245335      ENDIF 
    … …  
    247337      ioptio = 0  
    248338      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction) 
    249       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 
    250340      !      !--- dynamics without ridging/rafting and corr   (rheology + advection) 
    251341      IF( ln_dynRHGADV ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV    ;   ENDIF 
    252       !      !--- advection only with prescribed ice velocities (from namelist) 
    253       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 
    254346      ! 
    255347      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' ) 
    … …  
    261353      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip' 
    262354      ENDIF 
    263       !                                      !--- NO Landfast ice : set to zero once for all 
    264       IF( .NOT.ln_landfast )   tau_icebfr(:,:) = 0._wp 
     355      !                                      !--- Landfast ice 
     356      IF( .NOT.ln_landfast_L16 .AND. .NOT.ln_landfast_home )   tau_icebfr(:,:) = 0._wp 
     357      ! 
     358      IF ( ln_landfast_L16 .AND. ln_landfast_home ) THEN 
     359         CALL ctl_stop( 'ice_dyn_init: choose one and only one landfast parameterization (ln_landfast_L16 or ln_landfast_home)' ) 
     360      ENDIF 
    265361      ! 
    266362      CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters 

  • NEMO/trunk/src/ICE/icedyn_adv.F90

    r10069 r10413  
    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 
    … …  
    4947   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    5048   !! $Id$ 
    51    !! Software governed by the CeCILL license (see ./LICENSE) 
     49   !! Software governed by the CeCILL licence     (./LICENSE) 
    5250   !!---------------------------------------------------------------------- 
    5351CONTAINS 
    … …  
    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/trunk/src/ICE/icedyn_adv_umx.F90

    r10344 r10413  
    1414   !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders 
    1515   !!   macho             : ??? 
    16    !!   nonosc_2d         : compute monotonic tracer fluxes by a non-oscillatory algorithm  
     16   !!   nonosc            : compute monotonic tracer fluxes by a non-oscillatory algorithm  
    1717   !!---------------------------------------------------------------------- 
    1818   USE phycst         ! physical constant 
    … …  
    2020   USE sbc_oce , ONLY : nn_fsbc   ! update frequency of surface boundary condition 
    2121   USE ice            ! sea-ice variables 
     22   USE icevar         ! sea-ice: operations 
    2223   ! 
    2324   USE in_out_manager ! I/O manager 
    … …  
    3435   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120 
    3536 
     37   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_ai, amaxu, amaxv 
     38    
     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    
     71   ! alternate directions for upstream 
     72   LOGICAL :: ll_upsxy = .TRUE. 
     73 
     74   ! alternate directions for high order 
     75   LOGICAL :: ll_hoxy = .TRUE. 
     76    
     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) 
     80 
     81   ! iterate on the limiter (only nonosc) 
     82   LOGICAL :: ll_limiter_it2 = .FALSE. 
     83    
     84 
    3685   !! * Substitutions 
    3786#  include "vectopt_loop_substitute.h90" 
    … …  
    3988   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    4089   !! $Id$ 
    41    !! Software governed by the CeCILL license (see ./LICENSE) 
     90   !! Software governed by the CeCILL licence     (./LICENSE) 
    4291   !!---------------------------------------------------------------------- 
    4392CONTAINS 
    4493 
    45    SUBROUTINE ice_dyn_adv_umx( k_order, kt, pu_ice, pv_ice,  & 
    46       &                    pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     94   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice,  & 
     95      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    4796      !!---------------------------------------------------------------------- 
    4897      !!                  ***  ROUTINE ice_dyn_adv_umx  *** 
    … …  
    54103      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    55104      !!---------------------------------------------------------------------- 
    56       INTEGER                     , INTENT(in   ) ::   k_order    ! order of the scheme (1-5 or 20) 
     105      INTEGER                     , INTENT(in   ) ::   kn_umx     ! order of the scheme (1-5=UM or 20=CEN2) 
    57106      INTEGER                     , INTENT(in   ) ::   kt         ! time step 
    58107      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity 
    … …  
    70119      ! 
    71120      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
    72       INTEGER  ::   initad                  ! number of sub-timestep for the advection 
    73       REAL(wp) ::   zcfl , zusnit, zdt      !   -      - 
    74       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zudy, zvdx, zcu_box, zcv_box 
     121      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
     122      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
     123      REAL(wp) ::   zcfl , zdt 
     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 
    75127      !!---------------------------------------------------------------------- 
    76128      ! 
    77129      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 
    78130      ! 
    79       ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) ) 
    80131      ! 
    81132      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !         
    … …  
    84135      IF( lk_mpp )   CALL mpp_max( zcfl ) 
    85136 
    86       IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp 
    87       ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp 
    88       ENDIF 
    89  
    90       zdt = rdt_ice / REAL(initad) 
     137      IF( zcfl > 0.5 ) THEN   ;   icycle = 2  
     138      ELSE                    ;   icycle = 1  
     139      ENDIF 
     140       
     141      zdt = rdt_ice / REAL(icycle) 
    91142 
    92143      ! --- transport --- ! 
    … …  
    109160      END DO 
    110161 
     162      IF(.NOT. ALLOCATED(z1_ai))       ALLOCATE(z1_ai(jpi,jpj)) 
     163      IF( ll_zeroup2 ) THEN 
     164         IF(.NOT. ALLOCATED(amaxu))       ALLOCATE(amaxu (jpi,jpj)) 
     165         IF(.NOT. ALLOCATED(amaxv))       ALLOCATE(amaxv (jpi,jpj)) 
     166      ENDIF 
    111167      !---------------! 
    112168      !== advection ==! 
    113169      !---------------! 
    114       DO jt = 1, initad 
    115          CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) )             ! Open water area  
     170      DO jt = 1, icycle 
     171 
     172         IF( ll_ADVopw ) THEN 
     173            ll_dens=.FALSE. 
     174            zamsk = 1._wp 
     175            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) )   ! Open water area  
     176         ELSE 
     177            zai_b(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
     178         ENDIF 
     179          
    116180         DO jl = 1, jpl 
    117             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) )         ! Ice area 
    118             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) )         ! Ice  volume 
    119             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, psv_i(:,:,jl) )        ! Salt content 
    120             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i(:,:,jl) )        ! Age content 
     181            ! 
     182            WHERE( pa_i(:,:,jl) >= epsi20 )   ;   z1_ai_b(:,:) = 1._wp / pa_i(:,:,jl) 
     183            ELSEWHERE                         ;   z1_ai_b(:,:) = 0. 
     184            END WHERE 
     185            ! 
     186            IF( ll_zeroup2 ) THEN 
     187               DO jj = 2, jpjm1 
     188                  DO ji = fs_2, fs_jpim1 
     189                     amaxu(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji,jj-1,jl), pa_i(ji,jj+1,jl), & 
     190                        &                              pa_i(ji+1,jj,jl), pa_i(ji+1,jj-1,jl), pa_i(ji+1,jj+1,jl) ) 
     191                     amaxv(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji-1,jj,jl), pa_i(ji+1,jj,jl), & 
     192                        &                              pa_i(ji,jj+1,jl), pa_i(ji-1,jj+1,jl), pa_i(ji+1,jj+1,jl) ) 
     193                 END DO 
     194               END DO 
     195               CALL lbc_lnk_multi(amaxu, 'T', 1., amaxv, 'T', 1.) 
     196            ENDIF 
     197            ! 
     198            zamsk = 1._wp 
     199            ll_dens=.TRUE. 
     200            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 
     201            ll_dens=.FALSE. 
     202 
     203            WHERE( pa_i(:,:,jl) >= epsi20 )   ;   z1_ai(:,:) = 1._wp / pa_i(:,:,jl) 
     204            ELSEWHERE                         ;   z1_ai(:,:) = 0. 
     205            END WHERE               
     206 
     207            IF( ll_thickness ) THEN 
     208               zua_ho(:,:) = zudy(:,:) 
     209               zva_ho(:,:) = zvdx(:,:) 
     210            ENDIF 
     211             
     212            zamsk = 0._wp ; zhvar(:,:) = pv_i (:,:,jl) * z1_ai_b(:,:) 
     213            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 
     214            IF( ll_thickness )   pv_i(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl) 
     215             
     216            zamsk = 0._wp ; zhvar(:,:) = pv_s (:,:,jl) * z1_ai_b(:,:) 
     217            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 
     218            IF( ll_thickness )   pv_s(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl) 
     219 
     220            zamsk = 0._wp ; zhvar(:,:) = psv_i(:,:,jl) * z1_ai_b(:,:) 
     221            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 
     222 
     223            zamsk = 0._wp ; zhvar(:,:) = poa_i(:,:,jl) * z1_ai_b(:,:) 
     224            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 
     225 
    121226            DO jk = 1, nlay_i 
    122                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) )   ! Ice  heat content 
    123             END DO 
    124             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume 
     227               zamsk = 0._wp ; zhvar(:,:) = pe_i(:,:,jk,jl) * z1_ai_b(:,:) 
     228               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 
     229            END DO 
     230 
    125231            DO jk = 1, nlay_s 
    126                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,jk,jl) )   ! Snow heat content 
    127             END DO 
    128             IF ( ln_pnd_H12 ) THEN 
    129                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction 
    130                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume 
    131             ENDIF 
    132          END DO 
    133       END DO 
    134       ! 
    135       DEALLOCATE( zudy, zvdx, zcu_box, zcv_box ) 
     232               zamsk = 0._wp ; zhvar(:,:) = pe_s(:,:,jk,jl) * z1_ai_b(:,:) 
     233               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 
     234            END DO 
     235            ! 
     236            IF ( ln_pnd_H12 ) THEN   ! melt ponds (must be the last ones to be advected because of z1_ai_b...) 
     237               ! 
     238               WHERE( pa_ip(:,:,jl) >= epsi20 )   ;   z1_ai_b(:,:) = 1._wp / pa_ip(:,:,jl) 
     239               ELSEWHERE                          ;   z1_ai_b(:,:) = 0. 
     240               END WHERE 
     241               ! 
     242               zamsk = 1._wp 
     243               ll_dens=.TRUE. 
     244               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 
     245               ll_dens=.FALSE. 
     246 
     247               WHERE( pa_ip(:,:,jl) >= epsi20 )   ;   z1_ai(:,:) = 1._wp / pa_ip(:,:,jl) 
     248               ELSEWHERE                          ;   z1_ai(:,:) = 0. 
     249               END WHERE               
     250 
     251               zamsk = 0._wp ; zhvar(:,:) = pv_ip(:,:,jl) * z1_ai_b(:,:) 
     252               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 
     253            ENDIF 
     254            ! 
     255            ! 
     256         END DO 
     257 
     258         IF( .NOT. ll_ADVopw ) THEN 
     259            zai_a(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
     260            DO jj = 2, jpjm1 
     261               DO ji = fs_2, fs_jpim1 
     262                  pato_i(ji,jj) = pato_i(ji,jj) - ( zai_a(ji,jj) - zai_b(ji,jj) ) & 
     263                     &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
     264               END DO 
     265            END DO 
     266            CALL lbc_lnk( pato_i(:,:), 'T',  1. ) 
     267         ENDIF 
     268          
     269      END DO 
    136270      ! 
    137271   END SUBROUTINE ice_dyn_adv_umx 
    138272 
    139273    
    140    SUBROUTINE adv_umx( k_order, kt, pdt, puc, pvc, pubox, pvbox, ptc ) 
     274   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho ) 
    141275      !!---------------------------------------------------------------------- 
    142276      !!                  ***  ROUTINE adv_umx  *** 
    … …  
    151285      !! ** Action : - pt  the after advective tracer 
    152286      !!---------------------------------------------------------------------- 
    153       INTEGER                     , INTENT(in   ) ::   k_order        ! order of the ULTIMATE scheme 
    154       INTEGER                     , INTENT(in   ) ::   kt             ! number of iteration 
    155       REAL(wp)                    , INTENT(in   ) ::   pdt            ! tracer time-step 
    156       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc  , pvc     ! 2 ice velocity components => u*e2 
    157       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox   ! upstream velocity 
    158       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptc            ! tracer content field 
     287      REAL(wp)                    , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0) 
     288      INTEGER                     , INTENT(in   )           ::   kn_umx         ! order of the scheme (1-5=UM or 20=CEN2) 
     289      INTEGER                     , INTENT(in   )           ::   jt             ! number of sub-iteration 
     290      INTEGER                     , INTENT(in   )           ::   kt             ! number of iteration 
     291      REAL(wp)                    , INTENT(in   )           ::   pdt            ! tracer time-step 
     292      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2 
     293      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u 
     294      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity 
     295      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   pt             ! tracer field 
     296      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptc            ! tracer content field 
     297      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes 
    159298      ! 
    160299      INTEGER  ::   ji, jj           ! dummy loop indices   
    161300      REAL(wp) ::   ztra             ! local scalar 
    162       REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfu_ho, zt_u, zt_ups 
    163       REAL(wp), DIMENSION(jpi,jpj) ::   zfv_ups, zfv_ho, zt_v, ztrd 
    164       !!---------------------------------------------------------------------- 
    165       ! 
    166       !  upstream advection with initial mass fluxes & intermediate update 
    167       ! -------------------------------------------------------------------- 
    168       DO jj = 1, jpjm1         ! upstream tracer flux in the i and j direction 
    169          DO ji = 1, fs_jpim1   ! vector opt. 
    170             zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * ptc(ji,jj) + MIN( puc(ji,jj), 0._wp ) * ptc(ji+1,jj) 
    171             zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * ptc(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * ptc(ji,jj+1) 
    172          END DO 
    173       END DO 
     301      INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid) 
     302      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt 
     303      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups, zt_ups   ! only for nonosc  
     304      !!---------------------------------------------------------------------- 
     305      ! 
     306      !  upstream (_ups) advection with initial mass fluxes 
     307      ! --------------------------------------------------- 
     308      IF( ll_clem )    zfu_ups=0.; zfv_ups=0. 
     309 
     310      IF( ll_gurvan .AND. pamsk==0. ) THEN 
     311         DO jj = 2, jpjm1 
     312            DO ji = fs_2, fs_jpim1 
     313               pt(ji,jj) = ( pt (ji,jj) + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) & 
     314                  &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     315            END DO 
     316         END DO 
     317         CALL lbc_lnk( pt, 'T', 1. ) 
     318      ENDIF 
     319 
    174320       
    175       DO jj = 2, jpjm1            ! total intermediate advective trends 
    176          DO ji = fs_2, fs_jpim1   ! vector opt. 
    177             ztra = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   & 
    178                &       + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1)   ) * r1_e1e2t(ji,jj) 
    179             ! 
    180             ztrd(ji,jj) =                         ztra                         ! upstream trend [ -div(uh) or -div(uhT) ]   
    181             zt_ups (ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)      ! guess after content field with monotonic scheme 
    182          END DO 
    183       END DO 
    184       CALL lbc_lnk( zt_ups, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign) 
    185        
     321      IF( .NOT. ll_upsxy ) THEN 
     322 
     323         ! fluxes 
     324         DO jj = 1, jpjm1 
     325            DO ji = 1, fs_jpim1 
     326               IF( ll_clem ) THEN 
     327                  zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj) 
     328                  zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1) 
     329               ELSE 
     330                  zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj) 
     331                  zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1) 
     332               ENDIF 
     333            END DO 
     334         END DO 
     335 
     336      ELSE 
     337         ! 1 if advection of A 
     338         ! z1_ai already defined IF advection of other tracers 
     339         IF( pamsk == 1. )   z1_ai(:,:) = 1._wp 
     340         ! 
     341         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     342            ! flux in x-direction 
     343            DO jj = 1, jpjm1 
     344               DO ji = 1, fs_jpim1 
     345                  IF( ll_clem ) THEN 
     346                     zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj) 
     347                  ELSE 
     348                     zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj) 
     349                  ENDIF 
     350               END DO 
     351            END DO 
     352             
     353            ! first guess of tracer content from u-flux 
     354            DO jj = 2, jpjm1 
     355               DO ji = fs_2, fs_jpim1   ! vector opt. 
     356                  IF( ll_clem ) THEN 
     357                     IF( ll_gurvan ) THEN 
     358                        zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     359                     ELSE 
     360                        zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     361                           &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
     362                           &         * tmask(ji,jj,1) 
     363                     ENDIF 
     364                  ELSE 
     365                     zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) & 
     366                        &         * tmask(ji,jj,1) 
     367                  ENDIF 
     368!!                  IF( ji==26 .AND. jj==86) THEN 
     369!!                     WRITE(numout,*) '************************' 
     370!!                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 
     371!!                  ENDIF 
     372               END DO 
     373            END DO 
     374            CALL lbc_lnk( zpt, 'T', 1. ) 
     375            ! 
     376            ! flux in y-direction 
     377            DO jj = 1, jpjm1 
     378               DO ji = 1, fs_jpim1 
     379                  IF( ll_clem ) THEN 
     380                     zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1) 
     381                  ELSE 
     382                     zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1) 
     383                  ENDIF 
     384               END DO 
     385            END DO 
     386 
     387         ! 
     388         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     389            ! flux in y-direction 
     390            DO jj = 1, jpjm1 
     391               DO ji = 1, fs_jpim1 
     392                  IF( ll_clem ) THEN 
     393                     zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1) 
     394                  ELSE 
     395                     zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1) 
     396                  ENDIF 
     397               END DO 
     398            END DO 
     399 
     400            ! first guess of tracer content from v-flux 
     401            DO jj = 2, jpjm1 
     402               DO ji = fs_2, fs_jpim1   ! vector opt. 
     403                  IF( ll_clem ) THEN 
     404                     IF( ll_gurvan ) THEN 
     405                        zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     406                     ELSE 
     407                        zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 
     408                        &                        + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
     409                        &         * tmask(ji,jj,1) 
     410                     ENDIF 
     411                  ELSE 
     412                     zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) & 
     413                        &         * tmask(ji,jj,1) 
     414                  ENDIF 
     415!!                  IF( ji==26 .AND. jj==86) THEN 
     416!!                     WRITE(numout,*) '************************' 
     417!!                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 
     418!!                  ENDIF 
     419                END DO 
     420            END DO 
     421            CALL lbc_lnk( zpt, 'T', 1. ) 
     422            ! 
     423            ! flux in x-direction 
     424            DO jj = 1, jpjm1 
     425               DO ji = 1, fs_jpim1 
     426                  IF( ll_clem ) THEN 
     427                     zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj) 
     428                  ELSE 
     429                     zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj) 
     430                  ENDIF 
     431               END DO 
     432            END DO 
     433            ! 
     434         ENDIF 
     435          
     436      ENDIF 
     437 
     438      IF( ll_clem .AND. kn_limiter /= 1 )  & 
     439         & CALL ctl_stop( 'STOP', 'icedyn_adv_umx: ll_clem incompatible with limiters other than nonosc' ) 
     440 
     441      IF( ll_zeroup2 ) THEN 
     442         DO jj = 1, jpjm1 
     443            DO ji = 1, fs_jpim1   ! vector opt. 
     444               IF( amaxu(ji,jj) == 0._wp )   zfu_ups(ji,jj) = 0._wp 
     445               IF( amaxv(ji,jj) == 0._wp )   zfv_ups(ji,jj) = 0._wp 
     446            END DO 
     447         END DO 
     448      ENDIF 
     449 
     450      ! guess after content field with upstream scheme 
     451      DO jj = 2, jpjm1 
     452         DO ji = fs_2, fs_jpim1 
     453            ztra          = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   & 
     454               &                + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1) ) * r1_e1e2t(ji,jj) 
     455            IF( ll_clem ) THEN 
     456               IF( ll_gurvan ) THEN 
     457                  zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
     458               ELSE 
     459                  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) & 
     460                     &                                      + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     461               ENDIF 
     462            ELSE 
     463               zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
     464            ENDIF 
     465!!            IF( ji==26 .AND. jj==86) THEN 
     466!!               WRITE(numout,*) '**************************' 
     467!!               WRITE(numout,*) 'zt upstream',zt_ups(ji,jj) 
     468!!            ENDIF 
     469         END DO 
     470      END DO 
     471      CALL lbc_lnk( zt_ups, 'T', 1. ) 
     472 
    186473      ! High order (_ho) fluxes  
    187474      ! ----------------------- 
    188       SELECT CASE( k_order ) 
    189       CASE ( 20 )                          ! centered second order 
     475      SELECT CASE( kn_umx ) 
     476         ! 
     477      CASE ( 20 )                          !== centered second order ==! 
     478         ! 
     479         CALL cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  & 
     480            &       zt_ups, zfu_ups, zfv_ups ) 
     481         ! 
     482      CASE ( 1:5 )                         !== 1st to 5th order ULTIMATE-MACHO scheme ==! 
     483         ! 
     484         CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho,  & 
     485            &        zt_ups, zfu_ups, zfv_ups ) 
     486         ! 
     487      END SELECT 
     488 
     489      IF( ll_thickness ) THEN 
     490         ! final trend with corrected fluxes 
     491         ! ------------------------------------ 
     492         DO jj = 2, jpjm1 
     493            DO ji = fs_2, fs_jpim1 
     494               IF( ll_gurvan ) THEN 
     495                  ztra       = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj)  
     496               ELSE 
     497                  ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  &  
     498                     &           + ( pt(ji,jj) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) & 
     499                     &           + ( pt(ji,jj) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj) 
     500               ENDIF 
     501               pt(ji,jj) = ( pt(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
     502 
     503               IF( pt(ji,jj) < -epsi20 ) THEN 
     504                  WRITE(numout,*) 'T<0 ',pt(ji,jj) 
     505               ENDIF 
     506                
     507               IF( pt(ji,jj) < 0._wp .AND. pt(ji,jj) >= -epsi20 )   pt(ji,jj) = 0._wp 
     508                
     509!!               IF( ji==26 .AND. jj==86) THEN 
     510!!                  WRITE(numout,*) 'zt high order',pt(ji,jj) 
     511!!               ENDIF 
     512            END DO 
     513         END DO 
     514         CALL lbc_lnk( pt, 'T',  1. ) 
     515      ENDIF 
     516    
     517      ! Rachid trick 
     518      ! ------------ 
     519     IF( ll_clem ) THEN 
     520         IF( pamsk == 0. ) THEN 
     521            DO jj = 1, jpjm1 
     522               DO ji = 1, fs_jpim1 
     523                  IF( ABS( puc(ji,jj) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN 
     524                     zfu_ho (ji,jj) = zfu_ho (ji,jj) * puc(ji,jj) / pu(ji,jj) 
     525                     zfu_ups(ji,jj) = zfu_ups(ji,jj) * puc(ji,jj) / pu(ji,jj) 
     526                  ELSE 
     527                     zfu_ho (ji,jj) = 0._wp 
     528                     zfu_ups(ji,jj) = 0._wp 
     529                  ENDIF 
     530                  ! 
     531                  IF( ABS( pvc(ji,jj) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN 
     532                     zfv_ho (ji,jj) = zfv_ho (ji,jj) * pvc(ji,jj) / pv(ji,jj) 
     533                     zfv_ups(ji,jj) = zfv_ups(ji,jj) * pvc(ji,jj) / pv(ji,jj) 
     534                  ELSE 
     535                     zfv_ho (ji,jj) = 0._wp   
     536                     zfv_ups(ji,jj) = 0._wp   
     537                  ENDIF 
     538               ENDDO 
     539            ENDDO 
     540         ENDIF 
     541      ENDIF 
     542 
     543      IF( ll_zeroup5 ) THEN 
     544         DO jj = 2, jpjm1 
     545            DO ji = 2, fs_jpim1   ! vector opt. 
     546               zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     547                  &                      - ( zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  ) * tmask(ji,jj,1) 
     548               IF( zpt(ji,jj) < 0. ) THEN 
     549                  zfu_ho(ji,jj)   = zfu_ups(ji,jj) 
     550                  zfu_ho(ji-1,jj) = zfu_ups(ji-1,jj) 
     551                  zfv_ho(ji,jj)   = zfv_ups(ji,jj) 
     552                  zfv_ho(ji,jj-1) = zfv_ups(ji,jj-1) 
     553               ENDIF 
     554            END DO 
     555         END DO 
     556         CALL lbc_lnk_multi( zfu_ho, 'U',  -1., zfv_ho, 'V',  -1. ) 
     557      ENDIF 
     558 
     559      ! output upstream trend of concentration and high order fluxes 
     560      ! ------------------------------------------------------------ 
     561      IF( ll_dens ) THEN 
     562         ! high order u*a 
    190563         DO jj = 1, jpjm1 
    191             DO ji = 1, fs_jpim1   ! vector opt. 
    192                zfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( ptc(ji,jj) + ptc(ji+1,jj) ) 
    193                zfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj) + ptc(ji,jj+1) ) 
    194             END DO 
    195          END DO 
    196          ! 
    197       CASE ( 1:5 )                      ! 1st to 5th order ULTIMATE-MACHO scheme 
    198          CALL macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v ) 
    199          ! 
     564            DO ji = 1, fs_jpim1 
     565               pua_ho (ji,jj) = zfu_ho (ji,jj) 
     566               pva_ho (ji,jj) = zfv_ho (ji,jj) 
     567            END DO 
     568         END DO 
     569         !!CALL lbc_lnk( pua_ho, 'U',  -1. ) ! clem: not needed I think 
     570         !!CALL lbc_lnk( pva_ho, 'V',  -1. ) 
     571      ENDIF 
     572 
     573 
     574      IF( .NOT.ll_thickness ) THEN 
     575         ! final trend with corrected fluxes 
     576         ! ------------------------------------ 
    200577         DO jj = 2, jpjm1 
    201             DO ji = 1, fs_jpim1   ! vector opt. 
    202                zfu_ho(ji,jj) = puc(ji,jj) * zt_u(ji,jj) 
    203             END DO 
    204          END DO 
     578            DO ji = fs_2, fs_jpim1  
     579               ztra = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  &  !        Div(uaH) or        Div(ua) 
     580                  &     * r1_e1e2t(ji,jj) * pdt   
     581 
     582               !!IF( ptc(ji,jj)+ztra < 0._wp ) THEN 
     583               !!   ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) )  &  !        Div(uaH) or        Div(ua) 
     584               !!      &     * r1_e1e2t(ji,jj) * pdt                  
     585               !!ENDIF 
     586               !!IF( ptc(ji,jj)+ztra < 0._wp ) THEN 
     587               !!   WRITE(numout,*) 'Tc<0 ',ptc(ji,jj)+ztra 
     588               !!   ztra = 0._wp 
     589               !!ENDIF 
     590 
     591               ptc(ji,jj) = ( ptc(ji,jj) + ztra ) * tmask(ji,jj,1) 
     592                              
     593!!               IF( ji==26 .AND. jj==86) THEN 
     594!!                  WRITE(numout,*) 'ztc high order',ptc(ji,jj) 
     595!!               ENDIF 
     596                
     597            END DO 
     598         END DO 
     599         CALL lbc_lnk( ptc, 'T',  1. ) 
     600      ENDIF 
     601       
     602      ! 
     603   END SUBROUTINE adv_umx 
     604 
     605   SUBROUTINE cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, & 
     606      &             pt_ups, pfu_ups, pfv_ups ) 
     607      !!--------------------------------------------------------------------- 
     608      !!                    ***  ROUTINE macho  *** 
     609      !!      
     610      !! **  Purpose :   compute   
     611      !! 
     612      !! **  Method  :   ... ??? 
     613      !!                 TIM = transient interpolation Modeling  
     614      !! 
     615      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
     616      !!---------------------------------------------------------------------- 
     617      REAL(wp)                    , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     618      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter 
     619      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration 
     620      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration 
     621      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step 
     622      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt               ! tracer fields 
     623      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components 
     624      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components 
     625      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc              ! tracer content at before time step  
     626      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes  
     627      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content  
     628      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes  
     629      ! 
     630      INTEGER  ::   ji, jj    ! dummy loop indices 
     631      LOGICAL  ::   ll_xy = .TRUE. 
     632      REAL(wp), DIMENSION(jpi,jpj) ::   zzt 
     633      !!---------------------------------------------------------------------- 
     634      ! 
     635      IF( .NOT.ll_xy ) THEN   !-- no alternate directions --! 
     636         ! 
    205637         DO jj = 1, jpjm1 
    206             DO ji = fs_2, fs_jpim1   ! vector opt. 
    207                zfv_ho(ji,jj) = pvc(ji,jj) * zt_v(ji,jj) 
    208             END DO 
    209          END DO 
    210          ! 
    211       END SELECT 
     638            DO ji = 1, fs_jpim1 
     639               pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) ) 
     640               pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) ) 
     641            END DO 
     642         END DO 
     643         IF    ( kn_limiter == 1 ) THEN 
     644            IF( ll_clem ) THEN 
     645               CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     646            ELSE 
     647               CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     648            ENDIF 
     649         ELSEIF( kn_limiter == 2 ) THEN 
     650            CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     651            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     652         ELSEIF( kn_limiter == 3 ) THEN 
     653            CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     654            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     655         ENDIF 
     656         ! 
     657      ELSE                    !-- alternate directions --! 
     658         ! 
     659         IF( pamsk == 1. )   z1_ai(:,:) = 1._wp 
     660         ! 
     661         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     662            ! 
     663            ! flux in x-direction 
     664            DO jj = 1, jpjm1 
     665               DO ji = 1, fs_jpim1 
     666                  pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) ) 
     667               END DO 
     668            END DO 
     669            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     670            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     671 
     672            ! first guess of tracer content from u-flux 
     673            DO jj = 2, jpjm1 
     674               DO ji = fs_2, fs_jpim1   ! vector opt. 
     675                  IF( ll_clem ) THEN 
     676                     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) 
     677                  ELSE                      
     678                     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) 
     679                  ENDIF 
     680               END DO 
     681            END DO 
     682            CALL lbc_lnk( zzt, 'T', 1. ) 
     683 
     684            ! flux in y-direction 
     685            DO jj = 1, jpjm1 
     686               DO ji = 1, fs_jpim1 
     687                  pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) ) 
     688               END DO 
     689            END DO 
     690            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     691            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     692 
     693         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     694            ! 
     695            ! flux in y-direction 
     696            DO jj = 1, jpjm1 
     697               DO ji = 1, fs_jpim1 
     698                  pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) ) 
     699               END DO 
     700            END DO 
     701            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     702            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     703            ! 
     704            ! first guess of tracer content from v-flux 
     705            DO jj = 2, jpjm1 
     706               DO ji = fs_2, fs_jpim1   ! vector opt. 
     707                  IF( ll_clem ) THEN 
     708                     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) 
     709                  ELSE 
     710                     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) 
     711                  ENDIF 
     712               END DO 
     713            END DO 
     714            CALL lbc_lnk( zzt, 'T', 1. ) 
     715            ! 
     716            ! flux in x-direction 
     717            DO jj = 1, jpjm1 
     718               DO ji = 1, fs_jpim1 
     719                  pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) ) 
     720               END DO 
     721            END DO 
     722            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     723            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     724 
     725         ENDIF 
     726         IF( ll_clem ) THEN 
     727            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 ) 
     728         ELSE 
     729            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 ) 
     730         ENDIF 
    212731          
    213       ! antidiffusive flux : high order minus low order 
    214       ! -------------------------------------------------- 
    215       DO jj = 2, jpjm1 
    216          DO ji = 1, fs_jpim1   ! vector opt. 
    217             zfu_ho(ji,jj) = zfu_ho(ji,jj) - zfu_ups(ji,jj) 
    218          END DO 
    219       END DO 
    220       DO jj = 1, jpjm1 
    221          DO ji = fs_2, fs_jpim1   ! vector opt. 
    222             zfv_ho(ji,jj) = zfv_ho(ji,jj) - zfv_ups(ji,jj) 
    223          END DO 
    224       END DO 
    225        
    226       ! monotonicity algorithm 
    227       ! ------------------------- 
    228       CALL nonosc_2d( ptc, zfu_ho, zfv_ho, zt_ups, pdt ) 
    229        
    230       ! final trend with corrected fluxes 
    231       ! ------------------------------------ 
    232       DO jj = 2, jpjm1 
    233          DO ji = fs_2, fs_jpim1   ! vector opt.   
    234             ztra       = ztrd(ji,jj)  - (  zfu_ho(ji,jj) - zfu_ho(ji-1,jj  )   & 
    235                &                         + zfv_ho(ji,jj) - zfv_ho(ji  ,jj-1) ) * r1_e1e2t(ji,jj)   
    236             ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
    237          END DO 
    238       END DO 
    239       CALL lbc_lnk( ptc, 'T',  1. ) 
    240       ! 
    241    END SUBROUTINE adv_umx 
    242  
    243  
    244    SUBROUTINE macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 
     732      ENDIF 
     733    
     734   END SUBROUTINE cen2 
     735 
     736    
     737   SUBROUTINE macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, pt_u, pt_v, pfu_ho, pfv_ho, & 
     738      &              pt_ups, pfu_ups, pfv_ups ) 
     739      !!--------------------------------------------------------------------- 
     740      !!                    ***  ROUTINE macho  *** 
     741      !!      
     742      !! **  Purpose :   compute   
     743      !! 
     744      !! **  Method  :   ... ??? 
     745      !!                 TIM = transient interpolation Modeling  
     746      !! 
     747      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
     748      !!---------------------------------------------------------------------- 
     749      REAL(wp)                    , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     750      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter 
     751      INTEGER                     , INTENT(in   ) ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2) 
     752      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration 
     753      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration 
     754      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step 
     755      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt               ! tracer fields 
     756      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components 
     757      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components 
     758      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox     ! upstream velocity 
     759      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc              ! tracer content at before time step  
     760      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u, pt_v       ! tracer at u- and v-points  
     761      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes  
     762      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content  
     763      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes  
     764      ! 
     765      INTEGER  ::   ji, jj    ! dummy loop indices 
     766      REAL(wp) ::   ztra 
     767      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zzfu_ho, zzfv_ho 
     768      !!---------------------------------------------------------------------- 
     769      ! 
     770      IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     771         ! 
     772         !                                                        !--  ultimate interpolation of pt at u-point  --! 
     773         CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 
     774         !                                                        !--  limiter in x --! 
     775         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     776         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     777         !                                                        !--  advective form update in zzt  --! 
     778 
     779         IF( ll_1stguess_clem ) THEN 
     780 
     781            ! first guess of tracer content from u-flux 
     782            DO jj = 2, jpjm1 
     783               DO ji = fs_2, fs_jpim1   ! vector opt. 
     784                  IF( ll_clem ) THEN 
     785                     IF( ll_gurvan ) THEN 
     786                        zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     787                     ELSE 
     788                        zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     789                           &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     790                     ENDIF 
     791                  ELSE 
     792                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     793                  ENDIF 
     794               END DO 
     795            END DO 
     796            CALL lbc_lnk( zzt, 'T', 1. ) 
     797 
     798         ELSE 
     799 
     800            DO jj = 2, jpjm1 
     801               DO ji = fs_2, fs_jpim1   ! vector opt. 
     802                  IF( ll_gurvan ) THEN 
     803                     zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  & 
     804                        &                   - pt   (ji,jj) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) 
     805                  ELSE 
     806                     zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  & 
     807                        &                   - pt   (ji,jj) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk 
     808                  ENDIF 
     809                  zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 
     810               END DO 
     811            END DO 
     812            CALL lbc_lnk( zzt, 'T', 1. ) 
     813         ENDIF 
     814         ! 
     815         !                                                        !--  ultimate interpolation of pt at v-point  --! 
     816         IF( ll_hoxy ) THEN 
     817            CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho ) 
     818         ELSE 
     819            CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 
     820         ENDIF 
     821         !                                                        !--  limiter in y --! 
     822         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     823         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     824         !          
     825         ! 
     826      ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     827         ! 
     828         !                                                        !--  ultimate interpolation of pt at v-point  --! 
     829         CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 
     830         !                                                        !--  limiter in y --! 
     831         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     832         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     833         !                                                        !--  advective form update in zzt  --! 
     834         IF( ll_1stguess_clem ) THEN 
     835             
     836            ! first guess of tracer content from v-flux  
     837            DO jj = 2, jpjm1 
     838               DO ji = fs_2, fs_jpim1   ! vector opt. 
     839                  IF( ll_clem ) THEN 
     840                     IF( ll_gurvan ) THEN 
     841                        zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     842                     ELSE 
     843                        zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 
     844                           &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     845                     ENDIF 
     846                  ELSE 
     847                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) & 
     848                        &         * tmask(ji,jj,1) 
     849                  ENDIF 
     850                END DO 
     851            END DO 
     852            CALL lbc_lnk( zzt, 'T', 1. ) 
     853             
     854         ELSE 
     855             
     856            DO jj = 2, jpjm1 
     857               DO ji = fs_2, fs_jpim1 
     858                  IF( ll_gurvan ) THEN 
     859                     zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  & 
     860                        &                   - pt   (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) 
     861                  ELSE 
     862                     zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  & 
     863                        &                   - pt   (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk 
     864                  ENDIF 
     865                  zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 
     866               END DO 
     867            END DO 
     868            CALL lbc_lnk( zzt, 'T', 1. ) 
     869         ENDIF 
     870         ! 
     871         !                                                        !--  ultimate interpolation of pt at u-point  --! 
     872         IF( ll_hoxy ) THEN 
     873            CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho ) 
     874         ELSE 
     875            CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 
     876         ENDIF 
     877         !                                                        !--  limiter in x --! 
     878         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     879         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     880         ! 
     881         ! 
     882      ENDIF 
     883 
     884      
     885      IF( kn_limiter == 1 ) THEN 
     886         IF( .NOT. ll_limiter_it2 ) THEN 
     887            IF( ll_clem ) THEN 
     888               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     889            ELSE 
     890               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     891            ENDIF 
     892         ELSE 
     893            zzfu_ho(:,:) = pfu_ho(:,:) 
     894            zzfv_ho(:,:) = pfv_ho(:,:) 
     895            ! 1st iteration of nonosc (limit the flux with the upstream solution) 
     896            IF( ll_clem ) THEN 
     897               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 
     898            ELSE 
     899               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 
     900            ENDIF 
     901            ! guess after content field with high order 
     902            DO jj = 2, jpjm1 
     903               DO ji = fs_2, fs_jpim1 
     904                  ztra       = - ( zzfu_ho(ji,jj) - zzfu_ho(ji-1,jj) + zzfv_ho(ji,jj) - zzfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     905                  zzt(ji,jj) =   ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
     906               END DO 
     907            END DO 
     908            CALL lbc_lnk( zzt, 'T', 1. ) 
     909            ! 2nd iteration of nonosc (limit the flux with the limited high order solution) 
     910            IF( ll_clem ) THEN 
     911               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 
     912            ELSE 
     913               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 
     914            ENDIF 
     915         ENDIF 
     916      ENDIF 
     917      ! 
     918   END SUBROUTINE macho 
     919 
     920 
     921   SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 
    245922      !!--------------------------------------------------------------------- 
    246923      !!                    ***  ROUTINE ultimate_x  *** 
    … …  
    253930      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    254931      !!---------------------------------------------------------------------- 
    255       INTEGER                     , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme 
    256       INTEGER                     , INTENT(in   ) ::   kt         ! number of iteration 
    257       REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step 
    258       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc        ! tracer fields 
    259       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc   ! 2 ice velocity components 
    260       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox   ! upstream velocity 
    261       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u, pt_v ! tracer at u- and v-points  
    262       ! 
    263       INTEGER  ::   ji, jj    ! dummy loop indices 
    264       REAL(wp) ::   zc_box    !   -      - 
    265       REAL(wp), DIMENSION(jpi,jpj) :: zzt 
    266       !!---------------------------------------------------------------------- 
    267       ! 
    268       IF( MOD( (kt - 1) / nn_fsbc , 2 ) == 0 ) THEN         !==  odd ice time step:  adv_x then adv_y  ==! 
    269          ! 
    270          !                                                           !--  ultimate interpolation of pt at u-point  --! 
    271          CALL ultimate_x( k_order, pdt, ptc, puc, pt_u ) 
    272          ! 
    273          !                                                           !--  advective form update in zzt  --! 
    274          DO jj = 2, jpjm1 
    275             DO ji = fs_2, fs_jpim1   ! vector opt. 
    276                zzt(ji,jj) = ptc(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  & 
    277                   &                    - ptc  (ji,jj) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj) 
    278                zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 
    279             END DO 
    280          END DO 
    281          CALL lbc_lnk( zzt, 'T', 1. ) 
    282          ! 
    283          !                                                           !--  ultimate interpolation of pt at v-point  --! 
    284          CALL ultimate_y( k_order, pdt, zzt, pvc, pt_v ) 
    285          ! 
    286       ELSE                                                  !==  even ice time step:  adv_y then adv_x  ==! 
    287          ! 
    288          !                                                           !--  ultimate interpolation of pt at v-point  --! 
    289          CALL ultimate_y( k_order, pdt, ptc, pvc, pt_v ) 
    290          ! 
    291          !                                                           !--  advective form update in zzt  --! 
    292          DO jj = 2, jpjm1 
    293             DO ji = fs_2, fs_jpim1 
    294                zzt(ji,jj) = ptc(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  & 
    295                   &                    - ptc  (ji,jj) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj) 
    296                zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 
    297             END DO 
    298          END DO 
    299          CALL lbc_lnk( zzt, 'T', 1. ) 
    300          ! 
    301          !                                                           !--  ultimate interpolation of pt at u-point  --! 
    302          CALL ultimate_x( k_order, pdt, zzt, puc, pt_u ) 
    303          !       
    304       ENDIF       
    305       ! 
    306    END SUBROUTINE macho 
    307  
    308  
    309    SUBROUTINE ultimate_x( k_order, pdt, pt, puc, pt_u ) 
    310       !!--------------------------------------------------------------------- 
    311       !!                    ***  ROUTINE ultimate_x  *** 
    312       !!      
    313       !! **  Purpose :   compute   
    314       !! 
    315       !! **  Method  :   ... ??? 
    316       !!                 TIM = transient interpolation Modeling  
    317       !! 
    318       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    319       !!---------------------------------------------------------------------- 
    320       INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index 
     932      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
    321933      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step 
    322       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc       ! ice i-velocity component 
     934      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu        ! ice i-velocity component 
     935      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc       ! ice i-velocity * A component 
    323936      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields 
    324937      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u      ! tracer at u-point  
    325       ! 
    326       INTEGER  ::   ji, jj       ! dummy loop indices 
     938      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho    ! high order flux  
     939      ! 
     940      INTEGER  ::   ji, jj             ! dummy loop indices 
    327941      REAL(wp) ::   zcu, zdx2, zdx4    !   -      - 
    328       REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu2, ztu3, ztu4 
     942      REAL(wp), DIMENSION(jpi,jpj) ::   ztu1, ztu2, ztu3, ztu4 
    329943      !!---------------------------------------------------------------------- 
    330944      ! 
    … …  
    354968      ! 
    355969      ! 
    356       SELECT CASE (k_order ) 
     970      SELECT CASE (kn_umx ) 
    357971      ! 
    358972      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21) 
    359973         !         
    360          DO jj = 2, jpjm1 
     974         DO jj = 1, jpjm1 
    361975            DO ji = 1, fs_jpim1   ! vector opt. 
    362                pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj) + pt(ji,jj)   & 
    363                   &                                    - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
     976               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                              pt(ji+1,jj) + pt(ji,jj)   & 
     977                  &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
    364978            END DO 
    365979         END DO 
    … …  
    367981      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23) 
    368982         ! 
    369          DO jj = 2, jpjm1 
     983         DO jj = 1, jpjm1 
    370984            DO ji = 1, fs_jpim1   ! vector opt. 
    371                zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     985               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    372986               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                                   pt(ji+1,jj) + pt(ji,jj)   & 
    373987                  &                                               -              zcu   * ( pt(ji+1,jj) - pt(ji,jj) ) )  
    … …  
    377991      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24) 
    378992         ! 
    379          DO jj = 2, jpjm1 
     993         DO jj = 1, jpjm1 
    380994            DO ji = 1, fs_jpim1   ! vector opt. 
    381                zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     995               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    382996               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    383997!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    … …  
    3911005      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27) 
    3921006         ! 
    393          DO jj = 2, jpjm1 
     1007         DO jj = 1, jpjm1 
    3941008            DO ji = 1, fs_jpim1   ! vector opt. 
    395                zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     1009               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    3961010               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    3971011!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    … …  
    4051019      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29) 
    4061020         ! 
    407          DO jj = 2, jpjm1 
     1021         DO jj = 1, jpjm1 
    4081022            DO ji = 1, fs_jpim1   ! vector opt. 
    409                zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     1023               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    4101024               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    4111025!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    … …  
    4211035         ! 
    4221036      END SELECT 
     1037      !                                                     !-- High order flux in i-direction  --! 
     1038      IF( ll_neg ) THEN 
     1039         DO jj = 1, jpjm1 
     1040            DO ji = 1, fs_jpim1 
     1041               IF( pt_u(ji,jj) < 0._wp ) THEN 
     1042                  pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                              pt(ji+1,jj) + pt(ji,jj)   & 
     1043                     &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
     1044               ENDIF 
     1045            END DO 
     1046         END DO 
     1047      ENDIF 
     1048 
     1049      DO jj = 1, jpjm1 
     1050         DO ji = 1, fs_jpim1   ! vector opt. 
     1051            IF( ll_clem ) THEN 
     1052               pfu_ho(ji,jj) = pu(ji,jj) * pt_u(ji,jj) 
     1053            ELSE 
     1054               pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj) 
     1055            ENDIF 
     1056         END DO 
     1057      END DO 
    4231058      ! 
    4241059   END SUBROUTINE ultimate_x 
    4251060    
    4261061  
    427    SUBROUTINE ultimate_y( k_order, pdt, pt, pvc, pt_v ) 
     1062   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 
    4281063      !!--------------------------------------------------------------------- 
    4291064      !!                    ***  ROUTINE ultimate_y  *** 
    … …  
    4361071      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    4371072      !!---------------------------------------------------------------------- 
    438       INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index 
     1073      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
    4391074      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step 
    440       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvc       ! ice j-velocity component 
     1075      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pv        ! ice j-velocity component 
     1076      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvc       ! ice j-velocity*A component 
    4411077      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields 
    4421078      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_v      ! tracer at v-point  
     1079      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfv_ho    ! high order flux  
    4431080      ! 
    4441081      INTEGER  ::   ji, jj       ! dummy loop indices 
    4451082      REAL(wp) ::   zcv, zdy2, zdy4    !   -      - 
    446       REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv2, ztv3, ztv4 
     1083      REAL(wp), DIMENSION(jpi,jpj) ::   ztv1, ztv2, ztv3, ztv4 
    4471084      !!---------------------------------------------------------------------- 
    4481085      ! 
    … …  
    4741111      ! 
    4751112      ! 
    476       SELECT CASE (k_order ) 
     1113      SELECT CASE (kn_umx ) 
    4771114      ! 
    4781115      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21) 
    4791116         DO jj = 1, jpjm1 
    480             DO ji = fs_2, fs_jpim1 
    481                pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1) + pt(ji,jj) )  & 
    482                   &                                     - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 
     1117            DO ji = 1, fs_jpim1 
     1118               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  & 
     1119                  &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 
    4831120            END DO 
    4841121         END DO 
    … …  
    4861123      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23) 
    4871124         DO jj = 1, jpjm1 
    488             DO ji = fs_2, fs_jpim1 
    489                zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1125            DO ji = 1, fs_jpim1 
     1126               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    4901127               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (        ( pt(ji,jj+1) + pt(ji,jj) )  & 
    4911128                  &                                     - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) ) 
    … …  
    4961133      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24) 
    4971134         DO jj = 1, jpjm1 
    498             DO ji = fs_2, fs_jpim1 
    499                zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1135            DO ji = 1, fs_jpim1 
     1136               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    5001137               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    5011138!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    … …  
    5091146      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27) 
    5101147         DO jj = 1, jpjm1 
    511             DO ji = fs_2, fs_jpim1 
    512                zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1148            DO ji = 1, fs_jpim1 
     1149               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    5131150               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    5141151!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    … …  
    5221159      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29) 
    5231160         DO jj = 1, jpjm1 
    524             DO ji = fs_2, fs_jpim1 
    525                zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1161            DO ji = 1, fs_jpim1 
     1162               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    5261163               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    5271164!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    … …  
    5371174         ! 
    5381175      END SELECT 
     1176      !                                                     !-- High order flux in j-direction  --! 
     1177      IF( ll_neg ) THEN 
     1178         DO jj = 1, jpjm1 
     1179            DO ji = 1, fs_jpim1 
     1180               IF( pt_v(ji,jj) < 0._wp ) THEN 
     1181                  pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  & 
     1182                     &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 
     1183               ENDIF 
     1184            END DO 
     1185         END DO 
     1186      ENDIF 
     1187 
     1188      DO jj = 1, jpjm1 
     1189         DO ji = 1, fs_jpim1   ! vector opt. 
     1190            IF( ll_clem ) THEN 
     1191               pfv_ho(ji,jj) = pv(ji,jj) * pt_v(ji,jj) 
     1192            ELSE 
     1193               pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj) 
     1194            ENDIF 
     1195         END DO 
     1196      END DO 
    5391197      ! 
    5401198   END SUBROUTINE ultimate_y 
    541     
    542    
    543    SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, pdt ) 
     1199      
     1200 
     1201   SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho ) 
    5441202      !!--------------------------------------------------------------------- 
    5451203      !!                    ***  ROUTINE nonosc  *** 
    5461204      !!      
    547       !! **  Purpose :   compute monotonic tracer fluxes from the upstream  
     1205       !! **  Purpose :   compute monotonic tracer fluxes from the upstream  
    5481206      !!       scheme and the before field by a nonoscillatory algorithm  
    5491207      !! 
    5501208      !! **  Method  :   ... ??? 
    551       !!       warning : pbef and paft must be masked, but the boundaries 
     1209      !!       warning : pt and pt_low must be masked, but the boundaries 
    5521210      !!       conditions on the fluxes are not necessary zalezak (1979) 
    5531211      !!       drange (1995) multi-dimensional forward-in-time and upstream- 
    5541212      !!       in-space based differencing for fluid 
    5551213      !!---------------------------------------------------------------------- 
    556       REAL(wp)                     , INTENT(in   ) ::   pdt          ! tracer time-step 
    557       REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pbef, paft   ! before & after field 
    558       REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   paa, pbb     ! monotonic fluxes in the 2 directions 
     1214      REAL(wp)                     , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     1215      REAL(wp)                     , INTENT(in   ) ::   pdt              ! tracer time-step 
     1216      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pu               ! ice i-velocity => u*e2 
     1217      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   puc              ! ice i-velocity *A => u*e2*a 
     1218      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pv               ! ice j-velocity => v*e1 
     1219      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pvc              ! ice j-velocity *A => v*e1*a 
     1220      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt, pt_low  ! before field & upstream guess of after field 
     1221      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_low, pfu_low ! upstream flux 
     1222      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux 
    5591223      ! 
    5601224      INTEGER  ::   ji, jj    ! dummy loop indices 
    561       INTEGER  ::   ikm1      ! local integer 
    562       REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zsml, z1_dt   ! local scalars 
    563       REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      - 
    564       REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zmsk, zdiv 
    565       !!---------------------------------------------------------------------- 
    566       ! 
     1225      REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt, zpos2, zneg2   ! local scalars 
     1226      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign, zcoef        !   -      - 
     1227      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zti_low, ztj_low, zzt 
     1228      !!---------------------------------------------------------------------- 
    5671229      zbig = 1.e+40_wp 
    568       zsml = 1.e-15_wp 
    569  
    570       ! test on divergence 
    571       DO jj = 2, jpjm1 
    572          DO ji = fs_2, fs_jpim1   ! vector opt.   
    573             zdiv(ji,jj) =  - (  paa(ji,jj) - paa(ji-1,jj  )   & 
    574                &              + pbb(ji,jj) - pbb(ji  ,jj-1) )   
    575          END DO 
    576       END DO 
    577       CALL lbc_lnk( zdiv, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign) 
    578  
    579       ! Determine ice masks for before and after tracers  
    580       WHERE( pbef(:,:) == 0._wp .AND. paft(:,:) == 0._wp .AND. zdiv(:,:) == 0._wp )   ;   zmsk(:,:) = 0._wp 
    581       ELSEWHERE                                                                       ;   zmsk(:,:) = 1._wp * tmask(:,:,1) 
    582       END WHERE 
    583  
     1230      zsml = epsi20 
     1231 
     1232      IF( ll_zeroup2 ) THEN 
     1233         DO jj = 1, jpjm1 
     1234            DO ji = 1, fs_jpim1   ! vector opt. 
     1235               IF( amaxu(ji,jj) == 0._wp )   pfu_ho(ji,jj) = 0._wp 
     1236               IF( amaxv(ji,jj) == 0._wp )   pfv_ho(ji,jj) = 0._wp 
     1237            END DO 
     1238         END DO 
     1239      ENDIF 
     1240       
     1241      IF( ll_zeroup4 ) THEN 
     1242         DO jj = 1, jpjm1 
     1243            DO ji = 1, fs_jpim1   ! vector opt. 
     1244               IF( pfu_low(ji,jj) == 0._wp )   pfu_ho(ji,jj) = 0._wp 
     1245               IF( pfv_low(ji,jj) == 0._wp )   pfv_ho(ji,jj) = 0._wp 
     1246            END DO 
     1247         END DO 
     1248      ENDIF 
     1249 
     1250 
     1251      IF( ll_zeroup1 ) THEN 
     1252         DO jj = 2, jpjm1 
     1253            DO ji = fs_2, fs_jpim1 
     1254               IF( ll_gurvan ) THEN 
     1255                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1256                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1257               ELSE 
     1258                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1259                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  & 
     1260                     &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1261                     &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     1262               ENDIF 
     1263               IF( zzt(ji,jj) < 0._wp ) THEN 
     1264                  pfu_ho(ji,jj)   = pfu_low(ji,jj) 
     1265                  pfv_ho(ji,jj)   = pfv_low(ji,jj) 
     1266                  WRITE(numout,*) '*** 1 negative high order zzt ***',ji,jj,zzt(ji,jj) 
     1267               ENDIF 
     1268!!               IF( ji==26 .AND. jj==86) THEN 
     1269!!                  WRITE(numout,*) 'zzt high order',zzt(ji,jj) 
     1270!!                  WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1271!!                  WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1272!!                  WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt 
     1273!!                  WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt 
     1274!!               ENDIF 
     1275               IF( ll_gurvan ) THEN 
     1276                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1277                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1278               ELSE 
     1279                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1280                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  & 
     1281                     &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1282                     &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     1283               ENDIF 
     1284               IF( zzt(ji,jj) < 0._wp ) THEN 
     1285                  pfu_ho(ji-1,jj) = pfu_low(ji-1,jj) 
     1286                  pfv_ho(ji,jj-1) = pfv_low(ji,jj-1) 
     1287                  WRITE(numout,*) '*** 2 negative high order zzt ***',ji,jj,zzt(ji,jj) 
     1288               ENDIF 
     1289               IF( ll_gurvan ) THEN 
     1290                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1291                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1292               ELSE 
     1293                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1294                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  & 
     1295                     &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1296                     &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     1297               ENDIF 
     1298               IF( zzt(ji,jj) < 0._wp ) THEN 
     1299                  WRITE(numout,*) '*** 3 negative high order zzt ***',ji,jj,zzt(ji,jj) 
     1300               ENDIF 
     1301            END DO 
     1302         END DO 
     1303         CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) 
     1304      ENDIF 
     1305 
     1306       
     1307      ! antidiffusive flux : high order minus low order 
     1308      ! -------------------------------------------------- 
     1309      DO jj = 1, jpjm1 
     1310         DO ji = 1, fs_jpim1   ! vector opt. 
     1311            pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_low(ji,jj) 
     1312            pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_low(ji,jj) 
     1313          END DO 
     1314      END DO 
     1315 
     1316      ! extreme case where pfu_ho has to be zero 
     1317      ! ---------------------------------------- 
     1318      !                                    pfu_ho 
     1319      !                           *         ---> 
     1320      !                        |      |  *   |        |  
     1321      !                        |      |      |    *   |     
     1322      !                        |      |      |        |    * 
     1323      !            t_low :       i-1     i       i+1       i+2    
     1324      IF( ll_prelimiter_zalesak ) THEN 
     1325          
     1326         DO jj = 2, jpjm1 
     1327            DO ji = fs_2, fs_jpim1  
     1328               zti_low(ji,jj)= pt_low(ji+1,jj  ) 
     1329               ztj_low(ji,jj)= pt_low(ji  ,jj+1) 
     1330            END DO 
     1331         END DO 
     1332         CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. ) 
     1333 
     1334 
     1335         !! this does not work 
     1336         !!            DO jj = 2, jpjm1 
     1337         !!               DO ji = fs_2, fs_jpim1 
     1338         !!                  IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji+1,jj  ) - pt_low (ji  ,jj) ) .AND.     & 
     1339         !!               &      SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj+1) - pt_low (ji  ,jj) )           & 
     1340         !!               &    ) THEN 
     1341         !!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj ) - zti_low(ji  ,jj) ) .AND.   & 
     1342         !!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1 ) - ztj_low(ji  ,jj) )         & 
     1343         !!               &       ) THEN 
     1344         !!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0. 
     1345         !!                     ENDIF 
     1346         !!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji-1,jj  ) ) .AND.  & 
     1347         !!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji  ,jj-1) )        & 
     1348         !!               &       )  THEN 
     1349         !!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0. 
     1350         !!                     ENDIF 
     1351         !!                  ENDIF 
     1352         !!                END DO 
     1353         !!            END DO             
     1354 
     1355         DO jj = 2, jpjm1 
     1356            DO ji = fs_2, fs_jpim1 
     1357               IF ( pfu_ho(ji,jj) * ( pt_low(ji+1,jj) - pt_low(ji,jj) ) <= 0. .AND.  & 
     1358                  & pfv_ho(ji,jj) * ( pt_low(ji,jj+1) - pt_low(ji,jj) ) <= 0. ) THEN 
     1359                  ! 
     1360                  IF(  pfu_ho(ji,jj) * ( zti_low(ji+1,jj) - zti_low(ji,jj) ) <= 0 .AND.  & 
     1361                     & pfv_ho(ji,jj) * ( ztj_low(ji,jj+1) - ztj_low(ji,jj) ) <= 0)  pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0. 
     1362                  ! 
     1363                  IF(  pfu_ho(ji,jj) * ( pt_low(ji  ,jj) - pt_low(ji-1,jj) ) <= 0 .AND.  & 
     1364                     & pfv_ho(ji,jj) * ( pt_low(ji  ,jj) - pt_low(ji,jj-1) ) <= 0)  pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0. 
     1365                  ! 
     1366               ENDIF 
     1367            END DO 
     1368         END DO 
     1369         CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond. 
     1370 
     1371      ELSEIF( ll_prelimiter_devore ) THEN 
     1372         DO jj = 2, jpjm1 
     1373            DO ji = fs_2, fs_jpim1  
     1374               zti_low(ji,jj)= pt_low(ji+1,jj  ) 
     1375               ztj_low(ji,jj)= pt_low(ji  ,jj+1) 
     1376            END DO 
     1377         END DO 
     1378         CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. ) 
     1379 
     1380         z1_dt = 1._wp / pdt 
     1381         DO jj = 2, jpjm1 
     1382            DO ji = fs_2, fs_jpim1 
     1383               zsign = SIGN( 1., pt_low(ji+1,jj) - pt_low(ji,jj) ) 
     1384               pfu_ho(ji,jj) =  zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj)) , & 
     1385                  &                          zsign * ( pt_low (ji  ,jj) - pt_low (ji-1,jj) ) * e1e2t(ji  ,jj) * z1_dt , & 
     1386                  &                          zsign * ( zti_low(ji+1,jj) - zti_low(ji  ,jj) ) * e1e2t(ji+1,jj) * z1_dt ) ) 
     1387 
     1388               zsign = SIGN( 1., pt_low(ji,jj+1) - pt_low(ji,jj) ) 
     1389               pfv_ho(ji,jj) =  zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj)) , & 
     1390                  &                          zsign * ( pt_low (ji,jj  ) - pt_low (ji,jj-1) ) * e1e2t(ji,jj  ) * z1_dt , & 
     1391                  &                          zsign * ( ztj_low(ji,jj+1) - ztj_low(ji,jj  ) ) * e1e2t(ji,jj+1) * z1_dt ) ) 
     1392            END DO 
     1393         END DO 
     1394         CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond. 
     1395             
     1396      ENDIF 
     1397          
     1398       
    5841399      ! Search local extrema 
    5851400      ! -------------------- 
    586       ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 
    587 !      zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ),   & 
    588 !         &             paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) )  ) 
    589 !      zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ),   & 
    590 !         &             paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) )  ) 
    591       zbup(:,:) = MAX( pbef(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ),   & 
    592          &             paft(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) )  ) 
    593       zbdo(:,:) = MIN( pbef(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ),   & 
    594          &             paft(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) )  ) 
    595  
     1401      ! max/min of pt & pt_low with large negative/positive value (-/+zbig) outside ice cover 
     1402      DO jj = 1, jpj 
     1403         DO ji = 1, jpi 
     1404            IF    ( pt(ji,jj) <= 0._wp .AND. pt_low(ji,jj) <= 0._wp ) THEN 
     1405               zbup(ji,jj) = -zbig 
     1406               zbdo(ji,jj) =  zbig 
     1407            ELSEIF( pt(ji,jj) <= 0._wp .AND. pt_low(ji,jj) > 0._wp ) THEN 
     1408               zbup(ji,jj) = pt_low(ji,jj) 
     1409               zbdo(ji,jj) = pt_low(ji,jj) 
     1410            ELSEIF( pt(ji,jj) > 0._wp .AND. pt_low(ji,jj) <= 0._wp ) THEN 
     1411               zbup(ji,jj) = pt(ji,jj) 
     1412               zbdo(ji,jj) = pt(ji,jj) 
     1413            ELSE 
     1414               zbup(ji,jj) = MAX( pt(ji,jj) , pt_low(ji,jj) ) 
     1415               zbdo(ji,jj) = MIN( pt(ji,jj) , pt_low(ji,jj) ) 
     1416            ENDIF 
     1417         END DO 
     1418      END DO 
     1419 
     1420  
    5961421      z1_dt = 1._wp / pdt 
    5971422      DO jj = 2, jpjm1 
    5981423         DO ji = fs_2, fs_jpim1   ! vector opt. 
    5991424            ! 
    600             zup  = MAX(   zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ),   &        ! search max/min in neighbourhood 
    601                &                       zbup(ji  ,jj-1), zbup(ji  ,jj+1)    ) 
    602             zdo  = MIN(   zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ),   & 
    603                &                       zbdo(ji  ,jj-1), zbdo(ji  ,jj+1)    ) 
    604                ! 
    605             zpos = MAX( 0., paa(ji-1,jj  ) ) - MIN( 0., paa(ji  ,jj  ) )   &        ! positive/negative  part of the flux 
    606                & + MAX( 0., pbb(ji  ,jj-1) ) - MIN( 0., pbb(ji  ,jj  ) ) 
    607             zneg = MAX( 0., paa(ji  ,jj  ) ) - MIN( 0., paa(ji-1,jj  ) )   & 
    608                & + MAX( 0., pbb(ji  ,jj  ) ) - MIN( 0., pbb(ji  ,jj-1) ) 
    609                ! 
    610             zbt = e1e2t(ji,jj) * z1_dt                                   ! up & down beta terms 
    611             zbetup(ji,jj) = ( zup         - paft(ji,jj) ) / ( zpos + zsml ) * zbt 
    612             zbetdo(ji,jj) = ( paft(ji,jj) - zdo         ) / ( zneg + zsml ) * zbt 
     1425            IF( .NOT. ll_9points ) THEN 
     1426            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 
     1427            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1) ) 
     1428            ! 
     1429            ELSE 
     1430            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 
     1431               &                     zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1)  ) 
     1432            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1), & 
     1433               &                     zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1)  ) 
     1434            ENDIF 
     1435            ! 
     1436            zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji  ,jj) ) &  ! positive/negative part of the flux 
     1437               & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj  ) ) 
     1438            zneg = MAX( 0., pfu_ho(ji  ,jj) ) - MIN( 0., pfu_ho(ji-1,jj) ) & 
     1439               & + MAX( 0., pfv_ho(ji,jj  ) ) - MIN( 0., pfv_ho(ji,jj-1) ) 
     1440            ! 
     1441            IF( ll_HgradU .AND. .NOT.ll_gurvan ) THEN 
     1442               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 ) 
     1443               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 ) 
     1444            ELSE 
     1445               zneg2 = 0. ; zpos2 = 0. 
     1446            ENDIF 
     1447            ! 
     1448            !                                  ! up & down beta terms 
     1449!            zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt 
     1450!            zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt 
     1451 
     1452            IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj) = MAX( 0._wp, zup - pt_low(ji,jj) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt 
     1453            ELSE                         ; zbetup(ji,jj) = 0. ! zbig 
     1454            ENDIF 
     1455            ! 
     1456            IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj) = MAX( 0._wp, pt_low(ji,jj) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt 
     1457            ELSE                         ; zbetdo(ji,jj) = 0. ! zbig 
     1458            ENDIF 
     1459            ! 
     1460            ! if all the points are outside ice cover 
     1461            IF( zup == -zbig )   zbetup(ji,jj) = 0. ! zbig 
     1462            IF( zdo ==  zbig )   zbetdo(ji,jj) = 0. ! zbig             
     1463            ! 
     1464 
     1465!!            IF( ji==26 .AND. jj==86) THEN 
     1466!               WRITE(numout,*) '-----------------' 
     1467!               WRITE(numout,*) 'zpos',zpos,zpos2 
     1468!               WRITE(numout,*) 'zneg',zneg,zneg2 
     1469!               WRITE(numout,*) 'puc/pu',ABS(puc(ji,jj))/MAX(epsi20, ABS(pu(ji,jj))) 
     1470!               WRITE(numout,*) 'pvc/pv',ABS(pvc(ji,jj))/MAX(epsi20, ABS(pv(ji,jj))) 
     1471!               WRITE(numout,*) 'pucm1/pu',ABS(puc(ji-1,jj))/MAX(epsi20, ABS(pu(ji-1,jj))) 
     1472!               WRITE(numout,*) 'pvcm1/pv',ABS(pvc(ji,jj-1))/MAX(epsi20, ABS(pv(ji,jj-1))) 
     1473!               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)+pfu_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1474!               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)+pfv_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1475!               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)+pfu_low(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt 
     1476!               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)+pfv_low(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt 
     1477!               WRITE(numout,*) 'pfu_low',pfu_low(ji,jj) * r1_e1e2t(ji,jj) * pdt 
     1478!               WRITE(numout,*) 'pfv_low',pfv_low(ji,jj) * r1_e1e2t(ji,jj) * pdt 
     1479!               WRITE(numout,*) 'pfu_lowm1',pfu_low(ji-1,jj) * r1_e1e2t(ji,jj) * pdt 
     1480!               WRITE(numout,*) 'pfv_lowm1',pfv_low(ji,jj-1) * r1_e1e2t(ji,jj) * pdt 
     1481!                
     1482!               WRITE(numout,*) 'pt',pt(ji,jj) 
     1483!               WRITE(numout,*) 'ptim1',pt(ji-1,jj) 
     1484!               WRITE(numout,*) 'ptjm1',pt(ji,jj-1) 
     1485!               WRITE(numout,*) 'pt_low',pt_low(ji,jj) 
     1486!               WRITE(numout,*) 'zbetup',zbetup(ji,jj) 
     1487!               WRITE(numout,*) 'zbetdo',zbetdo(ji,jj) 
     1488!               WRITE(numout,*) 'zup',zup 
     1489!               WRITE(numout,*) 'zdo',zdo 
     1490!            ENDIF 
     1491            ! 
    6131492         END DO 
    6141493      END DO 
    6151494      CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign) 
    6161495 
    617       ! monotonic flux in the i & j direction (paa & pbb) 
    618       ! ------------------------------------- 
    619       DO jj = 2, jpjm1 
     1496       
     1497      ! monotonic flux in the y direction 
     1498      ! --------------------------------- 
     1499      DO jj = 1, jpjm1 
    6201500         DO ji = 1, fs_jpim1   ! vector opt. 
    6211501            zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) ) 
    6221502            zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) ) 
    623             zcu = 0.5  + SIGN( 0.5 , paa(ji,jj) ) 
    624             ! 
    625             paa(ji,jj) = paa(ji,jj) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 
    626          END DO 
    627       END DO 
    628       ! 
     1503            zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj) ) 
     1504            ! 
     1505            zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu ) 
     1506             
     1507            pfu_ho(ji,jj) = pfu_ho(ji,jj) * zcoef + pfu_low(ji,jj) 
     1508 
     1509!!            IF( ji==26 .AND. jj==86) THEN 
     1510!!               WRITE(numout,*) 'coefU',zcoef 
     1511!!               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1512!!               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt 
     1513!!            ENDIF 
     1514 
     1515         END DO 
     1516      END DO 
     1517 
    6291518      DO jj = 1, jpjm1 
    630          DO ji = fs_2, fs_jpim1   ! vector opt. 
     1519         DO ji = 1, fs_jpim1   ! vector opt. 
    6311520            zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) ) 
    6321521            zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) ) 
    633             zcv = 0.5  + SIGN( 0.5 , pbb(ji,jj) ) 
    634             ! 
    635             pbb(ji,jj) = pbb(ji,jj) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 
    636          END DO 
    637       END DO 
     1522            zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj) ) 
     1523            ! 
     1524            zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv ) 
     1525             
     1526            pfv_ho(ji,jj) = pfv_ho(ji,jj) * zcoef + pfv_low(ji,jj) 
     1527 
     1528!!            IF( ji==26 .AND. jj==86) THEN 
     1529!!               WRITE(numout,*) 'coefV',zcoef 
     1530!!               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1531!!               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt 
     1532!!            ENDIF 
     1533         END DO 
     1534      END DO 
     1535 
     1536      ! clem test 
     1537      DO jj = 2, jpjm1 
     1538         DO ji = 2, fs_jpim1   ! vector opt. 
     1539            IF( ll_gurvan ) THEN 
     1540               zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1541                  &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1542            ELSE 
     1543               zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  & 
     1544                  &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  & 
     1545                  &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1546                  &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     1547            ENDIF 
     1548            IF( zzt(ji,jj) < -epsi20 ) THEN 
     1549               WRITE(numout,*) 'T<0 nonosc',zzt(ji,jj) 
     1550            ENDIF 
     1551         END DO 
     1552      END DO 
     1553       
     1554      ! 
    6381555      ! 
    6391556   END SUBROUTINE nonosc_2d 
     1557 
     1558   SUBROUTINE limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     1559      !!--------------------------------------------------------------------- 
     1560      !!                    ***  ROUTINE limiter_x  *** 
     1561      !!      
     1562      !! **  Purpose :   compute flux limiter  
     1563      !!---------------------------------------------------------------------- 
     1564      REAL(wp)                     , INTENT(in   )           ::   pdt          ! tracer time-step 
     1565      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pu           ! ice i-velocity => u*e2 
     1566      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   puc          ! ice i-velocity *A => u*e2*a 
     1567      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pt           ! ice tracer 
     1568      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout)           ::   pfu_ho       ! high order flux 
     1569      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ), OPTIONAL ::   pfu_ups      ! upstream flux 
     1570      ! 
     1571      REAL(wp) ::   Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr 
     1572      INTEGER  ::   ji, jj    ! dummy loop indices 
     1573      REAL(wp), DIMENSION (jpi,jpj) ::   zslpx       ! tracer slopes  
     1574      !!---------------------------------------------------------------------- 
     1575      ! 
     1576      DO jj = 2, jpjm1 
     1577         DO ji = fs_2, fs_jpim1   ! vector opt. 
     1578            zslpx(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * umask(ji,jj,1) 
     1579         END DO 
     1580      END DO 
     1581      CALL lbc_lnk( zslpx, 'U', -1.)   ! lateral boundary cond. 
     1582       
     1583      DO jj = 2, jpjm1 
     1584         DO ji = fs_2, fs_jpim1   ! vector opt. 
     1585            uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) 
     1586 
     1587            Rjm = zslpx(ji-1,jj) 
     1588            Rj  = zslpx(ji  ,jj) 
     1589            Rjp = zslpx(ji+1,jj) 
     1590 
     1591            IF( PRESENT(pfu_ups) ) THEN 
     1592 
     1593               IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
     1594               ELSE                        ;   Rr = Rjp 
     1595               ENDIF 
     1596 
     1597               zh3 = pfu_ho(ji,jj) - pfu_ups(ji,jj)      
     1598               IF( Rj > 0. ) THEN 
     1599                  zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(puc(ji,jj)),  & 
     1600                     &        MIN( 2. * Rr * 0.5 * ABS(puc(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) ) 
     1601               ELSE 
     1602                  zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(puc(ji,jj)),  & 
     1603                     &        MIN(-2. * Rr * 0.5 * ABS(puc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) ) 
     1604               ENDIF 
     1605               pfu_ho(ji,jj) = pfu_ups(ji,jj) + zlimiter 
     1606                
     1607            ELSE 
     1608               IF( Rj /= 0. ) THEN 
     1609                  IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj 
     1610                  ELSE                        ;   Cr = Rjp / Rj 
     1611                  ENDIF 
     1612               ELSE 
     1613                  Cr = 0. 
     1614                  !IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20 
     1615                  !ELSE                        ;   Cr = Rjp * 1.e20 
     1616                  !ENDIF 
     1617               ENDIF 
     1618                
     1619               ! -- superbee -- 
     1620               zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 
     1621               ! -- van albada 2 -- 
     1622               !!zpsi = 2.*Cr / (Cr*Cr+1.) 
     1623 
     1624               ! -- sweby (with beta=1) -- 
     1625               !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 
     1626               ! -- van Leer -- 
     1627               !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 
     1628               ! -- ospre -- 
     1629               !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 
     1630               ! -- koren -- 
     1631               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 
     1632               ! -- charm -- 
     1633               !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 
     1634               !ELSE                 ;   zpsi = 0. 
     1635               !ENDIF 
     1636               ! -- van albada 1 -- 
     1637               !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 
     1638               ! -- smart -- 
     1639               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 
     1640               ! -- umist -- 
     1641               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 
     1642 
     1643               ! high order flux corrected by the limiter 
     1644               pfu_ho(ji,jj) = pfu_ho(ji,jj) - ABS( puc(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5 
     1645 
     1646            ENDIF 
     1647         END DO 
     1648      END DO 
     1649      CALL lbc_lnk( pfu_ho, 'U', -1.)   ! lateral boundary cond. 
     1650      ! 
     1651   END SUBROUTINE limiter_x 
     1652 
     1653   SUBROUTINE limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     1654      !!--------------------------------------------------------------------- 
     1655      !!                    ***  ROUTINE limiter_y  *** 
     1656      !!      
     1657      !! **  Purpose :   compute flux limiter  
     1658      !!---------------------------------------------------------------------- 
     1659      REAL(wp)                     , INTENT(in   )           ::   pdt          ! tracer time-step 
     1660      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pv           ! ice i-velocity => u*e2 
     1661      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pvc          ! ice i-velocity *A => u*e2*a 
     1662      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pt           ! ice tracer 
     1663      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout)           ::   pfv_ho       ! high order flux 
     1664      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ), OPTIONAL ::   pfv_ups      ! upstream flux 
     1665      ! 
     1666      REAL(wp) ::   Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr 
     1667      INTEGER  ::   ji, jj    ! dummy loop indices 
     1668      REAL(wp), DIMENSION (jpi,jpj) ::   zslpy       ! tracer slopes  
     1669      !!---------------------------------------------------------------------- 
     1670      ! 
     1671      DO jj = 2, jpjm1 
     1672         DO ji = fs_2, fs_jpim1   ! vector opt. 
     1673            zslpy(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * vmask(ji,jj,1) 
     1674         END DO 
     1675      END DO 
     1676      CALL lbc_lnk( zslpy, 'V', -1.)   ! lateral boundary cond. 
     1677       
     1678      DO jj = 2, jpjm1 
     1679         DO ji = fs_2, fs_jpim1   ! vector opt. 
     1680            vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) 
     1681 
     1682            Rjm = zslpy(ji,jj-1) 
     1683            Rj  = zslpy(ji,jj  ) 
     1684            Rjp = zslpy(ji,jj+1) 
     1685 
     1686            IF( PRESENT(pfv_ups) ) THEN 
     1687 
     1688               IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
     1689               ELSE                        ;   Rr = Rjp 
     1690               ENDIF 
     1691 
     1692               zh3 = pfv_ho(ji,jj) - pfv_ups(ji,jj)      
     1693               IF( Rj > 0. ) THEN 
     1694                  zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pvc(ji,jj)),  & 
     1695                     &        MIN( 2. * Rr * 0.5 * ABS(pvc(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) ) 
     1696               ELSE 
     1697                  zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pvc(ji,jj)),  & 
     1698                     &        MIN(-2. * Rr * 0.5 * ABS(pvc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) ) 
     1699               ENDIF 
     1700               pfv_ho(ji,jj) = pfv_ups(ji,jj) + zlimiter 
     1701 
     1702            ELSE 
     1703 
     1704               IF( Rj /= 0. ) THEN 
     1705                  IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj 
     1706                  ELSE                        ;   Cr = Rjp / Rj 
     1707                  ENDIF 
     1708               ELSE 
     1709                  Cr = 0. 
     1710                  !IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20 
     1711                  !ELSE                        ;   Cr = Rjp * 1.e20 
     1712                  !ENDIF 
     1713               ENDIF 
     1714 
     1715               ! -- superbee -- 
     1716               zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 
     1717               ! -- van albada 2 -- 
     1718               !!zpsi = 2.*Cr / (Cr*Cr+1.) 
     1719 
     1720               ! -- sweby (with beta=1) -- 
     1721               !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 
     1722               ! -- van Leer -- 
     1723               !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 
     1724               ! -- ospre -- 
     1725               !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 
     1726               ! -- koren -- 
     1727               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 
     1728               ! -- charm -- 
     1729               !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 
     1730               !ELSE                 ;   zpsi = 0. 
     1731               !ENDIF 
     1732               ! -- van albada 1 -- 
     1733               !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 
     1734               ! -- smart -- 
     1735               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 
     1736               ! -- umist -- 
     1737               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 
     1738 
     1739               ! high order flux corrected by the limiter 
     1740               pfv_ho(ji,jj) = pfv_ho(ji,jj) - ABS( pvc(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5 
     1741 
     1742            ENDIF 
     1743         END DO 
     1744      END DO 
     1745      CALL lbc_lnk( pfv_ho, 'V', -1.)   ! lateral boundary cond. 
     1746      ! 
     1747   END SUBROUTINE limiter_y 
    6401748 
    6411749#else 

  • NEMO/trunk/src/ICE/icedyn_rdgrft.F90

    r10069 r10413  
    3939 
    4040   ! Variables shared among ridging subroutines 
    41    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   closing_net     ! net rate at which area is removed    (1/s) 
    42       !                                                 ! (ridging ice area - area of new ridges) / dt 
    43    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   opning         ! rate of opening due to divergence/shear 
    44    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   closing_gross  ! rate at which area removed, not counting area of new ridges 
    45    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   apartf   ! participation function; fraction of ridging/closing associated w/ category n 
    46    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hrmin    ! minimum ridge thickness 
    47    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hrmax    ! maximum ridge thickness 
    48    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hraft    ! thickness of rafted ice 
    49    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hi_hrdg  ! thickness of ridging ice / mean ridge thickness 
    50    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   aridge   ! participating ice ridging 
    51    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   araft    ! participating ice rafting 
     41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   closing_net     ! net rate at which area is removed    (1/s) 
     42      !                                                               ! (ridging ice area - area of new ridges) / dt 
     43   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   opning          ! rate of opening due to divergence/shear 
     44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   closing_gross   ! rate at which area removed, not counting area of new ridges 
     45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   apartf          ! participation function; fraction of ridging/closing associated w/ category n 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hrmin           ! minimum ridge thickness 
     47   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hrmax           ! maximum ridge thickness 
     48   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hraft           ! thickness of rafted ice 
     49   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hi_hrdg         ! thickness of ridging ice / mean ridge thickness 
     50   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aridge          ! participating ice ridging 
     51   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   araft           ! participating ice rafting 
    5252   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_i_2d 
    5353   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_s_2d 
    … …  
    5959   LOGICAL  ::   ln_str_H79       ! ice strength parameterization (Hibler79) 
    6060   REAL(wp) ::   rn_pstar         ! determines ice strength, Hibler JPO79 
    61    REAL(wp) ::   rn_crhg          ! determines changes in ice strength 
    6261   REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging             
    6362   LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975)) 
    … …  
    7978   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    8079   !! $Id$ 
    81    !! Software governed by the CeCILL license (see ./LICENSE) 
     80   !! Software governed by the CeCILL licence     (./LICENSE) 
    8281   !!---------------------------------------------------------------------- 
    8382CONTAINS 
    … …  
    193192            ! divergence given by the advection scheme 
    194193            !   (which may not be equal to divu as computed from the velocity field) 
    195             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 
    196199            ! 
    197200            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/trunk/src/ICE/icedyn_rhg.F90

    r10069 r10413  
    4343   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    4444   !! $Id$ 
    45    !! Software governed by the CeCILL license (see ./LICENSE) 
     45   !! Software governed by the CeCILL licence     (./LICENSE) 
    4646   !!---------------------------------------------------------------------- 
    4747CONTAINS 

  • NEMO/trunk/src/ICE/icedyn_rhg_evp.F90

    r10332 r10413  
    119119      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity 
    120120      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017 
    121       REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                       ! ice/snow mass 
     121      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume 
    122122      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars 
    123123      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                           ! temporary scalars 
     124      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice 
     125      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast 
    124126      ! 
    125127      REAL(wp) ::   zresm                                               ! Maximal error on ice velocity 
    … …  
    137139      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points 
    138140      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points 
     141      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ib , ztauV_ib             ! ice-bottom stress at U-V points (landfast param) 
    139142      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points 
    140143      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                            
    … …  
    256259         END DO 
    257260      END DO 
    258              
     261 
     262      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 
     263      IF( ln_landfast_L16 .OR. ln_landfast_home ) THEN   ;   zkt = rn_tensile 
     264      ELSE                                               ;   zkt = 0._wp 
     265      ENDIF 
    259266      ! 
    260267      !------------------------------------------------------------------------------! 
    … …  
    317324      CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. ) 
    318325      ! 
     326      !                                  !== Landfast ice parameterization ==! 
     327      ! 
     328      IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016 
     329         DO jj = 2, jpjm1 
     330            DO ji = fs_2, fs_jpim1 
     331               ! ice thickness at U-V points 
     332               zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
     333               zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
     334               ! ice-bottom stress at U points 
     335               zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj) 
     336               zTauU_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
     337               ! ice-bottom stress at V points 
     338               zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj) 
     339               zTauV_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
     340               ! ice_bottom stress at T points 
     341               zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj) 
     342               tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
     343            END DO 
     344         END DO 
     345         CALL lbc_lnk( tau_icebfr(:,:), 'T', 1. ) 
     346         ! 
     347      ELSEIF( ln_landfast_home ) THEN          !-- Home made 
     348         DO jj = 2, jpjm1 
     349            DO ji = fs_2, fs_jpim1 
     350               zTauU_ib(ji,jj) = tau_icebfr(ji,jj) 
     351               zTauV_ib(ji,jj) = tau_icebfr(ji,jj) 
     352            END DO 
     353         END DO 
     354         ! 
     355      ELSE                                     !-- no landfast 
     356         DO jj = 2, jpjm1 
     357            DO ji = fs_2, fs_jpim1 
     358               zTauU_ib(ji,jj) = 0._wp 
     359               zTauV_ib(ji,jj) = 0._wp 
     360            END DO 
     361         END DO 
     362      ENDIF 
     363      IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) ) 
     364 
    319365      !------------------------------------------------------------------------------! 
    320366      ! 3) Solution of the momentum equation, iterative procedure 
    … …  
    380426               ENDIF 
    381427                
    382                ! stress at T points 
    383                zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1 
    384                zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2 
     428               ! stress at T points (zkt/=0 if landfast) 
     429               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 
     430               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 
    385431              
    386432            END DO 
    … …  
    401447               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 
    402448                
    403                ! stress at F points 
    404                zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2 
     449               ! stress at F points (zkt/=0 if landfast) 
     450               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 
    405451 
    406452            END DO 
    … …  
    450496                  ! 
    451497                  !                 !--- tau_bottom/v_ice 
    452                   zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
    453                   zTauB = - tau_icebfr(ji,jj) / zvel 
     498                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
     499                  zTauB = - zTauV_ib(ji,jj) / zvel 
    454500                  ! 
    455501                  !                 !--- Coriolis at V-points (energy conserving formulation) 
    … …  
    462508                  ! 
    463509                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    464                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     510                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    465511                  ! 
    466512                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    … …  
    498544                  ! 
    499545                  !                 !--- tau_bottom/u_ice 
    500                   zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
    501                   zTauB = - tau_icebfr(ji,jj) / zvel 
     546                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
     547                  zTauB = - zTauU_ib(ji,jj) / zvel 
    502548                  ! 
    503549                  !                 !--- Coriolis at U-points (energy conserving formulation) 
    … …  
    510556                  ! 
    511557                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    512                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     558                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    513559                  ! 
    514560                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    … …  
    548594                  ! 
    549595                  !                 !--- tau_bottom/u_ice 
    550                   zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
    551                   zTauB = - tau_icebfr(ji,jj) / zvel 
     596                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
     597                  zTauB = - zTauU_ib(ji,jj) / zvel 
    552598                  ! 
    553599                  !                 !--- Coriolis at U-points (energy conserving formulation) 
    … …  
    560606                  ! 
    561607                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    562                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     608                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    563609                  ! 
    564610                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    … …  
    596642                  ! 
    597643                  !                 !--- tau_bottom/v_ice 
    598                   zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
    599                   ztauB = - tau_icebfr(ji,jj) / zvel 
     644                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
     645                  zTauB = - zTauV_ib(ji,jj) / zvel 
    600646                  ! 
    601647                  !                 !--- Coriolis at v-points (energy conserving formulation) 
    … …  
    608654                  ! 
    609655                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    610                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     656                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    611657                  ! 
    612658                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 

  • NEMO/trunk/src/ICE/iceistate.F90

    r10069 r10413  
    6464   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    6565   !! $Id$ 
    66    !! Software governed by the CeCILL license (see ./LICENSE) 
     66   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    6767   !!---------------------------------------------------------------------- 
    6868CONTAINS 
    … …  
    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/trunk/src/ICE/icevar.F90

    r10332 r10413  
    557557      !!------------------------------------------------------------------- 
    558558      ! 
     559      ! 
     560      DO jl = 1, jpl       !==  loop over the categories  ==! 
     561         ! 
     562         !---------------------------------------- 
     563         ! zap ice energy and send it to the ocean 
     564         !---------------------------------------- 
     565         DO jk = 1, nlay_i 
     566            DO jj = 1 , jpj 
     567               DO ji = 1 , jpi 
     568                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
     569                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0 
     570                     pe_i(ji,jj,jk,jl) = 0._wp 
     571                  ENDIF 
     572               END DO 
     573            END DO 
     574         END DO 
     575         ! 
     576         DO jk = 1, nlay_s 
     577            DO jj = 1 , jpj 
     578               DO ji = 1 , jpi 
     579                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
     580                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0 
     581                     pe_s(ji,jj,jk,jl) = 0._wp 
     582                  ENDIF 
     583               END DO 
     584            END DO 
     585         END DO 
     586         ! 
     587         !----------------------------------------------------- 
     588         ! zap ice and snow volume, add water and salt to ocean 
     589         !----------------------------------------------------- 
     590         DO jj = 1 , jpj 
     591            DO ji = 1 , jpi 
     592               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
     593                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * r1_rdtice 
     594                  pv_i   (ji,jj,jl) = 0._wp 
     595               ENDIF 
     596               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
     597                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * r1_rdtice 
     598                  pv_s   (ji,jj,jl) = 0._wp 
     599               ENDIF 
     600               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
     601                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * r1_rdtice 
     602                  psv_i  (ji,jj,jl) = 0._wp 
     603               ENDIF 
     604            END DO 
     605         END DO 
     606         ! 
     607      END DO  
     608      ! 
    559609      WHERE( pato_i(:,:)   < 0._wp )   pato_i(:,:)   = 0._wp 
    560610      WHERE( poa_i (:,:,:) < 0._wp )   poa_i (:,:,:) = 0._wp 
    … …  
    563613      WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+) 
    564614      !                                                        but it does not change conservation, so keep it this way is ok 
    565       ! 
    566       DO jl = 1, jpl       !==  loop over the categories  ==! 
    567          ! 
    568          !---------------------------------------- 
    569          ! zap ice energy and send it to the ocean 
    570          !---------------------------------------- 
    571          DO jk = 1, nlay_i 
    572             DO jj = 1 , jpj 
    573                DO ji = 1 , jpi 
    574                   IF( pe_i(ji,jj,jk,jl) < 0._wp ) THEN 
    575                      hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0 
    576                      pe_i(ji,jj,jk,jl) = 0._wp 
    577                   ENDIF 
    578                END DO 
    579             END DO 
    580          END DO 
    581          ! 
    582          DO jk = 1, nlay_s 
    583             DO jj = 1 , jpj 
    584                DO ji = 1 , jpi 
    585                   IF( pe_s(ji,jj,jk,jl) < 0._wp ) THEN 
    586                      hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0 
    587                      pe_s(ji,jj,jk,jl) = 0._wp 
    588                   ENDIF 
    589                END DO 
    590             END DO 
    591          END DO 
    592          ! 
    593          !----------------------------------------------------- 
    594          ! zap ice and snow volume, add water and salt to ocean 
    595          !----------------------------------------------------- 
    596          DO jj = 1 , jpj 
    597             DO ji = 1 , jpi 
    598               IF( pv_i(ji,jj,jl) < 0._wp ) THEN 
    599                   wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * r1_rdtice 
    600                   pv_i   (ji,jj,jl) = 0._wp 
    601                ENDIF 
    602                IF( pv_s(ji,jj,jl) < 0._wp ) THEN 
    603                   wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * r1_rdtice 
    604                   pv_s   (ji,jj,jl) = 0._wp 
    605                ENDIF 
    606                IF( psv_i(ji,jj,jl) < 0._wp ) THEN 
    607                   sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * r1_rdtice 
    608                   psv_i  (ji,jj,jl) = 0._wp 
    609                ENDIF 
    610             END DO 
    611          END DO 
    612          ! 
    613       END DO  
    614615      ! 
    615616   END SUBROUTINE ice_var_zapneg 

  • NEMO/trunk/src/ICE/icewri.F90

    r10069 r10413  
    3838   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    3939   !! $Id$ 
    40    !! Software governed by the CeCILL license (see ./LICENSE) 
     40   !! Software governed by the CeCILL licence     (./LICENSE) 
    4141   !!---------------------------------------------------------------------- 
    4242CONTAINS 
    … …  
    5050      INTEGER  ::   ji, jj, jk, jl  ! dummy loop indices 
    5151      REAL(wp) ::   z2da, z2db, zrho1, zrho2 
    52       REAL(wp), DIMENSION(jpi,jpj)     ::   z2d !  2D workspace 
     52      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d, zfast !  2D workspace 
    5353      REAL(wp), DIMENSION(jpi,jpj)     ::   zmsk00, zmsk05, zmsk15, zmsksn ! O%, 5% and 15% concentration mask and snow mask 
    5454      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zmsk00l, zmsksnl               ! cat masks 
    … …  
    132132      IF( iom_use('vtau_ai'  ) )   CALL iom_put( "vtau_ai", vtau_ice * zmsk00     )   ! Wind stress term in force balance (y) 
    133133 
    134       IF( iom_use('icevel') ) THEN  
     134      IF( iom_use('icevel') .OR. iom_use('fasticepres') ) THEN  
     135        ! module of ice velocity 
    135136         DO jj = 2 , jpjm1 
    136137            DO ji = 2 , jpim1 
    … …  
    141142         END DO 
    142143         CALL lbc_lnk( z2d, 'T', 1. ) 
    143          IF( iom_use('icevel') )   CALL iom_put( "icevel" , z2d                   )   ! ice velocity module 
     144         IF( iom_use('icevel') )   CALL iom_put( "icevel" , z2d ) 
     145 
     146        ! record presence of fast ice 
     147         WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk00(:,:) == 1._wp ) ; zfast(:,:) = 1._wp 
     148         ELSEWHERE                                                ; zfast(:,:) = 0._wp 
     149         END WHERE 
     150         IF( iom_use('fasticepres') )   CALL iom_put( "fasticepres" , zfast ) 
    144151      ENDIF 
    145152 

  • NEMO/trunk/src/NST/agrif_oce.F90

    r10068 r10413  
    2222   !                                              !!* Namelist namagrif: AGRIF parameters 
    2323   LOGICAL , PUBLIC ::   ln_spc_dyn    = .FALSE.   !: 
    24    INTEGER , PUBLIC, PARAMETER ::   nn_sponge_len = 2  !: Sponge width (in number of parent grid points) 
     24   INTEGER , PUBLIC, PARAMETER ::   nn_sponge_len = 1  !: Sponge width (in number of parent grid points) 
    2525   REAL(wp), PUBLIC ::   rn_sponge_tra = 2800.     !: sponge coeff. for tracers 
    2626   REAL(wp), PUBLIC ::   rn_sponge_dyn = 2800.     !: sponge coeff. for dynamics 

  • NEMO/trunk/tests/ICEDYN/EXPREF/1_namelist_ice_cfg

    r9801 r10413  
    11!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    2 !! SI3 namelist:   
     2!! SI3 configuration namelist: Overwrites SHARED/namelist_ice_ref 
    33!!              1 - Generic parameters                 (nampar) 
    44!!              2 - Ice thickness discretization       (namitd) 
    … …  
    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 

  • NEMO/trunk/tests/ICEDYN/EXPREF/file_def_nemo-ice.xml

    r9740 r10413  
    2121        <field field_ref="snwvolu"          name="snvolu" /> 
    2222        <field field_ref="icethic"          name="sithic" /> 
     23        <field field_ref="icethic"          name="sithic_max" operation="maximum" /> 
     24        <field field_ref="icethic"          name="sithic_min" operation="minimum" /> 
     25        <field field_ref="iceneg_pres"      name="sineg_pres" /> 
     26        <field field_ref="iceneg_volu"      name="sineg_volu" /> 
     27        <field field_ref="fasticepres"      name="fasticepres" /> 
    2328        <field field_ref="icevolu"          name="sivolu" /> 
    2429        <field field_ref="iceconc"          name="siconc" /> 

  • NEMO/trunk/tests/ICEDYN/EXPREF/namelist_ice_cfg

    r10075 r10413  
    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 
    … …  
    9596   sn_tmi = 'initice'                 , -12 ,'tmi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    9697   sn_smi = 'initice'                 , -12 ,'smi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    97    cn_dir ='./' 
     98   cn_dir='./' 
    9899/ 
    99100!------------------------------------------------------------------------------ 

  • NEMO/trunk/tests/demo_cfgs.txt

    r9789 r10413  
    44OVERFLOW OCE 
    55ICEDYN OCE NST SAS ICE 
     6ICEADV OCE SAS ICE 
    67VORTEX OCE NST 
    78WAD OCE 
     9 
Note: See TracChangeset for help on using the changeset viewer.