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 10419 for NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE – NEMO

Ignore:
Timestamp:
2018-12-19T20:46:30+01:00 (5 years ago)
Author:
smasson
Message:

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10418, see #2133

Location:
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE
Files:
12 edited

Legend:

Unmodified
Added
Removed

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

    r10292 r10419  
    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) ** 
    … …  
    325334   !! * Old values of global variables 
    326335   !!---------------------------------------------------------------------- 
    327    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b    !: snow and ice volumes/thickness 
    328    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b, oa_i_b         !: 
    329    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                         !: snow heat content 
    330    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                         !: ice temperatures 
    331    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b              !: ice velocity 
    332    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                        !: ice concentration (total) 
     336   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b, h_ip_b    !: snow and ice volumes/thickness 
     337   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b, oa_i_b                 !: 
     338   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                                 !: snow heat content 
     339   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                                 !: ice temperatures 
     340   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b                      !: ice velocity 
     341   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                                !: ice concentration (total) 
    333342             
    334343   !!---------------------------------------------------------------------- 
    … …  
    434443      ! * Old values of global variables 
    435444      ii = ii + 1 
    436       ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl)        ,   & 
    437          &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,   & 
     445      ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl), h_ip_b(jpi,jpj,jpl),  & 
     446         &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,               & 
    438447         &      oa_i_b(jpi,jpj,jpl)                                                   , STAT=ierr(ii) ) 
    439448 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icectl.F90

    r10314 r10419  
    197197 
    198198      ! heat flux 
    199       zhfx  = glob_sum( 'icectl', ( qt_atm_oi - qt_oce_ai - diag_heat - diag_trp_ei - diag_trp_es   & 
    200       !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check) 
    201          &              ) * e1e2t ) * zconv 
     199      zhfx  = glob_sum( 'icectl', ( qt_atm_oi - qt_oce_ai - diag_heat ) * e1e2t ) * zconv 
    202200 
    203201      ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  

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

    r10170 r10419  
    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, zhip_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               zhip_max(ji,jj,jl) = MAX( epsi20, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj  ,jl), h_ip_b(ji  ,jj+1,jl), & 
     104                  &                                                h_ip_b(ji-1,jj  ,jl), h_ip_b(ji  ,jj-1,jl), & 
     105                  &                                                h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), & 
     106                  &                                                h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) ) 
     107               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), & 
     108                  &                                                h_i_b (ji-1,jj  ,jl), h_i_b (ji  ,jj-1,jl), & 
     109                  &                                                h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), & 
     110                  &                                                h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) ) 
     111               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), & 
     112                  &                                                h_s_b (ji-1,jj  ,jl), h_s_b (ji  ,jj-1,jl), & 
     113                  &                                                h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), & 
     114                  &                                                h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) ) 
     115            END DO 
     116         END DO 
     117      END DO 
     118      CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. ) 
     119      ! 
     120      ! 
     121      SELECT CASE( nice_dyn )           !-- Set which dynamics is running 
     122 
     123      CASE ( np_dynALL )           !==  all dynamical processes  ==! 
     124         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology   
     125         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
     126         CALL ice_dyn_rdgrft( kt )                                                 ! -- ridging/rafting  
     127         CALL ice_cor       ( kt , 1 )                                             ! -- Corrections 
     128 
     129      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==! 
     130         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology   
     131         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
     132         CALL Hpiling                                                              ! -- simple pile-up (replaces ridging/rafting) 
     133 
     134      CASE ( np_dynADV1D )         !==  pure advection ==!   (1D) 
     135         ALLOCATE( zdivu_i(jpi,jpj) ) 
     136         ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- ! 
     137         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length 
     138         ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s  
     139         DO jj = 1, jpj 
     140            DO ji = 1, jpi 
     141               zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. ) 
     142               zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. ) 
     143               u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1) 
     144               v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1) 
     145            END DO 
     146         END DO 
     147         ! --- 
     148         CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice + correction on ice thickness 
     149 
     150         ! diagnostics: divergence at T points  
     151         DO jj = 2, jpjm1 
    97152            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( 'icedyn', 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) 
     153               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     154                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     155            END DO 
     156         END DO 
     157         CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 
     158         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) ) 
     159 
     160         DEALLOCATE( zdivu_i ) 
     161 
     162      CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities) 
     163         ALLOCATE( zdivu_i(jpi,jpj) ) 
    120164         u_ice(:,:) = rn_uice * umask(:,:,1) 
    121165         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 
     166         !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1) 
     167         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 
     168         ! --- 
     169         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
     170 
     171         ! diagnostics: divergence at T points  
     172         DO jj = 2, jpjm1 
     173            DO ji = 2, jpim1 
     174               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     175                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     176            END DO 
     177         END DO 
     178         CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 
     179         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) ) 
     180 
     181         DEALLOCATE( zdivu_i ) 
    125182 
    126183      END SELECT 
    127       ! 
    128       IF( ln_timing )   CALL timing_stop('icedyn') 
     184       ! 
     185      ! controls 
     186      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     187      IF( ln_timing    )   CALL timing_stop ('icedyn')                                                             ! timing 
    129188      ! 
    130189   END SUBROUTINE ice_dyn 
    131190 
    132191 
    133    SUBROUTINE Hbig( phmax ) 
     192   SUBROUTINE Hbig( phi_max, phs_max, phip_max ) 
    134193      !!------------------------------------------------------------------- 
    135194      !!                  ***  ROUTINE Hbig  *** 
    136195      !! 
    137196      !! ** 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) 
     197      !!              abnormally tick ice or snow 
     198      !! 
     199      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points 
     200      !!                 (before advection) and reduce it by adapting ice concentration 
     201      !!              2- check whether snow thickness is larger than the surrounding 9-points 
     202      !!                 (before advection) and reduce it by sending the excess in the ocean 
     203      !!              3- check whether snow load deplets the snow-ice interface below sea level$ 
     204      !!                 and reduce it by sending the excess in the ocean 
     205      !!              4- correct pond fraction to avoid a_ip > a_i 
    144206      !! 
    145207      !! ** input   : Max thickness of the surrounding 9-points 
    146208      !!------------------------------------------------------------------- 
    147       REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phmax   ! max ice thick from surrounding 9-pts 
     209      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    148210      ! 
    149211      INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    150       REAL(wp) ::   zh, zdv 
     212      REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra 
    151213      !!------------------------------------------------------------------- 
    152214      ! 
    … …  
    156218         DO jj = 1, jpj 
    157219            DO ji = 1, jpi 
    158                IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to hmax 
     220               IF ( v_i(ji,jj,jl) > 0._wp ) THEN 
    159221                  ! 
    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) 
     222                  !                               ! -- check h_ip -- ! 
     223                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
     224                  IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN 
     225                     zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) ) 
     226                     IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN 
     227                        a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
     228                     ENDIF 
    166229                  ENDIF 
    167230                  ! 
     231                  !                               ! -- check h_i -- ! 
     232                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
     233                  zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
     234                  IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN 
     235                     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) 
     236                  ENDIF 
     237                  ! 
     238                  !                               ! -- check h_s -- ! 
     239                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
     240                  zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl) 
     241                  IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN 
     242                     zfra = a_i(ji,jj,jl) * phs_max(ji,jj,jl) / MAX( v_s(ji,jj,jl), epsi20 ) 
     243                     ! 
     244                     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 
     245                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
     246                     ! 
     247                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra 
     248                     v_s(ji,jj,jl)          = a_i(ji,jj,jl) * phs_max(ji,jj,jl) 
     249                  ENDIF            
     250                  ! 
     251                  !                               ! -- check snow load -- ! 
     252                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 
     253                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 
     254                  !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 
     255                  zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
     256                  IF( zvs_excess > 0._wp ) THEN 
     257                     zfra = zvs_excess / MAX( v_s(ji,jj,jl), epsi20 ) 
     258                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice 
     259                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
     260                     ! 
     261                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra 
     262                     v_s(ji,jj,jl)          = v_s(ji,jj,jl) - zvs_excess 
     263                  ENDIF 
     264                   
    168265               ENDIF 
    169266            END DO 
    … …  
    215312      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
    216313      !! 
    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 
     314      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  & 
     315         &             rn_ishlat ,                                                           & 
     316         &             ln_landfast_L16, ln_landfast_home, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile 
    220317      !!------------------------------------------------------------------- 
    221318      ! 
    … …  
    233330         WRITE(numout,*) '~~~~~~~~~~~~' 
    234331         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 
     332         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynALL       = ', ln_dynALL 
     333         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV    = ', ln_dynRHGADV 
     334         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)     ln_dynADV1D     = ', ln_dynADV1D 
     335         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                 ln_dynADV2D     = ', ln_dynADV2D 
     336         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)    = (', rn_uice,',', rn_vice,')' 
     337         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat       = ', rn_ishlat 
     338         WRITE(numout,*) '      Landfast: param from Lemieux 2016                       ln_landfast_L16 = ', ln_landfast_L16 
     339         WRITE(numout,*) '      Landfast: param from home made                          ln_landfast_home= ', ln_landfast_home 
     340         WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_depfra       = ', rn_depfra 
     341         WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr       = ', rn_icebfr 
     342         WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax      = ', rn_lfrelax 
     343         WRITE(numout,*) '         isotropic tensile strength                           rn_tensile      = ', rn_tensile 
    244344         WRITE(numout,*) 
    245345      ENDIF 
    … …  
    247347      ioptio = 0  
    248348      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction) 
    249       IF( ln_dynFULL   ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynFULL      ;   ENDIF 
     349      IF( ln_dynALL    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynALL       ;   ENDIF 
    250350      !      !--- dynamics without ridging/rafting and corr   (rheology + advection) 
    251351      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 
     352      !      !--- advection 1D only - test case from Schar & Smolarkiewicz 1996 
     353      IF( ln_dynADV1D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV1D     ;   ENDIF 
     354      !      !--- advection 2D only with prescribed ice velocities (from namelist) 
     355      IF( ln_dynADV2D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV2D     ;   ENDIF 
    254356      ! 
    255357      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' ) 
    … …  
    261363      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip' 
    262364      ENDIF 
    263       !                                      !--- NO Landfast ice : set to zero once for all 
    264       IF( .NOT.ln_landfast )   tau_icebfr(:,:) = 0._wp 
     365      !                                      !--- Landfast ice 
     366      IF( .NOT.ln_landfast_L16 .AND. .NOT.ln_landfast_home )   tau_icebfr(:,:) = 0._wp 
     367      ! 
     368      IF ( ln_landfast_L16 .AND. ln_landfast_home ) THEN 
     369         CALL ctl_stop( 'ice_dyn_init: choose one and only one landfast parameterization (ln_landfast_L16 or ln_landfast_home)' ) 
     370      ENDIF 
    265371      ! 
    266372      CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters 

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

    r10069 r10419  
    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/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90

    r10397 r10419  
    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 :: amaxu, amaxv 
     38    
     39   ! advect H all the way (and get V=H*A at the end) 
     40   LOGICAL :: ll_thickness = .FALSE. 
     41    
     42   ! look for 9 points around in nonosc limiter   
     43   LOGICAL :: ll_9points = .FALSE.  ! false=better h? 
     44 
     45   ! use HgradU in nonosc limiter   
     46   LOGICAL :: ll_HgradU = .TRUE.   ! no effect? 
     47 
     48   ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme 
     49   LOGICAL :: ll_neg = .TRUE.      ! keep TRUE 
     50    
     51   ! limit the fluxes 
     52   LOGICAL :: ll_zeroup1 = .FALSE. ! false ok if Hbig otherwise needed for 2D sinon on a des valeurs de H trop fortes !! 
     53   LOGICAL :: ll_zeroup2 = .FALSE. ! false ok for 1D, 2D, 3D 
     54   LOGICAL :: ll_zeroup4 = .FALSE. ! false ok for 1D, 2D, 3D 
     55   LOGICAL :: ll_zeroup5 = .FALSE. ! false ok for 1D, 2D 
     56 
     57   ! fluxes that are limited are u*H, then (u*H)*(ua/u) is used for V (only for nonosc) 
     58   LOGICAL :: ll_clem   = .TRUE.   ! simpler than rachid and works 
     59 
     60   ! First advect H as H*=Hdiv(u), then use H* for H(n+1)=H(n)-div(uH*) 
     61   LOGICAL :: ll_gurvan = .FALSE.  ! must be false for 1D case !! 
     62 
     63   ! First guess as div(uH) (-true-) or Hdiv(u)+ugradH (-false-) 
     64   LOGICAL :: ll_1stguess_clem = .FALSE. ! better negative values but less good h 
     65 
     66   ! advect (or not) open water. If not, retrieve it from advection of A 
     67   LOGICAL :: ll_ADVopw = .FALSE.  ! 
     68    
     69   ! alternate directions for upstream 
     70   LOGICAL :: ll_upsxy = .TRUE. 
     71 
     72   ! alternate directions for high order 
     73   LOGICAL :: ll_hoxy = .TRUE. 
     74    
     75   ! prelimiter: use it to avoid overshoot in H 
     76   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? 
     77   LOGICAL :: ll_prelimiter_devore  = .FALSE.  ! from: Devore eq. 11 (worth than zalesak) 
     78 
     79   ! iterate on the limiter (only nonosc) 
     80   LOGICAL :: ll_limiter_it2 = .FALSE. 
     81    
     82 
    3683   !! * Substitutions 
    3784#  include "vectopt_loop_substitute.h90" 
    … …  
    3986   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    4087   !! $Id$ 
    41    !! Software governed by the CeCILL license (see ./LICENSE) 
     88   !! Software governed by the CeCILL licence     (./LICENSE) 
    4289   !!---------------------------------------------------------------------- 
    4390CONTAINS 
    4491 
    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 ) 
     92   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice,  & 
     93      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    4794      !!---------------------------------------------------------------------- 
    4895      !!                  ***  ROUTINE ice_dyn_adv_umx  *** 
    … …  
    54101      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    55102      !!---------------------------------------------------------------------- 
    56       INTEGER                     , INTENT(in   ) ::   k_order    ! order of the scheme (1-5 or 20) 
     103      INTEGER                     , INTENT(in   ) ::   kn_umx     ! order of the scheme (1-5=UM or 20=CEN2) 
    57104      INTEGER                     , INTENT(in   ) ::   kt         ! time step 
    58105      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity 
    … …  
    70117      ! 
    71118      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
    72       INTEGER  ::   initad                  ! number of sub-timestep for the advection 
    73       INTEGER  ::   ipl                     ! third dimention of tracer array 
    74  
    75       REAL(wp) ::   zusnit, zdt 
    76       REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! send zcflnow and receive zcflprv 
    77       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zudy, zvdx, zcu_box, zcv_box 
    78       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zpato 
     119      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
     120      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
     121      REAL(wp) ::   zdt 
     122      REAL(wp), DIMENSION(1)       ::   zcflprv, zcflnow   ! send zcflnow and receive zcflprv 
     123      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx 
     124      REAL(wp), DIMENSION(jpi,jpj) ::   zati1, zati2 
     125 
     126 
     127 
     128      REAL(wp), DIMENSION(jpi,jpj)     ::   zcu_box, zcv_box 
     129      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho 
     130      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip 
     131      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhvar 
     132 
    79133      !!---------------------------------------------------------------------- 
    80134      ! 
    81135      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 
    82136      ! 
    83       ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) ) 
    84       ALLOCATE( zpato(jpi,jpj,1) ) 
    85137      ! 
    86138      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! 
    … …  
    93145      CALL mpp_delay_max( 'icedyn_adv_umx', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 ) 
    94146 
    95       IF( zcflprv(1) > .5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp 
    96       ELSE                         ;   initad = 1   ;   zusnit = 1.0_wp 
     147      IF( zcflprv(1) > .5 ) THEN   ;   icycle = 2 
     148      ELSE                         ;   icycle = 1 
    97149      ENDIF 
    98  
    99       zdt = rdt_ice / REAL(initad) 
     150       
     151      zdt = rdt_ice / REAL(icycle) 
    100152 
    101153      ! --- transport --- ! 
    … …  
    118170      END DO 
    119171 
    120       zpato (:,:,1) = pato_i(:,:) 
    121  
     172      IF( ll_zeroup2 ) THEN 
     173         IF(.NOT. ALLOCATED(amaxu))       ALLOCATE(amaxu (jpi,jpj,jpl)) 
     174         IF(.NOT. ALLOCATED(amaxv))       ALLOCATE(amaxv (jpi,jpj,jpl)) 
     175      ENDIF 
    122176      !---------------! 
    123177      !== advection ==! 
    124178      !---------------! 
    125       DO jt = 1, initad 
     179      DO jt = 1, icycle 
     180 
     181!!$         IF( ll_ADVopw ) THEN 
     182!!$            zamsk = 1._wp 
     183!!$            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) )   ! Open water area  
     184!!$            zamsk = 0._wp 
     185!!$         ELSE 
     186            zati1(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
     187!!$         ENDIF 
    126188          
    127          l_full_nf_update = .FALSE.   ! false: disable full North fold update (performances) 
    128          CALL adv_umx( k_order, kt,   1, zdt, zudy, zvdx, zcu_box, zcv_box, zpato(:,:,1) )        ! Open water area  
    129          CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,:) )         ! Ice area 
    130          CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,:) )         ! Ice  volume 
    131          CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, psv_i(:,:,:) )        ! Salt content 
    132          CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i(:,:,:) )        ! Age content 
     189         WHERE( pa_i(:,:,:) >= epsi20 )   ;   z1_ai(:,:,:) = 1._wp / pa_i(:,:,:) 
     190         ELSEWHERE                        ;   z1_ai(:,:,:) = 0. 
     191         END WHERE 
     192            ! 
     193         WHERE( pa_ip(:,:,:) >= epsi20 )  ;   z1_aip(:,:,:) = 1._wp / pa_ip(:,:,:) 
     194         ELSEWHERE                        ;   z1_aip(:,:,:) = 0. 
     195         END WHERE 
     196         ! 
     197         IF( ll_zeroup2 ) THEN 
     198            DO jl = 1, jpl 
     199               DO jj = 2, jpjm1 
     200                  DO ji = fs_2, fs_jpim1 
     201                     amaxu(ji,jj,jl)=MAX( pa_i(ji,jj,jl), pa_i(ji,jj-1,jl), pa_i(ji,jj+1,jl), & 
     202                        &                                 pa_i(ji+1,jj,jl), pa_i(ji+1,jj-1,jl), pa_i(ji+1,jj+1,jl) ) 
     203                     amaxv(ji,jj,jl)=MAX( pa_i(ji,jj,jl), pa_i(ji-1,jj,jl), pa_i(ji+1,jj,jl), & 
     204                        &                                 pa_i(ji,jj+1,jl), pa_i(ji-1,jj+1,jl), pa_i(ji+1,jj+1,jl) ) 
     205                  END DO 
     206               END DO 
     207            END DO 
     208            CALL lbc_lnk_multi('icedyn_adv_umx', amaxu, 'T', 1., amaxv, 'T', 1.) 
     209         ENDIF 
     210         ! 
     211         DO jl = 1, jpl 
     212            zua_ho(:,:,jl) = zudy(:,:) 
     213            zva_ho(:,:,jl) = zvdx(:,:) 
     214         END DO 
     215          
     216         zamsk = 1._wp 
     217         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) ! Ice area 
     218         zamsk = 0._wp 
     219         ! 
     220         IF( ll_thickness ) THEN 
     221            DO jl = 1, jpl 
     222               zua_ho(:,:,jl) = zudy(:,:) 
     223               zva_ho(:,:,jl) = zvdx(:,:) 
     224            END DO 
     225         ENDIF 
     226            ! 
     227         zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
     228         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_i )    ! Ice volume 
     229         IF( ll_thickness )   pv_i(:,:,:) = zhvar(:,:,:) * pa_i(:,:,:) 
     230         ! 
     231         zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
     232         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_s )    ! Snw volume 
     233         IF( ll_thickness )   pv_s(:,:,:) = zhvar(:,:,:) * pa_i(:,:,:) 
     234         ! 
     235         zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 
     236         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, psv_i )    ! Salt content 
     237         ! 
     238         zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:) 
     239         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, poa_i )    ! Age content 
     240         ! 
    133241         DO jk = 1, nlay_i 
    134             CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,:) )   ! Ice  heat content 
    135          END DO 
    136          CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,:) )         ! Snow volume 
     242            zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 
     243            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) ) ! Ice heat content 
     244         END DO 
     245         ! 
    137246         DO jk = 1, nlay_s 
    138             CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,jk,:) )   ! Snow heat content 
    139          END DO 
     247            zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 
     248            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) ) ! Snw heat content 
     249         END DO 
     250            ! 
    140251         IF ( ln_pnd_H12 ) THEN 
    141             CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,:) )     ! Melt pond fraction 
    142             CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,:) )     ! Melt pond volume 
     252            ! 
     253            DO jl = 1, jpl 
     254               zua_ho(:,:,jl) = zudy(:,:) 
     255               zva_ho(:,:,jl) = zvdx(:,:) 
     256            END DO 
     257             
     258            zamsk = 1._wp 
     259            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! mp fraction 
     260            zamsk = 0._wp 
     261            ! 
     262            zhvar(:,:,:) = pv_ip(:,:,:) * z1_ai(:,:,:) 
     263            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_ip ) ! mp volume 
    143264         ENDIF 
    144  
    145          l_full_nf_update = jt == initad   ! true: enable full North fold update (performances) when exiting the loop 
    146          CALL lbc_lnk( 'icedyn_adv_umx', zpato, 'T',  1. ) 
    147          CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T',  1. ) 
    148          CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T',  1. ) 
    149          IF ( ln_pnd_H12 ) THEN 
    150             CALL lbc_lnk_multi( 'icedyn_adv_umx',  pa_i, 'T',  1., pv_i, 'T',  1., psv_i, 'T',  1., & 
    151                                                 & poa_i, 'T',  1., pv_s, 'T',  1., pa_ip, 'T',  1., & 
    152                                                 & pv_ip, 'T',  1. ) 
    153          ELSE 
    154             CALL lbc_lnk_multi( 'icedyn_adv_umx',  pa_i, 'T',  1., pv_i, 'T',  1., psv_i, 'T',  1., & 
    155                                                 & poa_i, 'T',  1., pv_s, 'T',  1. ) 
    156          ENDIF 
    157           
     265         ! 
     266         ! 
     267!!$         IF( .NOT. ll_ADVopw ) THEN 
     268            zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
     269            DO jj = 2, jpjm1 
     270               DO ji = fs_2, fs_jpim1 
     271                  pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                                                  ! Open water area 
     272                     &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt 
     273               END DO 
     274            END DO 
     275            CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T',  1. ) 
     276!!$         ENDIF 
     277         ! 
    158278      END DO 
    159279      ! 
    160       pato_i(:,:) = zpato (:,:,1) 
    161       ! 
    162       DEALLOCATE( zudy, zvdx, zcu_box, zcv_box, zpato ) 
    163       ! 
    164280   END SUBROUTINE ice_dyn_adv_umx 
    165281 
    166282    
    167    SUBROUTINE adv_umx( k_order, kt, ipl, pdt, puc, pvc, pubox, pvbox, ptc ) 
     283   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho ) 
    168284      !!---------------------------------------------------------------------- 
    169285      !!                  ***  ROUTINE adv_umx  *** 
    … …  
    178294      !! ** Action : - pt  the after advective tracer 
    179295      !!---------------------------------------------------------------------- 
    180       INTEGER                         , INTENT(in   ) ::   k_order        ! order of the ULTIMATE scheme 
    181       INTEGER                         , INTENT(in   ) ::   kt             ! number of iteration 
    182       INTEGER                         , INTENT(in   ) ::   ipl            ! third dimension of tracer array 
    183       REAL(wp)                        , INTENT(in   ) ::   pdt            ! tracer time-step 
    184       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   puc  , pvc     ! 2 ice velocity components => u*e2 
    185       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pubox, pvbox   ! upstream velocity 
    186       REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(inout) ::   ptc            ! tracer content field 
     296      REAL(wp)                    , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0) 
     297      INTEGER                     , INTENT(in   )           ::   kn_umx         ! order of the scheme (1-5=UM or 20=CEN2) 
     298      INTEGER                     , INTENT(in   )           ::   jt             ! number of sub-iteration 
     299      INTEGER                     , INTENT(in   )           ::   kt             ! number of iteration 
     300      REAL(wp)                    , INTENT(in   )           ::   pdt            ! tracer time-step 
     301      REAL(wp), DIMENSION(:,:  ), INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2 
     302      REAL(wp), DIMENSION(:,:,:), INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u 
     303      REAL(wp), DIMENSION(:,:  ), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity 
     304      REAL(wp), DIMENSION(:,:,:), INTENT(inout)           ::   pt             ! tracer field 
     305      REAL(wp), DIMENSION(:,:,:), INTENT(inout)           ::   ptc            ! tracer content field 
     306      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes 
    187307      ! 
    188308      INTEGER  ::   ji, jj, jl       ! dummy loop indices   
    189  
    190309      REAL(wp) ::   ztra             ! local scalar 
    191       REAL(wp), DIMENSION(jpi,jpj,ipl) ::   zfu_ups, zfu_ho, zt_u, zt_ups 
    192       REAL(wp), DIMENSION(jpi,jpj,ipl) ::   zfv_ups, zfv_ho, zt_v, ztrd 
    193  
    194       DO jl = 1, ipl 
    195       !!---------------------------------------------------------------------- 
    196       ! 
    197       !  upstream advection with initial mass fluxes & intermediate update 
    198       ! -------------------------------------------------------------------- 
    199       DO jj = 1, jpjm1         ! upstream tracer flux in the i and j direction 
    200          DO ji = 1, fs_jpim1   ! vector opt. 
    201             zfu_ups(ji,jj,jl) = MAX( puc(ji,jj), 0._wp ) * ptc(ji,jj,jl) + MIN( puc(ji,jj), 0._wp ) * ptc(ji+1,jj,jl) 
    202             zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj), 0._wp ) * ptc(ji,jj,jl) + MIN( pvc(ji,jj), 0._wp ) * ptc(ji,jj+1,jl) 
     310      INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid) 
     311      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt 
     312      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ups, zfv_ups, zt_ups   ! only for nonosc  
     313      !!---------------------------------------------------------------------- 
     314      ! 
     315      !  upstream (_ups) advection with initial mass fluxes 
     316      ! --------------------------------------------------- 
     317 
     318      IF( ll_gurvan .AND. pamsk==0. ) THEN 
     319         DO jl = 1, jpl 
     320            DO jj = 2, jpjm1 
     321               DO ji = fs_2, fs_jpim1 
     322                  pt(ji,jj,jl) = ( pt (ji,jj,jl) + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj)     & 
     323                     &                           + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) )   & 
     324                     &           * tmask(ji,jj,1) 
     325               END DO 
     326            END DO 
     327         END DO 
     328         CALL lbc_lnk( 'icedyn_adv_umx', pt, 'T', 1. ) 
     329      ENDIF 
     330 
     331       
     332      IF( .NOT. ll_upsxy ) THEN 
     333 
     334         ! fluxes in both x-y directions 
     335         DO jl = 1, jpl 
     336            DO jj = 1, jpjm1 
     337               DO ji = 1, fs_jpim1 
     338                  IF( ll_clem ) THEN 
     339                     zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
     340                     zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
     341                  ELSE 
     342                     zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * pt(ji+1,jj,jl) 
     343                     zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj+1,jl) 
     344                  ENDIF 
     345               END DO 
     346            END DO 
     347         END DO 
     348 
     349      ELSE 
     350         ! 
     351         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     352            ! flux in x-direction 
     353            DO jl = 1, jpl 
     354               DO jj = 1, jpjm1 
     355                  DO ji = 1, fs_jpim1 
     356                     IF( ll_clem ) THEN 
     357                        zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
     358                     ELSE 
     359                        zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * pt(ji+1,jj,jl) 
     360                     ENDIF 
     361                  END DO 
     362               END DO 
     363            END DO 
     364             
     365            ! first guess of tracer content from u-flux 
     366            DO jl = 1, jpl 
     367               DO jj = 2, jpjm1 
     368                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     369                     IF( ll_clem ) THEN 
     370                        IF( ll_gurvan ) THEN 
     371                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     372                        ELSE 
     373                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     374                              &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     375                              &            ) * tmask(ji,jj,1) 
     376                        ENDIF 
     377                     ELSE 
     378                        zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) & 
     379                           &         * tmask(ji,jj,1) 
     380                     ENDIF 
     381                     !!                  IF( ji==26 .AND. jj==86) THEN 
     382                     !!                     WRITE(numout,*) '************************' 
     383                     !!                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 
     384                     !!                  ENDIF 
     385                  END DO 
     386               END DO 
     387            END DO 
     388            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
     389            ! 
     390            ! flux in y-direction 
     391            DO jl = 1, jpl 
     392               DO jj = 1, jpjm1 
     393                  DO ji = 1, fs_jpim1 
     394                     IF( ll_clem ) THEN 
     395                        zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 
     396                     ELSE 
     397                        zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * zpt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * zpt(ji,jj+1,jl) 
     398                     ENDIF 
     399                  END DO 
     400               END DO 
     401            END DO 
     402         ! 
     403         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     404            ! flux in y-direction 
     405            DO jl = 1, jpl 
     406               DO jj = 1, jpjm1 
     407                  DO ji = 1, fs_jpim1 
     408                     IF( ll_clem ) THEN 
     409                        zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
     410                     ELSE 
     411                        zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj+1,jl) 
     412                     ENDIF 
     413                  END DO 
     414               END DO 
     415            END DO 
     416 
     417            ! first guess of tracer content from v-flux 
     418            DO jl = 1, jpl 
     419               DO jj = 2, jpjm1 
     420                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     421                     IF( ll_clem ) THEN 
     422                        IF( ll_gurvan ) THEN 
     423                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     424                        ELSE 
     425                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 
     426                              &            + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
     427                              &            * tmask(ji,jj,1) 
     428                        ENDIF 
     429                     ELSE 
     430                        zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) & 
     431                           &            * tmask(ji,jj,1) 
     432                     ENDIF 
     433                     !!                  IF( ji==26 .AND. jj==86) THEN 
     434                     !!                     WRITE(numout,*) '************************' 
     435                     !!                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 
     436                     !!                  ENDIF 
     437                  END DO 
     438               END DO 
     439            END DO 
     440            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
     441            ! 
     442            ! flux in x-direction 
     443            DO jl = 1, jpl 
     444               DO jj = 1, jpjm1 
     445                  DO ji = 1, fs_jpim1 
     446                     IF( ll_clem ) THEN 
     447                        zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 
     448                     ELSE 
     449                        zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * zpt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * zpt(ji+1,jj,jl) 
     450                     ENDIF 
     451                  END DO 
     452               END DO 
     453            END DO 
     454            ! 
     455         ENDIF 
     456          
     457      ENDIF 
     458 
     459      IF( ll_clem .AND. kn_limiter /= 1 )  & 
     460         & CALL ctl_stop( 'STOP', 'icedyn_adv_umx: ll_clem incompatible with limiters other than nonosc' ) 
     461 
     462      IF( ll_zeroup2 ) THEN 
     463         DO jl = 1, jpl 
     464            DO jj = 1, jpjm1 
     465               DO ji = 1, fs_jpim1   ! vector opt. 
     466                  IF( amaxu(ji,jj,jl) == 0._wp )   zfu_ups(ji,jj,jl) = 0._wp 
     467                  IF( amaxv(ji,jj,jl) == 0._wp )   zfv_ups(ji,jj,jl) = 0._wp 
     468               END DO 
     469            END DO 
     470         END DO 
     471      ENDIF 
     472 
     473      ! guess after content field with upstream scheme 
     474      DO jl = 1, jpl 
     475         DO jj = 2, jpjm1 
     476            DO ji = fs_2, fs_jpim1 
     477               ztra          = - (   zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj  ,jl)   & 
     478                  &                + zfv_ups(ji,jj,jl) - zfv_ups(ji  ,jj-1,jl) ) * r1_e1e2t(ji,jj) 
     479               IF( ll_clem ) THEN 
     480                  IF( ll_gurvan ) THEN 
     481                     zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 
     482                  ELSE 
     483                     zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra + ( pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) )   & 
     484                        &                                            +   pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) ) & 
     485                        &                                              * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     486                  ENDIF 
     487               ELSE 
     488                  zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 
     489               ENDIF 
     490               !!            IF( ji==26 .AND. jj==86) THEN 
     491               !!               WRITE(numout,*) '**************************' 
     492               !!               WRITE(numout,*) 'zt upstream',zt_ups(ji,jj) 
     493               !!            ENDIF 
     494            END DO 
    203495         END DO 
    204496      END DO 
    205        
    206          DO jj = 2, jpjm1            ! total intermediate advective trends 
    207             DO ji = fs_2, fs_jpim1   ! vector opt. 
    208                ztra = - (   zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj  ,jl)   & 
    209                   &       + zfv_ups(ji,jj,jl) - zfv_ups(ji  ,jj-1,jl)   ) * r1_e1e2t(ji,jj) 
    210             ! 
    211                ztrd(ji,jj,jl) =                         ztra                         ! upstream trend [ -div(uh) or -div(uhT) ]   
    212                zt_ups (ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)   ! guess after content field with monotonic scheme 
    213             END DO 
    214          END DO 
    215       END DO 
    216        
     497      CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. ) 
     498 
    217499      ! High order (_ho) fluxes  
    218500      ! ----------------------- 
    219       SELECT CASE( k_order ) 
    220       CASE ( 20 )                          ! centered second order 
    221          DO jl = 1, ipl 
     501      SELECT CASE( kn_umx ) 
     502         ! 
     503      CASE ( 20 )                          !== centered second order ==! 
     504         ! 
     505         CALL cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  & 
     506            &       zt_ups, zfu_ups, zfv_ups ) 
     507         ! 
     508      CASE ( 1:5 )                         !== 1st to 5th order ULTIMATE-MACHO scheme ==! 
     509         ! 
     510         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,  & 
     511            &        zt_ups, zfu_ups, zfv_ups ) 
     512         ! 
     513      END SELECT 
     514 
     515      IF( ll_thickness ) THEN 
     516         ! final trend with corrected fluxes 
     517         ! ------------------------------------ 
     518         DO jl = 1, jpl 
     519            DO jj = 2, jpjm1 
     520               DO ji = fs_2, fs_jpim1 
     521                  IF( ll_gurvan ) THEN 
     522                     ztra       = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj)  
     523                  ELSE 
     524                     ztra       = ( - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )  &  
     525                        &           + ( pt(ji,jj,jl) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) & 
     526                        &           + ( pt(ji,jj,jl) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj) 
     527                  ENDIF 
     528                  pt(ji,jj,jl) = ( pt(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 
     529                   
     530                  IF( pt(ji,jj,jl) < -epsi20 ) THEN 
     531                     WRITE(numout,*) 'T<0 ',pt(ji,jj,jl) 
     532                  ENDIF 
     533                   
     534                  IF( pt(ji,jj,jl) < 0._wp .AND. pt(ji,jj,jl) >= -epsi20 )   pt(ji,jj,jl) = 0._wp 
     535                   
     536                  !!               IF( ji==26 .AND. jj==86) THEN 
     537                  !!                  WRITE(numout,*) 'zt high order',pt(ji,jj) 
     538                  !!               ENDIF 
     539               END DO 
     540            END DO 
     541         END DO 
     542         CALL lbc_lnk( 'icedyn_adv_umx', pt, 'T',  1. ) 
     543      ENDIF 
     544    
     545      ! Rachid trick 
     546      ! ------------ 
     547      IF( ll_clem ) THEN 
     548         IF( pamsk == 0. ) THEN 
     549            DO jl = 1, jpl 
     550               DO jj = 1, jpjm1 
     551                  DO ji = 1, fs_jpim1 
     552                     IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN 
     553                        zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 
     554                        zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 
     555                     ELSE 
     556                        zfu_ho (ji,jj,jl) = 0._wp 
     557                        zfu_ups(ji,jj,jl) = 0._wp 
     558                     ENDIF 
     559                     ! 
     560                     IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN 
     561                        zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 
     562                        zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 
     563                     ELSE 
     564                        zfv_ho (ji,jj,jl) = 0._wp   
     565                        zfv_ups(ji,jj,jl) = 0._wp   
     566                     ENDIF 
     567                  END DO 
     568               END DO 
     569            END DO 
     570         ENDIF 
     571      ENDIF 
     572 
     573      IF( ll_zeroup5 ) THEN 
     574         DO jl = 1, jpl 
     575            DO jj = 2, jpjm1 
     576               DO ji = 2, fs_jpim1   ! vector opt. 
     577                  zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     578                     &                            - ( zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  ) * tmask(ji,jj,1) 
     579                  IF( zpt(ji,jj,jl) < 0. ) THEN 
     580                     zfu_ho(ji  ,jj,jl) = zfu_ups(ji  ,jj,jl) 
     581                     zfu_ho(ji-1,jj,jl) = zfu_ups(ji-1,jj,jl) 
     582                     zfv_ho(ji  ,jj,jl) = zfv_ups(ji  ,jj,jl) 
     583                     zfv_ho(ji,jj-1,jl) = zfv_ups(ji,jj-1,jl) 
     584                  ENDIF 
     585               END DO 
     586            END DO 
     587         END DO 
     588         CALL lbc_lnk_multi( 'icedyn_adv_umx', zfu_ho, 'U',  -1., zfv_ho, 'V',  -1. ) 
     589      ENDIF 
     590 
     591      ! output high order fluxes u*a 
     592      ! ---------------------------- 
     593      IF( PRESENT( pua_ho ) ) THEN 
     594         DO jl = 1, jpl 
    222595            DO jj = 1, jpjm1 
    223                DO ji = 1, fs_jpim1   ! vector opt. 
    224                   zfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj) * ( ptc(ji,jj,jl) + ptc(ji+1,jj,jl) ) 
    225                   zfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj,jl) + ptc(ji,jj+1,jl) ) 
    226                END DO 
    227             END DO 
    228          END DO 
    229          ! 
    230       CASE ( 1:5 )                      ! 1st to 5th order ULTIMATE-MACHO scheme 
    231          CALL macho( k_order, kt, ipl, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v ) 
    232          ! 
    233          DO jl = 1, ipl 
     596               DO ji = 1, fs_jpim1 
     597                  pua_ho(ji,jj,jl) = zfu_ho(ji,jj,jl) 
     598                  pva_ho(ji,jj,jl) = zfv_ho(ji,jj,jl) 
     599               END DO 
     600            END DO 
     601         END DO 
     602      ENDIF 
     603 
     604 
     605      IF( .NOT.ll_thickness ) THEN 
     606         ! final trend with corrected fluxes 
     607         ! ------------------------------------ 
     608         DO jl = 1, jpl 
    234609            DO jj = 2, jpjm1 
    235                DO ji = 1, fs_jpim1   ! vector opt. 
    236                   zfu_ho(ji,jj,jl) = puc(ji,jj) * zt_u(ji,jj,jl) 
    237                END DO 
    238             END DO 
    239          END DO 
    240          DO jl = 1, ipl 
     610               DO ji = fs_2, fs_jpim1  
     611                  ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) * pdt   
     612                   
     613                  ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra ) * tmask(ji,jj,1) 
     614                   
     615                  !!               IF( ji==26 .AND. jj==86) THEN 
     616                  !!                  WRITE(numout,*) 'ztc high order',ptc(ji,jj) 
     617                  !!               ENDIF 
     618                   
     619               END DO 
     620            END DO 
     621         END DO 
     622         CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1. ) 
     623      ENDIF 
     624       
     625      ! 
     626   END SUBROUTINE adv_umx 
     627 
     628   SUBROUTINE cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, & 
     629      &             pt_ups, pfu_ups, pfv_ups ) 
     630      !!--------------------------------------------------------------------- 
     631      !!                    ***  ROUTINE macho  *** 
     632      !!      
     633      !! **  Purpose :   compute   
     634      !! 
     635      !! **  Method  :   ... ??? 
     636      !!                 TIM = transient interpolation Modeling  
     637      !! 
     638      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
     639      !!---------------------------------------------------------------------- 
     640      REAL(wp)                    , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     641      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter 
     642      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration 
     643      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration 
     644      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step 
     645      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt               ! tracer fields 
     646      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components 
     647      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components 
     648      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptc              ! tracer content at before time step  
     649      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes  
     650      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content  
     651      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes  
     652      ! 
     653      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     654      LOGICAL  ::   ll_xy = .TRUE. 
     655      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zzt 
     656      !!---------------------------------------------------------------------- 
     657      ! 
     658      IF( .NOT.ll_xy ) THEN   !-- no alternate directions --! 
     659         ! 
     660         DO jl = 1, jpl 
    241661            DO jj = 1, jpjm1 
    242                DO ji = fs_2, fs_jpim1   ! vector opt. 
    243                   zfv_ho(ji,jj,jl) = pvc(ji,jj) * zt_v(ji,jj,jl) 
    244                END DO 
    245             END DO 
    246          END DO 
    247          ! 
    248       END SELECT 
     662               DO ji = 1, fs_jpim1 
     663                  IF( ll_clem ) THEN 
     664                     pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
     665                     pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
     666                  ELSE 
     667                     pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
     668                     pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
     669                  ENDIF 
     670               END DO 
     671            END DO 
     672         END DO 
     673         IF    ( kn_limiter == 1 ) THEN 
     674            IF( ll_clem ) THEN 
     675               CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     676            ELSE 
     677               CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     678            ENDIF 
     679         ELSEIF( kn_limiter == 2 ) THEN 
     680            CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     681            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     682         ELSEIF( kn_limiter == 3 ) THEN 
     683            CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     684            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     685         ENDIF 
     686         ! 
     687      ELSE                    !-- alternate directions --! 
     688         ! 
     689         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     690            ! 
     691            ! flux in x-direction 
     692            DO jl = 1, jpl 
     693               DO jj = 1, jpjm1 
     694                  DO ji = 1, fs_jpim1 
     695                     IF( ll_clem ) THEN 
     696                        pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
     697                     ELSE 
     698                        pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
     699                     ENDIF 
     700                  END DO 
     701               END DO 
     702            END DO 
     703            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     704            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     705 
     706            ! first guess of tracer content from u-flux 
     707            DO jl = 1, jpl 
     708               DO jj = 2, jpjm1 
     709                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     710                     IF( ll_clem ) THEN 
     711                        zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)        & 
     712                           &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
     713                           &         * tmask(ji,jj,1) 
     714                     ELSE                      
     715                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     716                     ENDIF 
     717                  END DO 
     718               END DO 
     719            END DO 
     720            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     721 
     722            ! flux in y-direction 
     723            DO jl = 1, jpl 
     724               DO jj = 1, jpjm1 
     725                  DO ji = 1, fs_jpim1 
     726                     IF( ll_clem ) THEN 
     727                        pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) ) 
     728                     ELSE                      
     729                        pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) ) 
     730                     ENDIF 
     731                  END DO 
     732               END DO 
     733            END DO 
     734            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     735            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     736 
     737         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     738            ! 
     739            ! flux in y-direction 
     740            DO jl = 1, jpl 
     741               DO jj = 1, jpjm1 
     742                  DO ji = 1, fs_jpim1 
     743                     IF( ll_clem ) THEN 
     744                        pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
     745                     ELSE                      
     746                        pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
     747                     ENDIF 
     748                  END DO 
     749               END DO 
     750            END DO 
     751            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     752            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     753            ! 
     754            ! first guess of tracer content from v-flux 
     755            DO jl = 1, jpl 
     756               DO jj = 2, jpjm1 
     757                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     758                     IF( ll_clem ) THEN 
     759                        zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 
     760                           &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
     761                           &         * tmask(ji,jj,1) 
     762                     ELSE 
     763                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     764                     ENDIF 
     765                  END DO 
     766               END DO 
     767            END DO 
     768            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     769            ! 
     770            ! flux in x-direction 
     771            DO jl = 1, jpl 
     772               DO jj = 1, jpjm1 
     773                  DO ji = 1, fs_jpim1 
     774                     IF( ll_clem ) THEN 
     775                        pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) ) 
     776                     ELSE 
     777                        pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) ) 
     778                     ENDIF 
     779                  END DO 
     780               END DO 
     781            END DO 
     782            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     783            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     784 
     785         ENDIF 
     786         IF( ll_clem ) THEN 
     787            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 ) 
     788         ELSE 
     789            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 ) 
     790         ENDIF 
    249791          
    250       ! antidiffusive flux : high order minus low order 
    251       ! -------------------------------------------------- 
    252       DO jl = 1, ipl 
    253          DO jj = 2, jpjm1 
    254             DO ji = 1, fs_jpim1   ! vector opt. 
    255                zfu_ho(ji,jj,jl) = zfu_ho(ji,jj,jl) - zfu_ups(ji,jj,jl) 
    256             END DO 
    257          END DO 
    258       END DO 
    259       DO jl = 1, ipl 
    260          DO jj = 1, jpjm1 
    261             DO ji = fs_2, fs_jpim1   ! vector opt. 
    262                zfv_ho(ji,jj,jl) = zfv_ho(ji,jj,jl) - zfv_ups(ji,jj,jl) 
    263             END DO 
    264          END DO 
    265       END DO 
    266  
    267       CALL lbc_lnk("icedyn_adv_umx",zt_ups, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign) 
    268        
    269       ! monotonicity algorithm 
    270       ! ------------------------- 
    271       CALL nonosc_2d( ipl, ptc, zfu_ho, zfv_ho, zt_ups, pdt ) 
    272        
    273       ! final trend with corrected fluxes 
    274       ! ------------------------------------ 
    275       DO jl = 1, ipl 
    276          DO jj = 2, jpjm1 
    277             DO ji = fs_2, fs_jpim1   ! vector opt.   
    278                ztra       = ztrd(ji,jj,jl)  - (  zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj  ,jl)   & 
    279                   &                         +    zfv_ho(ji,jj,jl) - zfv_ho(ji  ,jj-1,jl) ) * r1_e1e2t(ji,jj)   
    280                ptc(ji,jj,jl) = ptc(ji,jj,jl) + pdt * ztra * tmask(ji,jj,1) 
    281             END DO 
    282          END DO 
    283       END DO 
    284       ! 
    285    END SUBROUTINE adv_umx 
    286  
    287    SUBROUTINE macho( k_order, kt, ipl, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 
     792      ENDIF 
     793    
     794   END SUBROUTINE cen2 
     795 
     796    
     797   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, & 
     798      &              pt_ups, pfu_ups, pfv_ups ) 
     799      !!--------------------------------------------------------------------- 
     800      !!                    ***  ROUTINE macho  *** 
     801      !!      
     802      !! **  Purpose :   compute   
     803      !! 
     804      !! **  Method  :   ... ??? 
     805      !!                 TIM = transient interpolation Modeling  
     806      !! 
     807      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
     808      !!---------------------------------------------------------------------- 
     809      REAL(wp)                    , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     810      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter 
     811      INTEGER                     , INTENT(in   ) ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2) 
     812      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration 
     813      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration 
     814      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step 
     815      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt               ! tracer fields 
     816      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components 
     817      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components 
     818      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pubox, pvbox     ! upstream velocity 
     819      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptc              ! tracer content at before time step  
     820      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_u, pt_v       ! tracer at u- and v-points  
     821      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes  
     822      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content  
     823      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes  
     824      ! 
     825      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     826      REAL(wp) ::   ztra 
     827      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zzt, zzfu_ho, zzfv_ho 
     828      !!---------------------------------------------------------------------- 
     829      ! 
     830      IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     831         ! 
     832         !                                                        !--  ultimate interpolation of pt at u-point  --! 
     833         CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 
     834         !                                                        !--  limiter in x --! 
     835         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     836         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     837         !                                                        !--  advective form update in zzt  --! 
     838 
     839         IF( ll_1stguess_clem ) THEN 
     840 
     841            ! first guess of tracer content from u-flux 
     842            DO jl = 1, jpl 
     843               DO jj = 2, jpjm1 
     844                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     845                     IF( ll_clem ) THEN 
     846                        IF( ll_gurvan ) THEN 
     847                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     848                        ELSE 
     849                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     850                              &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     851                              &         ) * tmask(ji,jj,1) 
     852                        ENDIF 
     853                     ELSE 
     854                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     855                     ENDIF 
     856                  END DO 
     857               END DO 
     858            END DO 
     859            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     860 
     861         ELSE 
     862 
     863            DO jl = 1, jpl 
     864               DO jj = 2, jpjm1 
     865                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     866                     IF( ll_gurvan ) THEN 
     867                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj   ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj)  & 
     868                           &                         - pt   (ji,jj,jl) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) 
     869                     ELSE 
     870                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj   ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj)  & 
     871                           &                         - pt   (ji,jj,jl) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk 
     872                     ENDIF 
     873                     zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 
     874                  END DO 
     875               END DO 
     876            END DO 
     877            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     878         ENDIF 
     879         ! 
     880         !                                                        !--  ultimate interpolation of pt at v-point  --! 
     881         IF( ll_hoxy ) THEN 
     882            CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho ) 
     883         ELSE 
     884            CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 
     885         ENDIF 
     886         !                                                        !--  limiter in y --! 
     887         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     888         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     889         !          
     890         ! 
     891      ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     892         ! 
     893         !                                                        !--  ultimate interpolation of pt at v-point  --! 
     894         CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 
     895         !                                                        !--  limiter in y --! 
     896         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     897         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     898         !                                                        !--  advective form update in zzt  --! 
     899         IF( ll_1stguess_clem ) THEN 
     900             
     901            ! first guess of tracer content from v-flux  
     902            DO jl = 1, jpl 
     903               DO jj = 2, jpjm1 
     904                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     905                     IF( ll_clem ) THEN 
     906                        IF( ll_gurvan ) THEN 
     907                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     908                        ELSE 
     909                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 
     910                              &            + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     911                              &         ) * tmask(ji,jj,1) 
     912                        ENDIF 
     913                     ELSE 
     914                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) & 
     915                           &         * tmask(ji,jj,1) 
     916                     ENDIF 
     917                  END DO 
     918               END DO 
     919            END DO 
     920            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     921             
     922         ELSE 
     923             
     924            DO jl = 1, jpl 
     925               DO jj = 2, jpjm1 
     926                  DO ji = fs_2, fs_jpim1 
     927                     IF( ll_gurvan ) THEN 
     928                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj   ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj)  & 
     929                           &                         - pt   (ji,jj,jl) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) 
     930                     ELSE 
     931                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj   ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj)  & 
     932                           &                         - pt   (ji,jj,jl) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk 
     933                     ENDIF 
     934                     zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 
     935                  END DO 
     936               END DO 
     937            END DO 
     938            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     939         ENDIF 
     940         ! 
     941         !                                                        !--  ultimate interpolation of pt at u-point  --! 
     942         IF( ll_hoxy ) THEN 
     943            CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho ) 
     944         ELSE 
     945            CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 
     946         ENDIF 
     947         !                                                        !--  limiter in x --! 
     948         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     949         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     950         ! 
     951         ! 
     952      ENDIF 
     953 
     954      
     955      IF( kn_limiter == 1 ) THEN 
     956         IF( .NOT. ll_limiter_it2 ) THEN 
     957            IF( ll_clem ) THEN 
     958               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     959            ELSE 
     960               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     961            ENDIF 
     962         ELSE 
     963            zzfu_ho(:,:,:) = pfu_ho(:,:,:) 
     964            zzfv_ho(:,:,:) = pfv_ho(:,:,:) 
     965            ! 1st iteration of nonosc (limit the flux with the upstream solution) 
     966            IF( ll_clem ) THEN 
     967               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 
     968            ELSE 
     969               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 
     970            ENDIF 
     971            ! guess after content field with high order 
     972            DO jl = 1, jpl 
     973               DO jj = 2, jpjm1 
     974                  DO ji = fs_2, fs_jpim1 
     975                     ztra          = - ( zzfu_ho(ji,jj,jl) - zzfu_ho(ji-1,jj,jl) + zzfv_ho(ji,jj,jl) - zzfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) 
     976                     zzt(ji,jj,jl) =   ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 
     977                  END DO 
     978               END DO 
     979            END DO 
     980            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     981            ! 2nd iteration of nonosc (limit the flux with the limited high order solution) 
     982            IF( ll_clem ) THEN 
     983               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 
     984            ELSE 
     985               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 
     986            ENDIF 
     987         ENDIF 
     988      ENDIF 
     989      ! 
     990   END SUBROUTINE macho 
     991 
     992 
     993   SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 
    288994      !!--------------------------------------------------------------------- 
    289995      !!                    ***  ROUTINE ultimate_x  *** 
    … …  
    2961002      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    2971003      !!---------------------------------------------------------------------- 
    298       INTEGER                         , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme 
    299       INTEGER                         , INTENT(in   ) ::   kt         ! number of iteration 
    300       INTEGER                         , INTENT(in   ) ::   ipl        ! third dimension of tracer array 
    301       REAL(wp)                        , INTENT(in   ) ::   pdt        ! tracer time-step 
    302       REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(in   ) ::   ptc        ! tracer fields 
    303       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   puc, pvc   ! 2 ice velocity components 
    304       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pubox, pvbox   ! upstream velocity 
    305       REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(  out) ::   pt_u, pt_v ! tracer at u- and v-points  
    306       ! 
    307       INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    308       REAL(wp) ::   zc_box    !   -      - 
    309       REAL(wp), DIMENSION(jpi,jpj,ipl) :: zzt 
    310       !!---------------------------------------------------------------------- 
    311       ! 
    312       IF( MOD( (kt - 1) / nn_fsbc , 2 ) == 0 ) THEN         !==  odd ice time step:  adv_x then adv_y  ==! 
    313          ! 
    314          !                                                           !--  ultimate interpolation of pt at u-point  --! 
    315          CALL ultimate_x( k_order, ipl, pdt, ptc, puc, pt_u ) 
    316          ! 
    317          !                                                           !--  advective form update in zzt  --! 
    318          DO jl = 1, ipl 
    319             DO jj = 2, jpjm1 
    320                DO ji = fs_2, fs_jpim1   ! vector opt. 
    321                   zzt(ji,jj,jl) = ptc(ji,jj,jl) - pubox(ji,jj) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj)  & 
    322                    &            - ptc(ji,jj,jl) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj) 
    323                   zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 
    324                END DO 
    325             END DO 
    326          END DO 
    327          CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
    328          ! 
    329          !                                                           !--  ultimate interpolation of pt at v-point  --! 
    330          CALL ultimate_y( k_order, ipl, pdt, zzt, pvc, pt_v ) 
    331          ! 
    332       ELSE                                                  !==  even ice time step:  adv_y then adv_x  ==! 
    333          ! 
    334          !                                                           !--  ultimate interpolation of pt at v-point  --! 
    335          CALL ultimate_y( k_order, ipl, pdt, ptc, pvc, pt_v ) 
    336          ! 
    337          !                                                           !--  advective form update in zzt  --! 
    338          DO jl = 1, ipl 
    339             DO jj = 2, jpjm1 
    340                DO ji = fs_2, fs_jpim1 
    341                   zzt(ji,jj,jl) = ptc(ji,jj,jl) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj)  & 
    342                      &                    - ptc  (ji,jj,jl) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj) 
    343                   zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 
    344                END DO 
    345             END DO 
    346          END DO 
    347          CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
    348          ! 
    349          !                                                           !--  ultimate interpolation of pt at u-point  --! 
    350          CALL ultimate_x( k_order, ipl, pdt, zzt, puc, pt_u ) 
    351          !       
    352       ENDIF       
    353       ! 
    354    END SUBROUTINE macho 
    355  
    356  
    357    SUBROUTINE ultimate_x( k_order, ipl, pdt, pt, puc, pt_u ) 
    358       !!--------------------------------------------------------------------- 
    359       !!                    ***  ROUTINE ultimate_x  *** 
    360       !!      
    361       !! **  Purpose :   compute   
    362       !! 
    363       !! **  Method  :   ... ??? 
    364       !!                 TIM = transient interpolation Modeling  
    365       !! 
    366       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    367       !!---------------------------------------------------------------------- 
    368       INTEGER                         , INTENT(in   ) ::   k_order   ! ocean time-step index 
    369       INTEGER                         , INTENT(in   ) ::   ipl        ! third dimension of tracer array 
    370       REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step 
    371       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   puc       ! ice i-velocity component 
    372       REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(in   ) ::   pt        ! tracer fields 
    373       REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(  out) ::   pt_u      ! tracer at u-point  
    374       ! 
    375       INTEGER  ::   ji, jj, jl       ! dummy loop indices 
     1004      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
     1005      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step 
     1006      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu        ! ice i-velocity component 
     1007      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   puc       ! ice i-velocity * A component 
     1008      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt        ! tracer fields 
     1009      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_u      ! tracer at u-point  
     1010      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho    ! high order flux  
     1011      ! 
     1012      INTEGER  ::   ji, jj, jl             ! dummy loop indices 
    3761013      REAL(wp) ::   zcu, zdx2, zdx4    !   -      - 
    377       REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu3 
    378       REAL(wp), DIMENSION(jpi,jpj,ipl) :: ztu2, ztu4 
     1014      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztu1, ztu2, ztu3, ztu4 
    3791015      !!---------------------------------------------------------------------- 
    3801016      ! 
    3811017      !                                                     !--  Laplacian in i-direction  --! 
    382       DO jl = 1, ipl 
     1018      DO jl = 1, jpl 
    3831019         DO jj = 2, jpjm1         ! First derivative (gradient) 
    3841020            DO ji = 1, fs_jpim1 
    385                ztu1(ji,jj) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
     1021               ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
    3861022            END DO 
    3871023            !                     ! Second derivative (Laplacian) 
    3881024            DO ji = fs_2, fs_jpim1 
    389                ztu2(ji,jj,jl) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj) 
     1025               ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj) 
    3901026            END DO 
    3911027         END DO 
    … …  
    3941030      ! 
    3951031      !                                                     !--  BiLaplacian in i-direction  --! 
    396       DO jl = 1, ipl 
     1032      DO jl = 1, jpl 
    3971033         DO jj = 2, jpjm1         ! Third derivative 
    3981034            DO ji = 1, fs_jpim1 
    399                ztu3(ji,jj) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
     1035               ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
    4001036            END DO 
    4011037            !                     ! Fourth derivative 
    4021038            DO ji = fs_2, fs_jpim1 
    403                ztu4(ji,jj,jl) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj) 
     1039               ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj) 
    4041040            END DO 
    4051041         END DO 
    … …  
    4081044      ! 
    4091045      ! 
    410       SELECT CASE (k_order ) 
     1046      SELECT CASE (kn_umx ) 
    4111047      ! 
    4121048      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21) 
    4131049         !         
    414          DO jl = 1, ipl 
    415             DO jj = 2, jpjm1 
     1050         DO jl = 1, jpl 
     1051            DO jj = 1, jpjm1 
    4161052               DO ji = 1, fs_jpim1   ! vector opt. 
    417                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    418                      &                                       - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
     1053                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                           pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     1054                     &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
    4191055               END DO 
    4201056            END DO 
    … …  
    4231059      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23) 
    4241060         ! 
    425          DO jl = 1, ipl 
    426             DO jj = 2, jpjm1 
     1061         DO jl = 1, jpl 
     1062            DO jj = 1, jpjm1 
    4271063               DO ji = 1, fs_jpim1   ! vector opt. 
    428                   zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    429                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    430                      &                                              -              zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )  
     1064                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     1065                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     1066                     &                                               -              zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )  
    4311067               END DO 
    4321068            END DO 
    … …  
    4351071      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24) 
    4361072         ! 
    437          DO jl = 1, ipl 
    438             DO jj = 2, jpjm1 
     1073         DO jl = 1, jpl 
     1074            DO jj = 1, jpjm1 
    4391075               DO ji = 1, fs_jpim1   ! vector opt. 
    440                   zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     1076                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    4411077                  zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    4421078!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    443                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                     pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        & 
    444                      &                                              -              zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   & 
    445                      &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                        ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        & 
    446                      &                                              - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   ) 
     1079                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        & 
     1080                     &                                               -              zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   & 
     1081                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        & 
     1082                     &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   ) 
    4471083               END DO 
    4481084            END DO 
    … …  
    4511087      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27) 
    4521088         ! 
    453          DO jl = 1, ipl 
    454             DO jj = 2, jpjm1 
     1089         DO jl = 1, jpl 
     1090            DO jj = 1, jpjm1 
    4551091               DO ji = 1, fs_jpim1   ! vector opt. 
    456                   zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     1092                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    4571093                  zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    4581094!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    459                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                     pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        & 
    460                      &                                                    -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   & 
    461                      &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                        ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        & 
    462                      &                                                    - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   ) 
     1095                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        & 
     1096                     &                                               -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   & 
     1097                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        & 
     1098                     &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   ) 
    4631099               END DO 
    4641100            END DO 
    … …  
    4671103      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29) 
    4681104         ! 
    469          DO jl = 1, ipl 
    470             DO jj = 2, jpjm1 
     1105         DO jl = 1, jpl 
     1106            DO jj = 1, jpjm1 
    4711107               DO ji = 1, fs_jpim1   ! vector opt. 
    472                   zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     1108                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    4731109                  zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    474 !!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
     1110                  !!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    4751111                  zdx4 = zdx2 * zdx2 
    476                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (           (                   pt  (ji+1,jj,jl) + pt  (ji,jj,jl)       & 
    477                      &                                                    -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) )   & 
    478                      &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *    (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)       & 
    479                      &                                                    - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) )   & 
    480                      &       + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)       & 
    481                      &                                              - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 
     1112                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (            (                   pt  (ji+1,jj,jl) + pt  (ji,jj,jl)       & 
     1113                     &                                                     -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) )   & 
     1114                     &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)       & 
     1115                     &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) )   & 
     1116                     &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)       & 
     1117                     &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 
    4821118               END DO 
    4831119            END DO 
    … …  
    4851121         ! 
    4861122      END SELECT 
     1123      !                                                     !-- High order flux in i-direction  --! 
     1124      IF( ll_neg ) THEN 
     1125         DO jl = 1, jpl 
     1126            DO jj = 1, jpjm1 
     1127               DO ji = 1, fs_jpim1 
     1128                  IF( pt_u(ji,jj,jl) < 0._wp ) THEN 
     1129                     pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                           pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     1130                        &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
     1131                  ENDIF 
     1132               END DO 
     1133            END DO 
     1134         END DO 
     1135      ENDIF 
     1136 
     1137      DO jl = 1, jpl 
     1138         DO jj = 1, jpjm1 
     1139            DO ji = 1, fs_jpim1   ! vector opt. 
     1140               IF( ll_clem ) THEN 
     1141                  pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 
     1142               ELSE 
     1143                  pfu_ho(ji,jj,jl) = puc(ji,jj,jl) * pt_u(ji,jj,jl) 
     1144               ENDIF 
     1145            END DO 
     1146         END DO 
     1147      END DO 
    4871148      ! 
    4881149   END SUBROUTINE ultimate_x 
    4891150    
    4901151  
    491    SUBROUTINE ultimate_y( k_order, ipl, pdt, pt, pvc, pt_v ) 
     1152   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 
    4921153      !!--------------------------------------------------------------------- 
    4931154      !!                    ***  ROUTINE ultimate_y  *** 
    … …  
    5001161      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    5011162      !!---------------------------------------------------------------------- 
    502       INTEGER                         , INTENT(in   ) ::   k_order   ! ocean time-step index 
    503       INTEGER                         , INTENT(in   ) ::   ipl       ! third dimension of tracer array 
    504       REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step 
    505       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pvc       ! ice j-velocity component 
    506       REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(in   ) ::   pt        ! tracer fields 
    507       REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(  out) ::   pt_v      ! tracer at v-point  
     1163      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
     1164      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step 
     1165      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pv        ! ice j-velocity component 
     1166      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pvc       ! ice j-velocity*A component 
     1167      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt        ! tracer fields 
     1168      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_v      ! tracer at v-point  
     1169      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfv_ho    ! high order flux  
    5081170      ! 
    5091171      INTEGER  ::   ji, jj, jl       ! dummy loop indices 
    5101172      REAL(wp) ::   zcv, zdy2, zdy4    !   -      - 
    511       REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv3 
    512       REAL(wp), DIMENSION(jpi,jpj,ipl) :: ztv2, ztv4 
     1173      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztv1, ztv2, ztv3, ztv4 
    5131174      !!---------------------------------------------------------------------- 
    5141175      ! 
    5151176      !                                                     !--  Laplacian in j-direction  --! 
    516       DO jl = 1, ipl 
     1177      DO jl = 1, jpl 
    5171178         DO jj = 1, jpjm1         ! First derivative (gradient) 
    5181179            DO ji = fs_2, fs_jpim1 
    519                ztv1(ji,jj) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
     1180               ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
    5201181            END DO 
    5211182         END DO 
    5221183         DO jj = 2, jpjm1         ! Second derivative (Laplacian) 
    5231184            DO ji = fs_2, fs_jpim1 
    524                ztv2(ji,jj,jl) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj) 
     1185               ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj) 
    5251186            END DO 
    5261187         END DO 
    … …  
    5291190      ! 
    5301191      !                                                     !--  BiLaplacian in j-direction  --! 
    531       DO jl = 1, ipl 
     1192      DO jl = 1, jpl 
    5321193         DO jj = 1, jpjm1         ! First derivative 
    5331194            DO ji = fs_2, fs_jpim1 
    534             ztv3(ji,jj) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
     1195               ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
    5351196            END DO 
    5361197         END DO 
    5371198         DO jj = 2, jpjm1         ! Second derivative 
    5381199            DO ji = fs_2, fs_jpim1 
    539                ztv4(ji,jj,jl) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj) 
     1200               ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj) 
    5401201            END DO 
    5411202         END DO 
    … …  
    5441205      ! 
    5451206      ! 
    546       SELECT CASE (k_order ) 
    547       ! 
     1207      SELECT CASE (kn_umx ) 
     1208         ! 
    5481209      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21) 
    549          DO jl = 1, ipl 
     1210         DO jl = 1, jpl 
    5501211            DO jj = 1, jpjm1 
    551                DO ji = fs_2, fs_jpim1 
    552                   pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                           ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
    553                      &                                     - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
     1212               DO ji = 1, fs_jpim1 
     1213                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                          ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
     1214                     &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
    5541215               END DO 
    5551216            END DO 
    … …  
    5571218         ! 
    5581219      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23) 
    559          DO jl = 1, ipl 
     1220         DO jl = 1, jpl 
    5601221            DO jj = 1, jpjm1 
    561                DO ji = fs_2, fs_jpim1 
    562                   zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1222               DO ji = 1, fs_jpim1 
     1223                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    5631224                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (     ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
    5641225                     &                                     - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
    … …  
    5691230         ! 
    5701231      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24) 
    571          DO jl = 1, ipl 
     1232         DO jl = 1, jpl 
    5721233            DO jj = 1, jpjm1 
    573                DO ji = fs_2, fs_jpim1 
    574                   zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1234               DO ji = 1, fs_jpim1 
     1235                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    5751236                  zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    5761237!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    … …  
    5841245         ! 
    5851246      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27) 
    586          DO jl = 1, ipl 
     1247         DO jl = 1, jpl 
    5871248            DO jj = 1, jpjm1 
    588                DO ji = fs_2, fs_jpim1 
    589                   zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1249               DO ji = 1, fs_jpim1 
     1250                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    5901251                  zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    5911252!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    … …  
    5991260         ! 
    6001261      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29) 
    601          DO jl = 1, ipl 
     1262         DO jl = 1, jpl 
    6021263            DO jj = 1, jpjm1 
    603                DO ji = fs_2, fs_jpim1 
    604                   zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1264               DO ji = 1, fs_jpim1 
     1265                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    6051266                  zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    6061267!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    … …  
    6171278         ! 
    6181279      END SELECT 
     1280      !                                                     !-- High order flux in j-direction  --! 
     1281      IF( ll_neg ) THEN 
     1282         DO jl = 1, jpl 
     1283            DO jj = 1, jpjm1 
     1284               DO ji = 1, fs_jpim1 
     1285                  IF( pt_v(ji,jj,jl) < 0._wp ) THEN 
     1286                     pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                          ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
     1287                        &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
     1288                  ENDIF 
     1289               END DO 
     1290            END DO 
     1291         END DO 
     1292      ENDIF 
     1293 
     1294      DO jl = 1, jpl 
     1295         DO jj = 1, jpjm1 
     1296            DO ji = 1, fs_jpim1   ! vector opt. 
     1297               IF( ll_clem ) THEN 
     1298                  pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 
     1299               ELSE 
     1300                  pfv_ho(ji,jj,jl) = pvc(ji,jj,jl) * pt_v(ji,jj,jl) 
     1301               ENDIF 
     1302            END DO 
     1303         END DO 
     1304      END DO 
    6191305      ! 
    6201306   END SUBROUTINE ultimate_y 
    621     
    622    
    623    SUBROUTINE nonosc_2d( ipl, pbef, paa, pbb, paft, pdt ) 
     1307      
     1308 
     1309   SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho ) 
    6241310      !!--------------------------------------------------------------------- 
    6251311      !!                    ***  ROUTINE nonosc  *** 
    6261312      !!      
    627       !! **  Purpose :   compute monotonic tracer fluxes from the upstream  
     1313       !! **  Purpose :   compute monotonic tracer fluxes from the upstream  
    6281314      !!       scheme and the before field by a nonoscillatory algorithm  
    6291315      !! 
    6301316      !! **  Method  :   ... ??? 
    631       !!       warning : pbef and paft must be masked, but the boundaries 
     1317      !!       warning : pt and pt_low must be masked, but the boundaries 
    6321318      !!       conditions on the fluxes are not necessary zalezak (1979) 
    6331319      !!       drange (1995) multi-dimensional forward-in-time and upstream- 
    6341320      !!       in-space based differencing for fluid 
    6351321      !!---------------------------------------------------------------------- 
    636       INTEGER                          , INTENT(in   ) ::   ipl          ! third dimension of tracer array 
    637       REAL(wp)                         , INTENT(in   ) ::   pdt          ! tracer time-step 
    638       REAL(wp), DIMENSION (jpi,jpj,ipl), INTENT(in   ) ::   pbef, paft   ! before & after field 
    639       REAL(wp), DIMENSION (jpi,jpj,ipl), INTENT(inout) ::   paa, pbb     ! monotonic fluxes in the 2 directions 
     1322      REAL(wp)                     , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     1323      REAL(wp)                     , INTENT(in   ) ::   pdt              ! tracer time-step 
     1324      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pu               ! ice i-velocity => u*e2 
     1325      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   puc              ! ice i-velocity *A => u*e2*a 
     1326      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pv               ! ice j-velocity => v*e1 
     1327      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pvc              ! ice j-velocity *A => v*e1*a 
     1328      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   ptc, pt, pt_low  ! before field & upstream guess of after field 
     1329      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_low, pfu_low ! upstream flux 
     1330      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux 
    6401331      ! 
    6411332      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    642       INTEGER  ::   ikm1      ! local integer 
    643       REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zsml, z1_dt   ! local scalars 
    644       REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      - 
    645       REAL(wp), DIMENSION(jpi,jpj) :: zbup, zbdo, zmsk 
    646       REAL(wp), DIMENSION(jpi,jpj,ipl) :: zbetup, zbetdo, zdiv 
    647       !!---------------------------------------------------------------------- 
    648       ! 
     1333      REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt, zpos2, zneg2   ! local scalars 
     1334      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign, zcoef        !   -      - 
     1335      REAL(wp), DIMENSION(jpi,jpj    ) :: zbup, zbdo 
     1336      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_low, ztj_low, zzt 
     1337      !!---------------------------------------------------------------------- 
    6491338      zbig = 1.e+40_wp 
    650       zsml = 1.e-15_wp 
    651  
    652       ! test on divergence 
    653       DO jl = 1, ipl 
     1339      zsml = epsi20 
     1340       
     1341      IF( ll_zeroup2 ) THEN 
     1342         DO jl = 1, jpl 
     1343            DO jj = 1, jpjm1 
     1344               DO ji = 1, fs_jpim1   ! vector opt. 
     1345                  IF( amaxu(ji,jj,jl) == 0._wp )   pfu_ho(ji,jj,jl) = 0._wp 
     1346                  IF( amaxv(ji,jj,jl) == 0._wp )   pfv_ho(ji,jj,jl) = 0._wp 
     1347               END DO 
     1348            END DO 
     1349         END DO 
     1350      ENDIF 
     1351 
     1352      IF( ll_zeroup4 ) THEN 
     1353         DO jl = 1, jpl 
     1354            DO jj = 1, jpjm1 
     1355               DO ji = 1, fs_jpim1   ! vector opt. 
     1356                  IF( pfu_low(ji,jj,jl) == 0._wp )   pfu_ho(ji,jj,jl) = 0._wp 
     1357                  IF( pfv_low(ji,jj,jl) == 0._wp )   pfv_ho(ji,jj,jl) = 0._wp 
     1358               END DO 
     1359            END DO 
     1360         END DO 
     1361      ENDIF 
     1362 
     1363 
     1364      IF( ll_zeroup1 ) THEN 
     1365         DO jl = 1, jpl 
     1366            DO jj = 2, jpjm1 
     1367               DO ji = fs_2, fs_jpim1 
     1368                  IF( ll_gurvan ) THEN 
     1369                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1370                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1371                  ELSE 
     1372                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1373                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1374                        &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1375                        &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1376                        &         ) * tmask(ji,jj,1) 
     1377                  ENDIF 
     1378                  IF( zzt(ji,jj,jl) < 0._wp ) THEN 
     1379                     pfu_ho(ji,jj,jl)   = pfu_low(ji,jj,jl) 
     1380                     pfv_ho(ji,jj,jl)   = pfv_low(ji,jj,jl) 
     1381                     WRITE(numout,*) '*** 1 negative high order zzt ***',ji,jj,zzt(ji,jj,jl) 
     1382                  ENDIF 
     1383                  !!               IF( ji==26 .AND. jj==86) THEN 
     1384                  !!                  WRITE(numout,*) 'zzt high order',zzt(ji,jj) 
     1385                  !!                  WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt 
     1386                  !!                  WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt 
     1387                  !!                  WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj,jl)) * r1_e1e2t(ji,jj) * pdt 
     1388                  !!                  WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1,jl)) * r1_e1e2t(ji,jj) * pdt 
     1389                  !!               ENDIF 
     1390                  IF( ll_gurvan ) THEN 
     1391                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1392                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1393                  ELSE 
     1394                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1395                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1396                        &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1397                        &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1398                        &         ) * tmask(ji,jj,1) 
     1399                  ENDIF 
     1400                  IF( zzt(ji,jj,jl) < 0._wp ) THEN 
     1401                     pfu_ho(ji-1,jj,jl) = pfu_low(ji-1,jj,jl) 
     1402                     pfv_ho(ji,jj-1,jl) = pfv_low(ji,jj-1,jl) 
     1403                     WRITE(numout,*) '*** 2 negative high order zzt ***',ji,jj,zzt(ji,jj,jl) 
     1404                  ENDIF 
     1405                  IF( ll_gurvan ) THEN 
     1406                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1407                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1408                  ELSE 
     1409                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1410                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1411                        &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1412                        &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1413                        &         ) * tmask(ji,jj,1) 
     1414                  ENDIF 
     1415                  IF( zzt(ji,jj,jl) < 0._wp ) THEN 
     1416                     WRITE(numout,*) '*** 3 negative high order zzt ***',ji,jj,zzt(ji,jj,jl) 
     1417                  ENDIF 
     1418               END DO 
     1419            END DO 
     1420         END DO 
     1421         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) 
     1422      ENDIF 
     1423 
     1424       
     1425      ! antidiffusive flux : high order minus low order 
     1426      ! -------------------------------------------------- 
     1427      DO jl = 1, jpl 
     1428         DO jj = 1, jpjm1 
     1429            DO ji = 1, fs_jpim1   ! vector opt. 
     1430               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_low(ji,jj,jl) 
     1431               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_low(ji,jj,jl) 
     1432            END DO 
     1433         END DO 
     1434      END DO 
     1435 
     1436      ! extreme case where pfu_ho has to be zero 
     1437      ! ---------------------------------------- 
     1438      !                                    pfu_ho 
     1439      !                           *         ---> 
     1440      !                        |      |  *   |        |  
     1441      !                        |      |      |    *   |     
     1442      !                        |      |      |        |    * 
     1443      !            t_low :       i-1     i       i+1       i+2    
     1444      IF( ll_prelimiter_zalesak ) THEN 
     1445          
     1446         DO jl = 1, jpl 
     1447            DO jj = 2, jpjm1 
     1448               DO ji = fs_2, fs_jpim1  
     1449                  zti_low(ji,jj,jl)= pt_low(ji+1,jj  ,jl) 
     1450                  ztj_low(ji,jj,jl)= pt_low(ji  ,jj+1,jl) 
     1451               END DO 
     1452            END DO 
     1453         END DO 
     1454         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_low, 'T', 1., ztj_low, 'T', 1. ) 
     1455 
     1456         !! this does not work ?? 
     1457         !!            DO jj = 2, jpjm1 
     1458         !!               DO ji = fs_2, fs_jpim1 
     1459         !!                  IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji+1,jj  ) - pt_low (ji  ,jj) ) .AND.     & 
     1460         !!               &      SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj+1) - pt_low (ji  ,jj) )           & 
     1461         !!               &    ) THEN 
     1462         !!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj ) - zti_low(ji  ,jj) ) .AND.   & 
     1463         !!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1 ) - ztj_low(ji  ,jj) )         & 
     1464         !!               &       ) THEN 
     1465         !!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0. 
     1466         !!                     ENDIF 
     1467         !!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji-1,jj  ) ) .AND.  & 
     1468         !!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji  ,jj-1) )        & 
     1469         !!               &       )  THEN 
     1470         !!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0. 
     1471         !!                     ENDIF 
     1472         !!                  ENDIF 
     1473         !!                END DO 
     1474         !!            END DO             
     1475 
     1476         DO jl = 1, jpl 
     1477            DO jj = 2, jpjm1 
     1478               DO ji = fs_2, fs_jpim1 
     1479                  IF ( pfu_ho(ji,jj,jl) * ( pt_low(ji+1,jj,jl) - pt_low(ji,jj,jl) ) <= 0. .AND.  & 
     1480                     & pfv_ho(ji,jj,jl) * ( pt_low(ji,jj+1,jl) - pt_low(ji,jj,jl) ) <= 0. ) THEN 
     1481                     ! 
     1482                     IF(  pfu_ho(ji,jj,jl) * ( zti_low(ji+1,jj,jl) - zti_low(ji,jj,jl) ) <= 0 .AND.  & 
     1483                        & pfv_ho(ji,jj,jl) * ( ztj_low(ji,jj+1,jl) - ztj_low(ji,jj,jl) ) <= 0)  pfu_ho(ji,jj,jl)=0. ; pfv_ho(ji,jj,jl)=0. 
     1484                     ! 
     1485                     IF(  pfu_ho(ji,jj,jl) * ( pt_low(ji  ,jj,jl) - pt_low(ji-1,jj,jl) ) <= 0 .AND.  & 
     1486                        & pfv_ho(ji,jj,jl) * ( pt_low(ji  ,jj,jl) - pt_low(ji,jj-1,jl) ) <= 0)  pfu_ho(ji,jj,jl)=0. ; pfv_ho(ji,jj,jl)=0. 
     1487                     ! 
     1488                  ENDIF 
     1489               END DO 
     1490            END DO 
     1491         END DO 
     1492         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond. 
     1493 
     1494      ELSEIF( ll_prelimiter_devore ) THEN 
     1495         DO jl = 1, jpl 
     1496            DO jj = 2, jpjm1 
     1497               DO ji = fs_2, fs_jpim1  
     1498                  zti_low(ji,jj,jl)= pt_low(ji+1,jj  ,jl) 
     1499                  ztj_low(ji,jj,jl)= pt_low(ji  ,jj+1,jl) 
     1500               END DO 
     1501            END DO 
     1502         END DO 
     1503         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_low, 'T', 1., ztj_low, 'T', 1. ) 
     1504 
     1505         z1_dt = 1._wp / pdt 
     1506         DO jl = 1, jpl 
     1507            DO jj = 2, jpjm1 
     1508               DO ji = fs_2, fs_jpim1 
     1509                  zsign = SIGN( 1., pt_low(ji+1,jj,jl) - pt_low(ji,jj,jl) ) 
     1510                  pfu_ho(ji,jj,jl) =  zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj,jl)) , & 
     1511                     &                          zsign * ( pt_low (ji  ,jj,jl) - pt_low (ji-1,jj,jl) ) * e1e2t(ji  ,jj) * z1_dt , & 
     1512                     &                          zsign * ( zti_low(ji+1,jj,jl) - zti_low(ji  ,jj,jl) ) * e1e2t(ji+1,jj) * z1_dt ) ) 
     1513 
     1514                  zsign = SIGN( 1., pt_low(ji,jj+1,jl) - pt_low(ji,jj,jl) ) 
     1515                  pfv_ho(ji,jj,jl) =  zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj,jl)) , & 
     1516                     &                          zsign * ( pt_low (ji,jj  ,jl) - pt_low (ji,jj-1,jl) ) * e1e2t(ji,jj  ) * z1_dt , & 
     1517                     &                          zsign * ( ztj_low(ji,jj+1,jl) - ztj_low(ji,jj  ,jl) ) * e1e2t(ji,jj+1) * z1_dt ) ) 
     1518               END DO 
     1519            END DO 
     1520         END DO 
     1521         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond. 
     1522 
     1523      ENDIF 
     1524 
     1525 
     1526      ! Search local extrema 
     1527      ! -------------------- 
     1528      ! max/min of pt & pt_low with large negative/positive value (-/+zbig) outside ice cover 
     1529      z1_dt = 1._wp / pdt 
     1530      DO jl = 1, jpl 
     1531          
     1532         DO jj = 1, jpj 
     1533            DO ji = 1, jpi 
     1534               IF    ( pt(ji,jj,jl) <= 0._wp .AND. pt_low(ji,jj,jl) <= 0._wp ) THEN 
     1535                  zbup(ji,jj) = -zbig 
     1536                  zbdo(ji,jj) =  zbig 
     1537               ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_low(ji,jj,jl) > 0._wp ) THEN 
     1538                  zbup(ji,jj) = pt_low(ji,jj,jl) 
     1539                  zbdo(ji,jj) = pt_low(ji,jj,jl) 
     1540               ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_low(ji,jj,jl) <= 0._wp ) THEN 
     1541                  zbup(ji,jj) = pt(ji,jj,jl) 
     1542                  zbdo(ji,jj) = pt(ji,jj,jl) 
     1543               ELSE 
     1544                  zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_low(ji,jj,jl) ) 
     1545                  zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_low(ji,jj,jl) ) 
     1546               ENDIF 
     1547            END DO 
     1548         END DO 
     1549 
    6541550         DO jj = 2, jpjm1 
    655             DO ji = fs_2, fs_jpim1   ! vector opt.   
    656                zdiv(ji,jj,jl) =  - (  paa(ji,jj,jl) - paa(ji-1,jj  ,jl)   & 
    657                   &              +    pbb(ji,jj,jl) - pbb(ji  ,jj-1,jl) )   
    658             END DO 
    659          END DO 
    660       END DO 
    661       CALL lbc_lnk( 'icedyn_adv_umx', zdiv, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign) 
    662  
    663       DO jl = 1, ipl 
    664          ! Determine ice masks for before and after tracers  
    665          WHERE( pbef(:,:,jl) == 0._wp .AND. paft(:,:,jl) == 0._wp .AND. zdiv(:,:,jl) == 0._wp )  
    666             zmsk(:,:) = 0._wp 
    667          ELSEWHERE                                                                                    
    668             zmsk(:,:) = 1._wp * tmask(:,:,1) 
    669          END WHERE 
    670  
    671          ! Search local extrema 
    672          ! -------------------- 
    673          ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 
    674 !         zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ),   & 
    675 !            &             paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) )  ) 
    676 !         zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ),   & 
    677 !            &             paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) )  ) 
    678          zbup(:,:) = MAX( pbef(:,:,jl) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ),   & 
    679             &             paft(:,:,jl) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) )  ) 
    680          zbdo(:,:) = MIN( pbef(:,:,jl) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ),   & 
    681             &             paft(:,:,jl) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) )  ) 
    682  
    683       z1_dt = 1._wp / pdt 
    684       DO jj = 2, jpjm1 
    685          DO ji = fs_2, fs_jpim1   ! vector opt. 
     1551            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1552               ! 
     1553               IF( .NOT. ll_9points ) THEN 
     1554                  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 
     1555                  zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1) ) 
     1556                  ! 
     1557               ELSE 
     1558                  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 
     1559                     &                     zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1)  ) 
     1560                  zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1), & 
     1561                     &                     zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1)  ) 
     1562               ENDIF 
     1563               ! 
     1564               zpos = MAX( 0., pfu_ho(ji-1,jj,jl) ) - MIN( 0., pfu_ho(ji  ,jj,jl) ) &  ! positive/negative part of the flux 
     1565                  & + MAX( 0., pfv_ho(ji,jj-1,jl) ) - MIN( 0., pfv_ho(ji,jj  ,jl) ) 
     1566               zneg = MAX( 0., pfu_ho(ji  ,jj,jl) ) - MIN( 0., pfu_ho(ji-1,jj,jl) ) & 
     1567                  & + MAX( 0., pfv_ho(ji,jj  ,jl) ) - MIN( 0., pfv_ho(ji,jj-1,jl) ) 
     1568               ! 
     1569               IF( ll_HgradU .AND. .NOT.ll_gurvan ) THEN 
     1570                  zneg2 = (   pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) & 
     1571                     &    ) * ( 1. - pamsk ) 
     1572                  zpos2 = ( - pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) - pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) & 
     1573                     &    ) * ( 1. - pamsk ) 
     1574               ELSE 
     1575                  zneg2 = 0. ; zpos2 = 0. 
     1576               ENDIF 
     1577               ! 
     1578               !                                  ! up & down beta terms 
     1579               IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_low(ji,jj,jl) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt 
     1580               ELSE                         ; zbetup(ji,jj,jl) = 0. ! zbig 
     1581               ENDIF 
     1582               ! 
     1583               IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_low(ji,jj,jl) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt 
     1584               ELSE                         ; zbetdo(ji,jj,jl) = 0. ! zbig 
     1585               ENDIF 
     1586               ! 
     1587               ! if all the points are outside ice cover 
     1588               IF( zup == -zbig )   zbetup(ji,jj,jl) = 0. ! zbig 
     1589               IF( zdo ==  zbig )   zbetdo(ji,jj,jl) = 0. ! zbig             
     1590               ! 
     1591 
     1592!!            IF( ji==26 .AND. jj==86) THEN 
     1593!               WRITE(numout,*) '-----------------' 
     1594!               WRITE(numout,*) 'zpos',zpos,zpos2 
     1595!               WRITE(numout,*) 'zneg',zneg,zneg2 
     1596!               WRITE(numout,*) 'puc/pu',ABS(puc(ji,jj))/MAX(epsi20, ABS(pu(ji,jj))) 
     1597!               WRITE(numout,*) 'pvc/pv',ABS(pvc(ji,jj))/MAX(epsi20, ABS(pv(ji,jj))) 
     1598!               WRITE(numout,*) 'pucm1/pu',ABS(puc(ji-1,jj))/MAX(epsi20, ABS(pu(ji-1,jj))) 
     1599!               WRITE(numout,*) 'pvcm1/pv',ABS(pvc(ji,jj-1))/MAX(epsi20, ABS(pv(ji,jj-1))) 
     1600!               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)+pfu_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1601!               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)+pfv_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1602!               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)+pfu_low(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt 
     1603!               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)+pfv_low(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt 
     1604!               WRITE(numout,*) 'pfu_low',pfu_low(ji,jj) * r1_e1e2t(ji,jj) * pdt 
     1605!               WRITE(numout,*) 'pfv_low',pfv_low(ji,jj) * r1_e1e2t(ji,jj) * pdt 
     1606!               WRITE(numout,*) 'pfu_lowm1',pfu_low(ji-1,jj) * r1_e1e2t(ji,jj) * pdt 
     1607!               WRITE(numout,*) 'pfv_lowm1',pfv_low(ji,jj-1) * r1_e1e2t(ji,jj) * pdt 
     1608!                
     1609!               WRITE(numout,*) 'pt',pt(ji,jj) 
     1610!               WRITE(numout,*) 'ptim1',pt(ji-1,jj) 
     1611!               WRITE(numout,*) 'ptjm1',pt(ji,jj-1) 
     1612!               WRITE(numout,*) 'pt_low',pt_low(ji,jj) 
     1613!               WRITE(numout,*) 'zbetup',zbetup(ji,jj) 
     1614!               WRITE(numout,*) 'zbetdo',zbetdo(ji,jj) 
     1615!               WRITE(numout,*) 'zup',zup 
     1616!               WRITE(numout,*) 'zdo',zdo 
     1617!            ENDIF 
    6861618            ! 
    687             zup  = MAX(   zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ),   &        ! search max/min in neighbourhood 
    688                &                       zbup(ji  ,jj-1), zbup(ji  ,jj+1)    ) 
    689             zdo  = MIN(   zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ),   & 
    690                &                       zbdo(ji  ,jj-1), zbdo(ji  ,jj+1)    ) 
    691                ! 
    692                zpos = MAX( 0., paa(ji-1,jj  ,jl) ) - MIN( 0., paa(ji  ,jj  ,jl) )   &        ! positive/negative  part of the flux 
    693                   & + MAX( 0., pbb(ji  ,jj-1,jl) ) - MIN( 0., pbb(ji  ,jj  ,jl) ) 
    694                zneg = MAX( 0., paa(ji  ,jj  ,jl) ) - MIN( 0., paa(ji-1,jj  ,jl) )   & 
    695                   & + MAX( 0., pbb(ji  ,jj  ,jl) ) - MIN( 0., pbb(ji  ,jj-1,jl) ) 
    696                ! 
    697                zbt = e1e2t(ji,jj) * z1_dt                                   ! up & down beta terms 
    698                zbetup(ji,jj,jl) = ( zup            - paft(ji,jj,jl) ) / ( zpos + zsml ) * zbt 
    699                zbetdo(ji,jj,jl) = ( paft(ji,jj,jl) - zdo         )    / ( zneg + zsml ) * zbt 
    7001619            END DO 
    7011620         END DO 
    … …  
    7031622      CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign) 
    7041623 
    705       DO jl = 1, ipl 
    706          ! monotonic flux in the i & j direction (paa & pbb) 
    707          ! ------------------------------------- 
    708          DO jj = 2, jpjm1 
     1624       
     1625      ! monotonic flux in the y direction 
     1626      ! --------------------------------- 
     1627      DO jl = 1, jpl 
     1628         DO jj = 1, jpjm1 
    7091629            DO ji = 1, fs_jpim1   ! vector opt. 
    7101630               zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 
    7111631               zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 
    712                zcu = 0.5  + SIGN( 0.5 , paa(ji,jj,jl) ) 
     1632               zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj,jl) ) 
    7131633               ! 
    714                paa(ji,jj,jl) = paa(ji,jj,jl) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 
    715             END DO 
    716          END DO 
    717          ! 
     1634               zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu ) 
     1635 
     1636               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_low(ji,jj,jl) 
     1637 
     1638!!            IF( ji==26 .AND. jj==86) THEN 
     1639!!               WRITE(numout,*) 'coefU',zcoef 
     1640!!               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt 
     1641!!               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj,jl)) * r1_e1e2t(ji,jj) * pdt 
     1642!!            ENDIF 
     1643 
     1644            END DO 
     1645         END DO 
     1646 
    7181647         DO jj = 1, jpjm1 
    719             DO ji = fs_2, fs_jpim1   ! vector opt. 
     1648            DO ji = 1, fs_jpim1   ! vector opt. 
    7201649               zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 
    7211650               zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 
    722                zcv = 0.5  + SIGN( 0.5 , pbb(ji,jj,jl) ) 
     1651               zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj,jl) ) 
    7231652               ! 
    724                pbb(ji,jj,jl) = pbb(ji,jj,jl) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 
    725             END DO 
    726          END DO 
     1653               zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv ) 
     1654 
     1655               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_low(ji,jj,jl) 
     1656 
     1657!!            IF( ji==26 .AND. jj==86) THEN 
     1658!!               WRITE(numout,*) 'coefV',zcoef 
     1659!!               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt 
     1660!!               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1,jl)) * r1_e1e2t(ji,jj) * pdt 
     1661!!            ENDIF 
     1662            END DO 
     1663         END DO 
     1664 
     1665         ! clem test 
     1666         DO jj = 2, jpjm1 
     1667            DO ji = 2, fs_jpim1   ! vector opt. 
     1668               IF( ll_gurvan ) THEN 
     1669                  zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1670                     &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1671               ELSE 
     1672                  zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1673                     &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1674                     &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1675                     &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1676                     &         ) * tmask(ji,jj,1) 
     1677               ENDIF 
     1678               IF( zzt(ji,jj,jl) < -epsi20 ) THEN 
     1679                  WRITE(numout,*) 'T<0 nonosc',zzt(ji,jj,jl) 
     1680               ENDIF 
     1681            END DO 
     1682         END DO 
     1683 
    7271684      END DO 
     1685 
     1686      ! 
    7281687      ! 
    7291688   END SUBROUTINE nonosc_2d 
     1689 
     1690   SUBROUTINE limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     1691      !!--------------------------------------------------------------------- 
     1692      !!                    ***  ROUTINE limiter_x  *** 
     1693      !!      
     1694      !! **  Purpose :   compute flux limiter  
     1695      !!---------------------------------------------------------------------- 
     1696      REAL(wp)                   , INTENT(in   )           ::   pdt          ! tracer time-step 
     1697      REAL(wp), DIMENSION (:,:  ), INTENT(in   )           ::   pu           ! ice i-velocity => u*e2 
     1698      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   puc          ! ice i-velocity *A => u*e2*a 
     1699      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pt           ! ice tracer 
     1700      REAL(wp), DIMENSION (:,:,:), INTENT(inout)           ::   pfu_ho       ! high order flux 
     1701      REAL(wp), DIMENSION (:,:,:), INTENT(in   ), OPTIONAL ::   pfu_ups      ! upstream flux 
     1702      ! 
     1703      REAL(wp) ::   Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr 
     1704      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     1705      REAL(wp), DIMENSION (jpi,jpj,jpl) ::   zslpx       ! tracer slopes  
     1706      !!---------------------------------------------------------------------- 
     1707      ! 
     1708      DO jl = 1, jpl 
     1709         DO jj = 2, jpjm1 
     1710            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1711               zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1) 
     1712            END DO 
     1713         END DO 
     1714      END DO 
     1715      CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.)   ! lateral boundary cond. 
     1716       
     1717      DO jl = 1, jpl 
     1718         DO jj = 2, jpjm1 
     1719            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1720               uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) 
     1721                
     1722               Rjm = zslpx(ji-1,jj,jl) 
     1723               Rj  = zslpx(ji  ,jj,jl) 
     1724               Rjp = zslpx(ji+1,jj,jl) 
     1725 
     1726               IF( PRESENT(pfu_ups) ) THEN 
     1727 
     1728                  IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
     1729                  ELSE                        ;   Rr = Rjp 
     1730                  ENDIF 
     1731 
     1732                  zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)      
     1733                  IF( Rj > 0. ) THEN 
     1734                     zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)),  & 
     1735                        &        MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 
     1736                  ELSE 
     1737                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)),  & 
     1738                        &        MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 
     1739                  ENDIF 
     1740                  pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 
     1741 
     1742               ELSE 
     1743                  IF( Rj /= 0. ) THEN 
     1744                     IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj 
     1745                     ELSE                        ;   Cr = Rjp / Rj 
     1746                     ENDIF 
     1747                  ELSE 
     1748                     Cr = 0. 
     1749                     !IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20 
     1750                     !ELSE                        ;   Cr = Rjp * 1.e20 
     1751                     !ENDIF 
     1752                  ENDIF 
     1753 
     1754                  ! -- superbee -- 
     1755                  zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 
     1756                  ! -- van albada 2 -- 
     1757                  !!zpsi = 2.*Cr / (Cr*Cr+1.) 
     1758 
     1759                  ! -- sweby (with beta=1) -- 
     1760                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 
     1761                  ! -- van Leer -- 
     1762                  !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 
     1763                  ! -- ospre -- 
     1764                  !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 
     1765                  ! -- koren -- 
     1766                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 
     1767                  ! -- charm -- 
     1768                  !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 
     1769                  !ELSE                 ;   zpsi = 0. 
     1770                  !ENDIF 
     1771                  ! -- van albada 1 -- 
     1772                  !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 
     1773                  ! -- smart -- 
     1774                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 
     1775                  ! -- umist -- 
     1776                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 
     1777 
     1778                  ! high order flux corrected by the limiter 
     1779                  pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5 
     1780 
     1781               ENDIF 
     1782            END DO 
     1783         END DO 
     1784      END DO 
     1785      CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.)   ! lateral boundary cond. 
     1786      ! 
     1787   END SUBROUTINE limiter_x 
     1788 
     1789   SUBROUTINE limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     1790      !!--------------------------------------------------------------------- 
     1791      !!                    ***  ROUTINE limiter_y  *** 
     1792      !!      
     1793      !! **  Purpose :   compute flux limiter  
     1794      !!---------------------------------------------------------------------- 
     1795      REAL(wp)                   , INTENT(in   )           ::   pdt          ! tracer time-step 
     1796      REAL(wp), DIMENSION (:,:  ), INTENT(in   )           ::   pv           ! ice i-velocity => u*e2 
     1797      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pvc          ! ice i-velocity *A => u*e2*a 
     1798      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pt           ! ice tracer 
     1799      REAL(wp), DIMENSION (:,:,:), INTENT(inout)           ::   pfv_ho       ! high order flux 
     1800      REAL(wp), DIMENSION (:,:,:), INTENT(in   ), OPTIONAL ::   pfv_ups      ! upstream flux 
     1801      ! 
     1802      REAL(wp) ::   Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr 
     1803      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     1804      REAL(wp), DIMENSION (jpi,jpj,jpl) ::   zslpy       ! tracer slopes  
     1805      !!---------------------------------------------------------------------- 
     1806      ! 
     1807      DO jl = 1, jpl 
     1808         DO jj = 2, jpjm1 
     1809            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1810               zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1) 
     1811            END DO 
     1812         END DO 
     1813      END DO 
     1814      CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.)   ! lateral boundary cond. 
     1815 
     1816      DO jl = 1, jpl 
     1817         DO jj = 2, jpjm1 
     1818            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1819               vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) 
     1820 
     1821               Rjm = zslpy(ji,jj-1,jl) 
     1822               Rj  = zslpy(ji,jj  ,jl) 
     1823               Rjp = zslpy(ji,jj+1,jl) 
     1824 
     1825               IF( PRESENT(pfv_ups) ) THEN 
     1826 
     1827                  IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
     1828                  ELSE                        ;   Rr = Rjp 
     1829                  ENDIF 
     1830 
     1831                  zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)      
     1832                  IF( Rj > 0. ) THEN 
     1833                     zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)),  & 
     1834                        &        MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 
     1835                  ELSE 
     1836                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)),  & 
     1837                        &        MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 
     1838                  ENDIF 
     1839                  pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 
     1840 
     1841               ELSE 
     1842 
     1843                  IF( Rj /= 0. ) THEN 
     1844                     IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj 
     1845                     ELSE                        ;   Cr = Rjp / Rj 
     1846                     ENDIF 
     1847                  ELSE 
     1848                     Cr = 0. 
     1849                     !IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20 
     1850                     !ELSE                        ;   Cr = Rjp * 1.e20 
     1851                     !ENDIF 
     1852                  ENDIF 
     1853 
     1854                  ! -- superbee -- 
     1855                  zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 
     1856                  ! -- van albada 2 -- 
     1857                  !!zpsi = 2.*Cr / (Cr*Cr+1.) 
     1858 
     1859                  ! -- sweby (with beta=1) -- 
     1860                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 
     1861                  ! -- van Leer -- 
     1862                  !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 
     1863                  ! -- ospre -- 
     1864                  !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 
     1865                  ! -- koren -- 
     1866                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 
     1867                  ! -- charm -- 
     1868                  !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 
     1869                  !ELSE                 ;   zpsi = 0. 
     1870                  !ENDIF 
     1871                  ! -- van albada 1 -- 
     1872                  !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 
     1873                  ! -- smart -- 
     1874                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 
     1875                  ! -- umist -- 
     1876                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 
     1877 
     1878                  ! high order flux corrected by the limiter 
     1879                  pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5 
     1880 
     1881               ENDIF 
     1882            END DO 
     1883         END DO 
     1884      END DO 
     1885      CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.)   ! lateral boundary cond. 
     1886      ! 
     1887   END SUBROUTINE limiter_y 
    7301888 
    7311889#else 

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

    r10402 r10419  
    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/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rhg.F90

    r10069 r10419  
    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/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rhg_evp.F90

    r10402 r10419  
    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                            
    … …  
    144147      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components 
    145148!!$      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence 
    146       REAL(wp), DIMENSION(jpi,jpj) ::   zssh_lead_m                     ! array used for the calculation of ice surface slope: 
     149      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope: 
    147150      !                                                                 !    ocean surface (ssh_m) if ice is not embedded 
    148       !                                                                 !    ice top surface if ice is embedded    
     151      !                                                                 !    ice bottom surface if ice is embedded    
    149152      REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array 
    150153      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array 
    … …  
    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      !------------------------------------------------------------------------------! 
    261268      ! 2) Wind / ocean stress, mass terms, coriolis terms 
    262269      !------------------------------------------------------------------------------! 
    263  
    264       !== embedded sea ice: compute representative ice top surface      ==! 
    265       !== non-embedded sea ice: use ocean surface for slope calculation ==! 
    266       zssh_lead_m(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 
     270      ! sea surface height 
     271      !    embedded sea ice: compute representative ice top surface 
     272      !    non-embedded sea ice: use ocean surface for slope calculation 
     273      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 
    267274 
    268275      DO jj = 2, jpjm1 
    … …  
    302309 
    303310            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 
    304             zspgU(ji,jj)    = - zmassU * grav * ( zssh_lead_m(ji+1,jj) - zssh_lead_m(ji,jj) ) * r1_e1u(ji,jj) 
    305             zspgV(ji,jj)    = - zmassV * grav * ( zssh_lead_m(ji,jj+1) - zssh_lead_m(ji,jj) ) * r1_e2v(ji,jj) 
     311            zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj) 
     312            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj) 
    306313 
    307314            ! masks 
    … …  
    317324      CALL lbc_lnk_multi( 'icedyn_rhg_evp', 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( 'icedyn_rhg_evp', 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 
    … …  
    382428               ENDIF 
    383429                
    384                ! stress at T points 
    385                zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1 
    386                zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2 
     430               ! stress at T points (zkt/=0 if landfast) 
     431               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 
     432               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 
    387433              
    388434            END DO 
    … …  
    403449               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) ) 
    404450                
    405                ! stress at F points 
    406                zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2 
     451               ! stress at F points (zkt/=0 if landfast) 
     452               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 
    407453 
    408454            END DO 
    … …  
    452498                  ! 
    453499                  !                 !--- tau_bottom/v_ice 
    454                   zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
    455                   zTauB = - tau_icebfr(ji,jj) / zvel 
     500                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
     501                  zTauB = - zTauV_ib(ji,jj) / zvel 
    456502                  ! 
    457503                  !                 !--- Coriolis at V-points (energy conserving formulation) 
    … …  
    464510                  ! 
    465511                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    466                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     512                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    467513                  ! 
    468514                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    … …  
    500546                  ! 
    501547                  !                 !--- tau_bottom/u_ice 
    502                   zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
    503                   zTauB = - tau_icebfr(ji,jj) / zvel 
     548                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
     549                  zTauB = - zTauU_ib(ji,jj) / zvel 
    504550                  ! 
    505551                  !                 !--- Coriolis at U-points (energy conserving formulation) 
    … …  
    512558                  ! 
    513559                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    514                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     560                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    515561                  ! 
    516562                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    … …  
    550596                  ! 
    551597                  !                 !--- tau_bottom/u_ice 
    552                   zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
    553                   zTauB = - tau_icebfr(ji,jj) / zvel 
     598                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
     599                  zTauB = - zTauU_ib(ji,jj) / zvel 
    554600                  ! 
    555601                  !                 !--- Coriolis at U-points (energy conserving formulation) 
    … …  
    562608                  ! 
    563609                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    564                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     610                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    565611                  ! 
    566612                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    … …  
    598644                  ! 
    599645                  !                 !--- tau_bottom/v_ice 
    600                   zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
    601                   ztauB = - tau_icebfr(ji,jj) / zvel 
     646                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
     647                  zTauB = - zTauV_ib(ji,jj) / zvel 
    602648                  ! 
    603649                  !                 !--- Coriolis at v-points (energy conserving formulation) 
    … …  
    610656                  ! 
    611657                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    612                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     658                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    613659                  ! 
    614660                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 

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

    r10358 r10419  
    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/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icestp.F90

    r10402 r10419  
    194194                                        CALL ice_var_agg( 2 )         ! necessary calls (at least for coupling) 
    195195         ! 
    196 !! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work 
    197 !!       moreover it should only be called at the update frequency (as in agrif_ice_update.F90) 
    198 !# if defined key_agrif 
    199 !         IF( .NOT. Agrif_Root() )     CALL Agrif_ChildGrid_To_ParentGrid() 
    200 !# endif 
    201196                                        CALL ice_update_flx( kt )     ! -- Update ocean surface mass, heat and salt fluxes 
    202 !# if defined key_agrif 
    203 !         IF( .NOT. Agrif_Root() )     CALL Agrif_ParentGrid_To_ChildGrid() 
    204 !# endif 
     197         ! 
    205198         IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs  
    206199         ! 
    … …  
    368361      e_i_b (:,:,:,:) = e_i (:,:,:,:)   ! ice thermal energy 
    369362      WHERE( a_i_b(:,:,:) >= epsi20 ) 
    370          h_i_b(:,:,:) = v_i_b (:,:,:) / a_i_b(:,:,:)   ! ice thickness 
    371          h_s_b(:,:,:) = v_s_b (:,:,:) / a_i_b(:,:,:)   ! snw thickness 
     363         h_i_b(:,:,:) = v_i_b(:,:,:) / a_i_b(:,:,:)   ! ice thickness 
     364         h_s_b(:,:,:) = v_s_b(:,:,:) / a_i_b(:,:,:)   ! snw thickness 
    372365      ELSEWHERE 
    373366         h_i_b(:,:,:) = 0._wp 
    374367         h_s_b(:,:,:) = 0._wp 
     368      END WHERE 
     369       
     370      WHERE( a_ip(:,:,:) >= epsi20 ) 
     371         h_ip_b(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)   ! ice pond thickness 
     372      ELSEWHERE 
     373         h_ip_b(:,:,:) = 0._wp 
    375374      END WHERE 
    376375      ! 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icevar.F90

    r10345 r10419  
    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 
    … …  
    953954   FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b) 
    954955      !!--------------------------------------------------------------------- 
    955       !!                   ***  ROUTINE rhg_evp_rst  *** 
     956      !!                   ***  ROUTINE ice_var_sshdyn  *** 
    956957      !!                      
    957958      !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded 
    … …  
    962963      !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*,  
    963964      !!                Ocean Modelling, Volume 24, Issues 1-2, 2008 
    964       !! 
    965965      !!---------------------------------------------------------------------- 
    966966      ! 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icewri.F90

    r10358 r10419  
    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( 'icewri', 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 
Note: See TracChangeset for help on using the changeset viewer.