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 10910 – NEMO

Changeset 10910


Ignore:
Timestamp:
2019-04-29T18:10:38+02:00 (5 years ago)
Author:
clem
Message:

Major change: the advection scheme UMx has been revisited to clean all the unphysical values which occured. Minor change: the ref config SPITZ12 has been slightly modified to test more parameterizations that are available in the code (melt ponds and landfast)

Location:
NEMO/releases/release-4.0
Files:
11 edited

Legend:

Unmodified
Added
Removed
  • NEMO/releases/release-4.0/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml

    r9943 r10910  
    7979       <field field_ref="vfxsnw"           name="vfxsnw" /> 
    8080        
    81        <!-- diag error for negative ice volume after advection --> 
    82        <field field_ref="iceneg_pres"      name="sineg_pres" /> 
    83        <field field_ref="iceneg_volu"      name="sineg_volu" /> 
    84        <field field_ref="iceneg_hfx"       name="sineg_hfx"  /> 
    85  
    8681       <!-- categories --> 
    8782       <field field_ref="icemask_cat"      name="simskcat"/> 
  • NEMO/releases/release-4.0/cfgs/SHARED/field_def_nemo-ice.xml

    r10578 r10910  
    160160          <field id="hfxdhc"       long_name="Heat content variation in snow and ice (neg = ice cooling)"   unit="W/m2" /> 
    161161 
    162      <!-- diagnostics of the negative values resulting from the advection scheme --> 
    163      <field id="iceneg_pres"  long_name="Fraction of time steps with negative sea ice volume"                    unit=""     /> 
    164      <field id="iceneg_volu"  long_name="Negative sea ice volume per area arising from advection"                unit="m"    /> 
    165           <field id="iceneg_hfx"   long_name="Negative sea ice heat content (eq. heat flux) arising from advection"   unit="W/m2" /> 
    166  
    167162     <!-- sbcssm variables --> 
    168163          <field id="sst_m"    unit="degC" /> 
  • NEMO/releases/release-4.0/cfgs/SHARED/namelist_ice_ref

    r10612 r10910  
    154154&namthd_do      !   Ice growth in open water 
    155155!------------------------------------------------------------------------------ 
    156    rn_hinew         =   0.1           !  thickness for new ice formation in open water (m), must be larger than rn_hnewice 
     156   rn_hinew         =   0.1           !  thickness for new ice formation in open water (m), must be larger than rn_himin 
    157157   ln_frazil        = .false.         !  Frazil ice parameterization (ice collection as a function of wind) 
    158158      rn_maxfraz    =   1.0           !     maximum fraction of frazil ice collecting at the ice base 
  • NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml

    r9943 r10910  
    7979       <field field_ref="vfxsnw"           name="vfxsnw" /> 
    8080 
    81        <!-- diag error for negative ice volume after advection --> 
    82        <field field_ref="iceneg_pres"      name="sineg_pres" /> 
    83        <field field_ref="iceneg_volu"      name="sineg_volu" /> 
    84        <field field_ref="iceneg_hfx"       name="sineg_hfx"  /> 
    85         
    8681       <!-- categories --> 
    8782       <field field_ref="icemask_cat"      name="simskcat"/> 
  • NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/namelist_cfg

    r10075 r10910  
    8080                     ! Sea-ice : 
    8181   nn_ice      = 2         !  SI3 
     82   ln_ice_embd = .false.   !  =T embedded sea-ice (pressure + mass and salt exchanges) 
     83      !                    !  =F levitating ice (no pressure, mass and salt exchanges) 
    8284                     ! Misc. options of sbc : 
    8385   ln_traqsr   = .true.    !  Light penetration in the ocean            (T => fill namtra_qsr ) 
     
    9092      ! 
    9193      ln_Cd_L12   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2012) 
    92       ln_Cd_L15   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2015) 
     94      ln_Cd_L15   = .true.    !  air-ice drags = F(ice concentration) (Lupkes et al. 2015) 
    9395      ! 
    9496   cn_dir = './'  !  root directory for the bulk data location 
     
    9698   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ !          weights filename            ! rotation ! land/sea mask ! 
    9799   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                                      ! pairing  !    filename   ! 
    98 !! ERAI 
    99 !   sn_wndi     = 'u10_era_spitz'               ,         6         , 'u10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Uwnd' , '' 
    100 !   sn_wndj     = 'v10_era_spitz'               ,         6         , 'v10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Vwnd' , '' 
    101 !   sn_qsr      = 'ssrd_era_spitz'              ,         6         , 'ssrd'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    102 !   sn_qlw      = 'strd_era_spitz'              ,         6         , 'strd'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    103 !   sn_tair     = 't2_era_spitz'                ,         6         , 't2'      ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    104 !   sn_humi     = 'humi_era_spitz'              ,         6         , 'humi'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    105 !   sn_prec     = 'precip_era_spitz'            ,         6         , 'precip'  ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    106 !   sn_snow     = 'snow_era_spitz'              ,         6         , 'snow'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    107 !   sn_tdif     = 'taudif'                      ,         6         , 'taudif'  ,   .true.    , .false. , 'yearly'  , ''             , '' , '' 
    108 !      rn_zqt      = 10.       !  Air temperature & humidity reference height (m) 
    109 !! MAR 
    110100   sn_wndi     = 'MARv3.6-9km-Svalbard-2hourly_spitz' ,   2 ,  'u10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Uwnd' , '' 
    111101   sn_wndj     = 'MARv3.6-9km-Svalbard-2hourly_spitz' ,   2 ,  'v10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Vwnd' , '' 
     
    129119   !                       ! type of penetration                        (default: NO selection) 
    130120   ln_qsr_rgb  = .true.       !  RGB light penetration (Red-Green-Blue) 
    131 / 
    132 !----------------------------------------------------------------------- 
    133 &namsbc_ssr    !   surface boundary condition : sea surface restoring   (ln_ssr =T) 
    134 !----------------------------------------------------------------------- 
    135 / 
    136 !----------------------------------------------------------------------- 
    137 &namsbc_rnf    !   runoffs                                              (ln_rnf =T) 
    138 !----------------------------------------------------------------------- 
    139121/ 
    140122!!====================================================================== 
     
    236218!----------------------------------------------------------------------- 
    237219   ln_loglayer = .true.   !  logarithmic drag: Cd = vkarmn/log(z/z0) |U| 
     220   ln_drgimp   = .true.   !  implicit top/bottom friction flag 
    238221/ 
    239222!----------------------------------------------------------------------- 
    240223&namdrg_bot    !   BOTTOM friction                                      (ln_OFF =F) 
    241224!----------------------------------------------------------------------- 
    242    rn_Cd0      =  2.5e-3   !!1.e-3    !  drag coefficient [-] 
    243    rn_Cdmax    =  0.01     !!0.1      !  drag value maximum [-] (logarithmic drag) 
    244    rn_ke0      =  0.       !!2.5e-3   !  background kinetic energy  [m2/s2] (non-linear cases) 
     225   rn_Cd0      =  2.5e-3   !  drag coefficient [-] 
     226   rn_Cdmax    =  0.1      !  drag value maximum [-] (logarithmic drag) 
     227   rn_ke0      =  0.       !  background kinetic energy  [m2/s2] (non-linear cases) 
    245228   rn_z0       =  3.e-3    !  roughness [m] (ln_loglayer=T) 
    246    ln_boost    = .false.   !  =T regional boost of Cd0 ; =F constant 
    247       rn_boost =  50.         !  local boost factor  [-] 
    248 / 
    249 !----------------------------------------------------------------------- 
    250 &nambbc        !   bottom temperature boundary condition                (default: OFF) 
    251 !----------------------------------------------------------------------- 
    252    ln_trabbc   = .false.    !  Apply a geothermal heating at the ocean bottom 
    253229/ 
    254230!----------------------------------------------------------------------- 
     
    258234      nn_bbl_ldf  =  1        !  diffusive bbl (=1)   or not (=0) 
    259235      nn_bbl_adv  =  0        !  advective bbl (=1/2) or not (=0) 
    260       rn_ahtbbl   =  1000.    !  lateral mixing coefficient in the bbl  [m2/s] 
    261       rn_gambbl   =  10.      !  advective bbl coefficient                 [s] 
    262236/ 
    263237!!====================================================================== 
     
    293267   !                       !  Coefficients: 
    294268   nn_aht_ijk_t    = 31        !  space/time variation of eddy coefficient: 
    295       !                             !   = 31 F(i,j,k,t)=F(local velocity and grid-spacing) 
     269   !                                !   = 31 F(i,j,k,t)=F(local velocity and grid-spacing) 
    296270/ 
    297271!!====================================================================== 
     
    320294&namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection) 
    321295!----------------------------------------------------------------------- 
    322    ln_dynvor_enT = .true. !  energy conserving scheme (T-point) 
     296   ln_dynvor_eeT = .true.  !  energy conserving scheme (een using e3t) 
    323297/ 
    324298!----------------------------------------------------------------------- 
     
    352326&namzdf        !   vertical physics manager                             (default: NO selection) 
    353327!----------------------------------------------------------------------- 
     328   !                       ! adaptive-implicit vertical advection 
     329   ln_zad_Aimp = .true.      !  Courant number dependent scheme (Shchepetkin 2015) 
    354330   !                       ! type of vertical closure (required) 
    355331   ln_zdftke   = .true.       !  Turbulent Kinetic Energy closure       (T =>   fill namzdf_tke) 
    356332   ln_zdfgls   = .false.      !  Generic Length Scale closure           (T =>   fill namzdf_gls) 
     333   !                       ! convection 
     334   ln_zdfevd   = .true.       !  enhanced vertical diffusion 
     335   ! 
    357336   ln_zdfddm   = .true.    ! double diffusive mixing 
     337   ! 
    358338   !                       !  Coefficients 
    359339   rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F) 
     
    384364/ 
    385365!----------------------------------------------------------------------- 
    386 &nam_diaharm   !   Harmonic analysis of tidal constituents              ('key_diaharm') 
    387 !----------------------------------------------------------------------- 
    388 / 
    389 !----------------------------------------------------------------------- 
    390366&namnc4        !   netcdf4 chunking and compression settings            ("key_netcdf4") 
    391367!----------------------------------------------------------------------- 
  • NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/namelist_ice_cfg

    r10535 r10910  
    3535&namdyn         !   Ice dynamics 
    3636!------------------------------------------------------------------------------ 
     37   ln_landfast_L16  = .true.          !  landfast: parameterization from Lemieux 2016 
    3738/ 
    3839!------------------------------------------------------------------------------ 
     
    4344&namdyn_rhg     !   Ice rheology 
    4445!------------------------------------------------------------------------------ 
     46   ln_rhg_EVP       = .true.          !  EVP rheology 
     47      ln_aEVP       = .true.          !     adaptive rheology (Kimmritz et al. 2016 & 2017) 
    4548/ 
    4649!------------------------------------------------------------------------------ 
     
    6770&namthd_do      !   Ice growth in open water 
    6871!------------------------------------------------------------------------------ 
    69    rn_hinew         =   0.02          !  thickness for new ice formation in open water (m), must be larger than rn_hnewice 
     72   rn_hinew         =   0.02          !  thickness for new ice formation in open water (m), must be larger than rn_himin 
    7073   ln_frazil        = .true.          !  Frazil ice parameterization (ice collection as a function of wind) 
    7174/ 
     
    7780&namthd_pnd     !   Melt ponds 
    7881!------------------------------------------------------------------------------ 
    79    ln_pnd_H12       = .false.         !  activate evolutive melt ponds (from Holland et al 2012) 
    80    ln_pnd_CST       = .false.         !  activate constant melt ponds 
    81    ln_pnd_alb       = .false.         !  melt ponds affect albedo or not 
     82   ln_pnd_H12       = .true.          !  activate evolutive melt ponds (from Holland et al 2012) 
     83   ln_pnd_alb       = .true.          !  melt ponds affect albedo or not 
    8284/ 
    8385 
  • NEMO/releases/release-4.0/src/ICE/icedyn.F90

    r10564 r10910  
    7474      INTEGER, INTENT(in) ::   kt     ! ice time step 
    7575      !! 
    76       INTEGER  ::   ji, jj, jl        ! dummy loop indices 
     76      INTEGER  ::   ji, jj        ! dummy loop indices 
    7777      REAL(wp) ::   zcoefu, zcoefv 
    78       REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max 
    79       REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i 
     78      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdivu_i 
    8079      !!-------------------------------------------------------------------- 
    8180      ! 
     
    8988      ENDIF 
    9089      !                       
    91       IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization 
    92          tau_icebfr(:,:) = 0._wp 
    93          DO jl = 1, jpl 
    94             WHERE( h_i_b(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
    95          END DO 
    96       ENDIF 
    97       ! 
    98       !                                !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig) 
    99       DO jl = 1, jpl 
    100          DO jj = 2, jpjm1 
    101             DO ji = fs_2, fs_jpim1 
    102                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), & 
    103                   &                                                h_ip_b(ji-1,jj  ,jl), h_ip_b(ji  ,jj-1,jl), & 
    104                   &                                                h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), & 
    105                   &                                                h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) ) 
    106                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), & 
    107                   &                                                h_i_b (ji-1,jj  ,jl), h_i_b (ji  ,jj-1,jl), & 
    108                   &                                                h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), & 
    109                   &                                                h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) ) 
    110                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), & 
    111                   &                                                h_s_b (ji-1,jj  ,jl), h_s_b (ji  ,jj-1,jl), & 
    112                   &                                                h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), & 
    113                   &                                                h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) ) 
    114             END DO 
    115          END DO 
    116       END DO 
    117       CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. ) 
    118       ! 
    119       ! 
    120       SELECT CASE( nice_dyn )           !-- Set which dynamics is running 
     90      ! retrieve thickness from volume for landfast param. and UMx advection scheme 
     91      WHERE( a_i(:,:,:) >= epsi20 ) 
     92         h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:) 
     93         h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:) 
     94      ELSEWHERE 
     95         h_i(:,:,:) = 0._wp 
     96         h_s(:,:,:) = 0._wp 
     97      END WHERE 
     98      ! 
     99      WHERE( a_ip(:,:,:) >= epsi20 ) 
     100         h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) 
     101      ELSEWHERE 
     102         h_ip(:,:,:) = 0._wp 
     103      END WHERE 
     104      ! 
     105      ! 
     106      SELECT CASE( nice_dyn )          !-- Set which dynamics is running 
    121107 
    122108      CASE ( np_dynALL )           !==  all dynamical processes  ==! 
    123          CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology   
    124          CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
    125          CALL ice_dyn_rdgrft( kt )                                                 ! -- ridging/rafting  
    126          CALL ice_cor       ( kt , 1 )                                             ! -- Corrections 
    127  
     109         ! 
     110         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology   
     111         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice 
     112         CALL ice_dyn_rdgrft( kt )                                          ! -- ridging/rafting  
     113         CALL ice_cor       ( kt , 1 )                                      ! -- Corrections 
     114         ! 
    128115      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==! 
    129          CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology   
    130          CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
    131          CALL Hpiling                                                              ! -- simple pile-up (replaces ridging/rafting) 
    132  
     116         ! 
     117         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology   
     118         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice 
     119         CALL Hpiling                                                       ! -- simple pile-up (replaces ridging/rafting) 
     120         ! 
    133121      CASE ( np_dynADV1D )         !==  pure advection ==!   (1D) 
    134          ALLOCATE( zdivu_i(jpi,jpj) ) 
     122         ! 
    135123         ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- ! 
    136124         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length 
     
    145133         END DO 
    146134         ! --- 
    147          CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice 
    148  
    149          ! diagnostics: divergence at T points  
    150          DO jj = 2, jpjm1 
    151             DO ji = 2, jpim1 
    152                zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    153                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
    154             END DO 
    155          END DO 
    156          CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 
    157          IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) ) 
    158  
    159          DEALLOCATE( zdivu_i ) 
    160  
     135         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice 
     136         ! 
    161137      CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities) 
    162          ALLOCATE( zdivu_i(jpi,jpj) ) 
     138         ! 
    163139         u_ice(:,:) = rn_uice * umask(:,:,1) 
    164140         v_ice(:,:) = rn_vice * vmask(:,:,1) 
     
    166142         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 
    167143         ! --- 
    168          CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice 
    169  
    170          ! diagnostics: divergence at T points  
    171          DO jj = 2, jpjm1 
    172             DO ji = 2, jpim1 
    173                zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    174                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     144         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice 
     145 
     146      END SELECT 
     147      ! 
     148      ! 
     149      ! diagnostics: divergence at T points  
     150      IF( iom_use('icediv') ) THEN 
     151         ! 
     152         SELECT CASE( nice_dyn ) 
     153 
     154         CASE ( np_dynADV1D , np_dynADV2D ) 
     155 
     156            ALLOCATE( zdivu_i(jpi,jpj) ) 
     157            DO jj = 2, jpjm1 
     158               DO ji = 2, jpim1 
     159                  zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     160                     &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     161               END DO 
    175162            END DO 
    176          END DO 
    177          CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 
    178          IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) ) 
    179  
    180          DEALLOCATE( zdivu_i ) 
    181  
    182       END SELECT 
     163            CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 
     164            CALL iom_put( "icediv" , zdivu_i(:,:) ) 
     165            DEALLOCATE( zdivu_i ) 
     166 
     167         END SELECT 
     168         ! 
     169      ENDIF 
    183170      ! 
    184171      ! controls 
     
    186173      ! 
    187174   END SUBROUTINE ice_dyn 
    188  
    189  
    190    SUBROUTINE Hbig( phi_max, phs_max, phip_max ) 
    191       !!------------------------------------------------------------------- 
    192       !!                  ***  ROUTINE Hbig  *** 
    193       !! 
    194       !! ** Purpose : Thickness correction in case advection scheme creates 
    195       !!              abnormally tick ice or snow 
    196       !! 
    197       !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points 
    198       !!                 (before advection) and reduce it by adapting ice concentration 
    199       !!              2- check whether snow thickness is larger than the surrounding 9-points 
    200       !!                 (before advection) and reduce it by sending the excess in the ocean 
    201       !!              3- check whether snow load deplets the snow-ice interface below sea level$ 
    202       !!                 and reduce it by sending the excess in the ocean 
    203       !!              4- correct pond fraction to avoid a_ip > a_i 
    204       !! 
    205       !! ** input   : Max thickness of the surrounding 9-points 
    206       !!------------------------------------------------------------------- 
    207       REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    208       ! 
    209       INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    210       REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra 
    211       !!------------------------------------------------------------------- 
    212       ! controls 
    213       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    214       ! 
    215       CALL ice_var_zapsmall                       !-- zap small areas 
    216       ! 
    217       DO jl = 1, jpl 
    218          DO jj = 1, jpj 
    219             DO ji = 1, jpi 
    220                IF ( v_i(ji,jj,jl) > 0._wp ) THEN 
    221                   ! 
    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 
    229                   ENDIF 
    230                   ! 
    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 = phs_max(ji,jj,jl) / MAX( zhs, 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 = ( v_s(ji,jj,jl) - 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                    
    265                ENDIF 
    266             END DO 
    267          END DO 
    268       END DO  
    269       !                                           !-- correct pond fraction to avoid a_ip > a_i 
    270       WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:) 
    271       ! 
    272       ! controls 
    273       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    274       ! 
    275    END SUBROUTINE Hbig 
    276175 
    277176 
     
    337236         WRITE(numout,*) '~~~~~~~~~~~~' 
    338237         WRITE(numout,*) '   Namelist namdyn:' 
    339          WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynALL       = ', ln_dynALL 
    340          WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV    = ', ln_dynRHGADV 
    341          WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)     ln_dynADV1D     = ', ln_dynADV1D 
    342          WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                 ln_dynADV2D     = ', ln_dynADV2D 
    343          WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)    = (', rn_uice,',', rn_vice,')' 
    344          WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat       = ', rn_ishlat 
    345          WRITE(numout,*) '      Landfast: param from Lemieux 2016                       ln_landfast_L16 = ', ln_landfast_L16 
    346          WRITE(numout,*) '      Landfast: param from home made                          ln_landfast_home= ', ln_landfast_home 
    347          WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_depfra       = ', rn_depfra 
    348          WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr       = ', rn_icebfr 
    349          WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax      = ', rn_lfrelax 
    350          WRITE(numout,*) '         isotropic tensile strength                           rn_tensile      = ', rn_tensile 
     238         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr) ln_dynALL       = ', ln_dynALL 
     239         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                     ln_dynRHGADV    = ', ln_dynRHGADV 
     240         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)    ln_dynADV1D     = ', ln_dynADV1D 
     241         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                ln_dynADV2D     = ', ln_dynADV2D 
     242         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)   = (', rn_uice,',',rn_vice,')' 
     243         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat 
     244         WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16 
     245         WRITE(numout,*) '      Landfast: param from home made                         ln_landfast_home= ', ln_landfast_home 
     246         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra 
     247         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr 
     248         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax      = ', rn_lfrelax 
     249         WRITE(numout,*) '         isotropic tensile strength                          rn_tensile      = ', rn_tensile 
    351250         WRITE(numout,*) 
    352251      ENDIF 
  • NEMO/releases/release-4.0/src/ICE/icedyn_adv.F90

    r10413 r10910  
    6464      !!---------------------------------------------------------------------- 
    6565      INTEGER, INTENT(in) ::   kt   ! number of iteration 
    66       ! 
    67       INTEGER ::   jl   ! dummy loop indice 
    68       REAL(wp), DIMENSION(jpi,jpj) ::   zmask  ! fraction of time step with v_i < 0 
    6966      !!--------------------------------------------------------------------- 
    7067      ! 
    71       IF( ln_timing )   CALL timing_start('icedyn_adv') 
     68      ! controls 
     69      IF( ln_timing    )   CALL timing_start('icedyn_adv')                                                             ! timing 
     70      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    7271      ! 
    7372      IF( kt == nit000 .AND. lwp ) THEN 
     
    7675         WRITE(numout,*) '~~~~~~~~~~~' 
    7776      ENDIF 
    78        
    79       ! conservation test 
    80       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    81                       
    82       !---------- 
    83       ! Advection 
    84       !---------- 
     77      ! 
     78      !---------------! 
     79      !== Advection ==! 
     80      !---------------! 
    8581      SELECT CASE( nice_adv ) 
    8682      !                                !-----------------------! 
    8783      CASE( np_advUMx )                ! ULTIMATE-MACHO scheme ! 
    8884         !                             !-----------------------! 
    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 ) 
    90       !                                !-----------------------! 
     85         CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, & 
     86            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     87         !                             !-----------------------! 
    9188      CASE( np_advPRA )                ! PRATHER scheme        ! 
    9289         !                             !-----------------------! 
    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 ) 
     90         CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, & 
     91            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    9492      END SELECT 
    95  
    96       !---------------------------- 
    97       ! Debug the advection schemes 
    98       !---------------------------- 
    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. 
    102       !                   These values are usually very small but in some occasions they can also be non-negligible 
    103       !                   Therefore one needs to bound the advected fields by 0 (though this is not a clean fix) 
    104       ! 
    105       ! record the negative values resulting from UMx 
    106       zmask(:,:) = 0._wp ! keep the init to 0 here 
    107       DO jl = 1, jpl 
    108          WHERE( v_i(:,:,jl) < 0._wp )   zmask(:,:) = 1._wp 
    109       END DO 
    110       IF( iom_use('iceneg_pres') )   CALL iom_put("iceneg_pres", zmask                                      )  ! fraction of time step with v_i < 0 
    111       IF( iom_use('iceneg_volu') )   CALL iom_put("iceneg_volu", SUM(MIN( v_i, 0. ), dim=3 )                )  ! negative ice volume (only) 
    112       IF( iom_use('iceneg_hfx' ) )   CALL iom_put("iceneg_hfx" , ( SUM(SUM( MIN( e_i(:,:,1:nlay_i,:), 0. )  &  ! negative ice heat content (only) 
    113          &                                                                  , dim=4 ), dim=3 ) )* r1_rdtice )  ! -- eq. heat flux [W/m2] 
    114       ! 
    115       ! ==> conservation is ensured by calling this routine below, 
    116       !     however the global ice volume is then changed by advection (but errors are small)  
    117       CALL ice_var_zapneg( ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    11893 
    11994      !------------ 
  • NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90

    r10785 r10910  
    1111   !!   'key_si3'                                       SI3 sea-ice model 
    1212   !!---------------------------------------------------------------------- 
    13    !!   ice_dyn_adv_umx   : update the tracer trend with the 3D advection trends using a TVD scheme 
     13   !!   ice_dyn_adv_umx   : update the tracer fields 
    1414   !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders 
    15    !!   macho             : ??? 
    16    !!   nonosc_ice        : compute monotonic tracer fluxes by a non-oscillatory algorithm  
     15   !!   macho             : compute the fluxes 
     16   !!   nonosc_ice        : limit the fluxes using a non-oscillatory algorithm  
    1717   !!---------------------------------------------------------------------- 
    1818   USE phycst         ! physical constant 
     
    3939   INTEGER ::   kn_limiter = 1 
    4040 
    41    ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme 
    42    !   clem: if set to true, the 2D test case "diagonal advection" does not work (I do not understand why) 
    43    !         but in realistic cases, it avoids having very negative ice temperature (-50) at low ice concentration  
     41   ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge), 
     42   !    interpolated T at u/v points can be non-zero while it should 
     43   !    (because of the high order of the advection scheme). Thus set it to 0 in this case 
     44   LOGICAL ::   ll_icedge = .TRUE. 
     45 
     46   ! if T interpolated at u/v points is negative or v_i < 1.e-6 
     47   !    then interpolate T at u/v points using the upstream scheme 
    4448   LOGICAL ::   ll_neg = .TRUE. 
    4549    
     
    5155    
    5256   ! prelimiter: use it to avoid overshoot in H 
    53    !   clem: if set to true, the 2D test case "diagnoal advection" does not work (I do not understand why) 
    5457   LOGICAL ::   ll_prelimiter_zalesak = .FALSE.  ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D 
    5558 
     59   ! advection for S and T:    dVS/dt = -div( uA * uHS / u ) => ll_advS = F 
     60   !                        or dVS/dt = -div( uV * uS  / u ) => ll_advS = T 
     61   LOGICAL ::   ll_advS = .FALSE. 
    5662 
    5763   !! * Substitutions 
     
    6470CONTAINS 
    6571 
    66    SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, & 
     72   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, & 
    6773      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    6874      !!---------------------------------------------------------------------- 
     
    7985      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity 
    8086      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity 
     87      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_i       ! ice thickness 
     88      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_s       ! snw thickness 
     89      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_ip      ! ice pond thickness 
    8190      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area 
    8291      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume 
     
    94103      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
    95104      REAL(wp) ::   zdt 
    96       REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! send zcflnow and receive zcflprv 
     105      REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! for global communication 
    97106      REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box  
    98107      REAL(wp), DIMENSION(jpi,jpj)     ::   zati1, zati2 
    99       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho 
    100       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip 
     108      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho, zua_ups, zva_ups 
     109      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups  
     110      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip, z1_vi, z1_vs  
    101111      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhvar 
     112      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max 
    102113      !!---------------------------------------------------------------------- 
    103114      ! 
    104115      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 
    105116      ! 
    106       ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! 
    107       !     When needed, the advection split is applied at the next time-step in order to avoid blocking global comm. 
    108       !     ...this should not affect too much the stability... Was ok on the tests we did... 
     117      ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 
     118      DO jl = 1, jpl 
     119         DO jj = 2, jpjm1 
     120            DO ji = fs_2, fs_jpim1 
     121               zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
     122                  &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
     123                  &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
     124                  &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
     125               zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
     126                  &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
     127                  &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
     128                  &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
     129               zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
     130                  &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
     131                  &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
     132                  &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
     133            END DO 
     134         END DO 
     135      END DO 
     136      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. ) 
     137      ! 
     138      ! 
     139      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- ! 
     140      !        Note: the advection split is applied at the next time-step in order to avoid blocking global comm. 
     141      !              this should not affect too much the stability 
    109142      zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 
    110143      zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 
     
    116149      ELSE                         ;   icycle = 1 
    117150      ENDIF 
    118        
    119151      zdt = rdt_ice / REAL(icycle) 
    120152 
     
    154186         END WHERE 
    155187         ! 
    156          ! set u*a=u for advection of A only  
     188         ! set u*A=u for advection of A only  
    157189         DO jl = 1, jpl 
    158190            zua_ho(:,:,jl) = zudy(:,:) 
    159191            zva_ho(:,:,jl) = zvdx(:,:) 
    160192         END DO 
    161           
     193 
     194         !== Ice area ==! 
    162195         zamsk = 1._wp 
    163          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 
     196         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     197            &                                      pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho ) 
    164198         zamsk = 0._wp 
    165          ! 
    166          zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
    167          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_i )                !-- Ice volume 
    168          ! 
    169          zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
    170          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_s )                !-- Snw volume 
    171          ! 
    172          zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 
    173          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, psv_i )               !-- Salt content 
    174          ! 
    175          DO jk = 1, nlay_i 
    176             zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 
    177             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 
    178          END DO 
    179          ! 
    180          DO jk = 1, nlay_s 
    181             zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 
    182             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 
    183          END DO 
    184          ! 
    185          IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN                                                              !-- Ice Age 
     199 
     200         IF( .NOT. ll_advS ) THEN   !-- advection form: uA * uHS / u 
     201            !== Ice volume ==! 
     202            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
     203            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     204               &                                      zhvar, pv_i, zua_ups, zva_ups ) 
     205            ! 
     206            !== Snw volume ==!          
     207            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
     208            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     209               &                                      zhvar, pv_s, zua_ups, zva_ups ) 
     210            ! 
     211            !== Salt content ==! 
     212            zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 
     213            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     214               &                                      zhvar, psv_i, zua_ups, zva_ups ) 
     215            ! 
     216            !== Ice heat content ==! 
     217            DO jk = 1, nlay_i 
     218               zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 
     219               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, & 
     220                  &                                      zhvar, pe_i(:,:,jk,:), zua_ups, zva_ups ) 
     221            END DO 
     222            ! 
     223            !== Snw heat content ==! 
     224            DO jk = 1, nlay_s 
     225               zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 
     226               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, & 
     227                  &                                      zhvar, pe_s(:,:,jk,:), zua_ups, zva_ups ) 
     228            END DO 
     229            ! 
     230         ELSE                       !-- advection form: uV * uS / u 
     231            ! inverse of Vi 
     232            WHERE( pv_i(:,:,:) >= epsi20 )   ;   z1_vi(:,:,:) = 1._wp / pv_i(:,:,:) 
     233            ELSEWHERE                        ;   z1_vi(:,:,:) = 0. 
     234            END WHERE 
     235            ! inverse of Vs 
     236            WHERE( pv_s(:,:,:) >= epsi20 )   ;   z1_vs(:,:,:) = 1._wp / pv_s(:,:,:) 
     237            ELSEWHERE                        ;   z1_vs(:,:,:) = 0. 
     238            END WHERE 
     239            ! 
     240            ! It is important to first calculate the ice fields and then the snow fields (because we use the same arrays) 
     241            ! 
     242            !== Ice volume ==! 
     243            zuv_ups = zua_ups 
     244            zvv_ups = zva_ups 
     245            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
     246            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     247               &                                      zhvar, pv_i, zuv_ups, zvv_ups, zuv_ho , zvv_ho ) 
     248            ! 
     249            !== Salt content ==! 
     250            zhvar(:,:,:) = psv_i(:,:,:) * z1_vi(:,:,:) 
     251            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, & 
     252               &                                      zhvar, psv_i, zuv_ups, zvv_ups ) 
     253            ! 
     254            !== Ice heat content ==! 
     255            DO jk = 1, nlay_i 
     256               zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_vi(:,:,:) 
     257               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, & 
     258                  &                                      zhvar, pe_i(:,:,jk,:), zuv_ups, zvv_ups ) 
     259            END DO 
     260            ! 
     261            !== Snow volume ==!          
     262            zuv_ups = zua_ups 
     263            zvv_ups = zva_ups 
     264            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
     265            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     266               &                                      zhvar, pv_s, zuv_ups, zvv_ups, zuv_ho , zvv_ho ) 
     267            ! 
     268            !== Snw heat content ==! 
     269            DO jk = 1, nlay_s 
     270               zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_vs(:,:,:) 
     271               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, & 
     272                  &                                      zhvar, pe_s(:,:,jk,:), zuv_ups, zvv_ups ) 
     273            END DO 
     274            ! 
     275         ENDIF 
     276         ! 
     277         !== Ice age ==! 
     278         IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN 
    186279            ! clem: in theory we should use the formulation below to advect the ice age, but the code is unable to deal with 
    187280            !       fields that do not depend on volume (here oa_i depends on concentration). It creates abnormal ages that 
     
    189282            !!zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:) 
    190283            !!CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i ) 
    191             ! set u*a=u for advection of ice age 
     284            ! set u*A=u for advection of ice age 
    192285            DO jl = 1, jpl 
    193286               zua_ho(:,:,jl) = zudy(:,:) 
     
    195288            END DO 
    196289            zamsk = 1._wp 
    197             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, poa_i, poa_i ) 
     290            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho, zva_ho, zcu_box, zcv_box, & 
     291               &                                      poa_i, poa_i ) 
    198292            zamsk = 0._wp 
    199293         ENDIF 
    200294         ! 
    201          IF ( ln_pnd_H12 ) THEN                                                                                               !-- melt ponds 
    202             ! set u*a=u for advection of Ap only  
     295         !== melt ponds ==! 
     296         IF ( ln_pnd_H12 ) THEN 
     297            ! set u*A=u for advection of Ap only  
    203298            DO jl = 1, jpl 
    204299               zua_ho(:,:,jl) = zudy(:,:) 
    205300               zva_ho(:,:,jl) = zvdx(:,:) 
    206301            END DO 
    207             ! 
     302            ! fraction 
    208303            zamsk = 1._wp 
    209             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 ) ! fraction 
     304            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     305               &                                      pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho ) 
    210306            zamsk = 0._wp 
    211             ! 
     307            ! volume 
    212308            zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:) 
    213             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_ip )                 ! volume 
     309            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     310               &                                      zhvar, pv_ip, zua_ups, zva_ups ) 
    214311         ENDIF 
    215312         ! 
     313         !== Open water area ==! 
    216314         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
    217315         DO jj = 2, jpjm1 
    218316            DO ji = fs_2, fs_jpim1 
    219                pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                                              !-- Open water area 
     317               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &  
    220318                  &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
    221319            END DO 
     
    223321         CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T',  1. ) 
    224322         ! 
     323         ! 
     324         ! --- Ensure non-negative fields and in-bound thicknesses --- ! 
     325         ! 
     326         ! Remove negative values (conservation is ensured) 
     327         !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
     328         CALL ice_var_zapneg( pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     329         ! 
     330         ! Make sure ice thickness is not too big 
     331         !    (because ice thickness can be too large where ice concentration is very small) 
     332         CALL Hbig( zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     333 
    225334      END DO 
    226335      ! 
     
    228337 
    229338    
    230    SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho ) 
     339   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox,  & 
     340      &                                            pt, ptc, pua_ups, pva_ups, pua_ho, pva_ho ) 
    231341      !!---------------------------------------------------------------------- 
    232342      !!                  ***  ROUTINE adv_umx  *** 
     
    235345      !!                 tracers and add it to the general trend of tracer equations 
    236346      !! 
    237       !! **  Method  :   - calculate upstream fluxes and upstream solution for tracer H 
     347      !! **  Method  :   - calculate upstream fluxes and upstream solution for tracers V/A(=H) etc 
    238348      !!                 - calculate tracer H at u and v points (Ultimate) 
    239       !!                 - calculate the high order fluxes using alterning directions (Macho?) 
     349      !!                 - calculate the high order fluxes using alterning directions (Macho) 
    240350      !!                 - apply a limiter on the fluxes (nonosc_ice) 
    241       !!                 - convert this tracer flux to a tracer content flux (uH -> uV) 
    242       !!                 - calculate the high order solution for tracer content V 
     351      !!                 - convert this tracer flux to a "volume" flux (uH -> uV) 
     352      !!                 - apply a limiter a second time on the volumes fluxes (nonosc_ice) 
     353      !!                 - calculate the high order solution for V 
    243354      !! 
    244       !! ** Action : solve 2 equations => a) da/dt = -div(ua) 
    245       !!                                  b) dV/dt = -div(uV) using dH/dt = -u.grad(H) 
    246       !!             in eq. b), - fluxes uH are evaluated (with UMx) and limited (with nonosc_ice). This step is necessary to get a good H. 
    247       !!                        - then we convert this flux to a "volume" flux this way => uH*ua/u 
    248       !!                             where ua is the flux from eq. a) 
    249       !!                        - at last we estimate dV/dt = -div(uH*ua/u) 
     355      !! ** Action : solve 3 equations => a) dA/dt  = -div(uA) 
     356      !!                                  b) dV/dt  = -div(uV)  using dH/dt = -u.grad(H) 
     357      !!                                  c) dVS/dt = -div(uVS) using either dHS/dt = -u.grad(HS) or dS/dt = -u.grad(S) 
    250358      !! 
    251       !! ** Note : - this method can lead to small negative V (since we only limit H) => corrected in icedyn_adv.F90 conserving mass etc. 
    252       !!           - negative tracers at u-v points can also occur from the Ultimate scheme (usually at the ice edge) and the solution for now 
    253       !!             is to apply an upstream scheme when it occurs. A better solution would be to degrade the order of 
    254       !!             the scheme automatically by applying a mask of the ice cover inside Ultimate (not done). 
     359      !!             in eq. b), - fluxes uH are evaluated (with UMx) and limited with nonosc_ice. This step is necessary to get a good H. 
     360      !!                        - then we convert this flux to a "volume" flux this way => uH * uA / u 
     361      !!                             where uA is the flux from eq. a) 
     362      !!                             this "volume" flux is also limited with nonosc_ice (otherwise overshoots can occur) 
     363      !!                        - at last we estimate dV/dt = -div(uH * uA / u) 
     364      !! 
     365      !!             in eq. c), one can solve the equation for  S (ln_advS=T), then dVS/dt = -div(uV * uS  / u) 
     366      !!                                                or for HS (ln_advS=F), then dVS/dt = -div(uA * uHS / u)  
     367      !! 
     368      !! ** Note : - this method can lead to tiny negative V (-1.e-20) => set it to 0 while conserving mass etc. 
     369      !!           - At the ice edge, Ultimate scheme can lead to: 
     370      !!                              1) negative interpolated tracers at u-v points 
     371      !!                              2) non-zero interpolated tracers at u-v points eventhough there is no ice and velocity is outward 
     372      !!                              Solution for 1): apply an upstream scheme when it occurs. A better solution would be to degrade the order of 
     373      !!                                               the scheme automatically by applying a mask of the ice cover inside Ultimate (not done). 
     374      !!                              Solution for 2): we set it to 0 in this case 
    255375      !!           - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good. 
    256376      !!             Large values of H can appear for very small ice concentration, and when it does it messes the things up since we 
    257       !!             work on H (and not V). It probably comes from the prelimiter of zalesak which is coded for 1D and not 2D. 
     377      !!             work on H (and not V). It is partly related to the multi-category approach 
    258378      !!             Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 
    259       !!             concentration is small). 
    260       !! To-do: expand the prelimiter from zalesak to make it work in 2D 
    261       !!---------------------------------------------------------------------- 
    262       REAL(wp)                        , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0) 
    263       INTEGER                         , INTENT(in   )           ::   kn_umx         ! order of the scheme (1-5=UM or 20=CEN2) 
    264       INTEGER                         , INTENT(in   )           ::   jt             ! number of sub-iteration 
    265       INTEGER                         , INTENT(in   )           ::   kt             ! number of iteration 
    266       REAL(wp)                        , INTENT(in   )           ::   pdt            ! tracer time-step 
    267       REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2 
    268       REAL(wp), DIMENSION(:,:,:)      , INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u 
    269       REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pubox, pvbox   ! upstream velocity 
    270       REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   pt             ! tracer field 
    271       REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   ptc            ! tracer content field 
    272       REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes 
     379      !!             concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 
     380      !!             since sv_i and e_i are still good. 
     381      !!---------------------------------------------------------------------- 
     382      REAL(wp)                        , INTENT(in   )           ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     383      INTEGER                         , INTENT(in   )           ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2) 
     384      INTEGER                         , INTENT(in   )           ::   jt               ! number of sub-iteration 
     385      INTEGER                         , INTENT(in   )           ::   kt               ! number of iteration 
     386      REAL(wp)                        , INTENT(in   )           ::   pdt              ! tracer time-step 
     387      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pu   , pv        ! 2 ice velocity components => u*e2 
     388      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   )           ::   puc  , pvc       ! 2 ice velocity components => u*e2 or u*a*e2u 
     389      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pubox, pvbox     ! upstream velocity 
     390      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   pt               ! tracer field 
     391      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   ptc              ! tracer content field 
     392      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout), OPTIONAL ::   pua_ups, pva_ups ! upstream u*a fluxes 
     393      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho   ! high order u*a fluxes 
    273394      ! 
    274395      INTEGER  ::   ji, jj, jl       ! dummy loop indices   
     
    303424            DO jj = 1, jpjm1 
    304425               DO ji = 1, fs_jpim1 
    305                   IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN 
    306                      zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 
    307                      zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 
     426                  IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 
     427                     zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj) 
     428                     zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj) 
    308429                  ELSE 
    309430                     zfu_ho (ji,jj,jl) = 0._wp 
     
    311432                  ENDIF 
    312433                  ! 
    313                   IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN 
    314                      zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 
    315                      zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 
     434                  IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 
     435                     zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj) 
     436                     zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj) 
    316437                  ELSE 
    317438                     zfv_ho (ji,jj,jl) = 0._wp   
     
    321442            END DO 
    322443         END DO 
     444 
     445         ! the new "volume" fluxes must also be "flux corrected" 
     446         ! thus we calculate the upstream solution and apply a limiter again 
     447         DO jl = 1, jpl 
     448            DO jj = 2, jpjm1 
     449               DO ji = fs_2, fs_jpim1 
     450                  ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) 
     451                  ! 
     452                  zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 
     453               END DO 
     454            END DO 
     455         END DO 
     456         CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T',  1. ) 
     457         ! 
     458         IF    ( kn_limiter == 1 ) THEN 
     459            CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho ) 
     460         ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN 
     461            CALL limiter_x( pdt, pu, ptc, zfu_ups, zfu_ho ) 
     462            CALL limiter_y( pdt, pv, ptc, zfv_ups, zfv_ho ) 
     463         ENDIF 
     464         ! 
    323465      ENDIF 
    324       !                                   --ho 
    325       ! in case of advection of A: output u*a 
    326       ! ------------------------------------- 
     466      !                                   --ho    --ups 
     467      ! in case of advection of A: output u*a and u*a 
     468      ! ----------------------------------------------- 
    327469      IF( PRESENT( pua_ho ) ) THEN 
    328470         DO jl = 1, jpl 
    329471            DO jj = 1, jpjm1 
    330472               DO ji = 1, fs_jpim1 
    331                   pua_ho(ji,jj,jl) = zfu_ho(ji,jj,jl) 
    332                   pva_ho(ji,jj,jl) = zfv_ho(ji,jj,jl) 
    333                END DO 
     473                  pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 
     474                  pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 
     475              END DO 
    334476            END DO 
    335477         END DO 
     
    609751         ! 
    610752         !                                                        !--  ultimate interpolation of pt at u-point  --! 
    611          CALL ultimate_x( kn_umx, pdt, pt, pu, zt_u, pfu_ho ) 
     753         CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho ) 
    612754         !                                                        !--  limiter in x --! 
    613755         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     
    619761                     &                              + pt   (ji,jj,jl) * ( pu  (ji,jj   ) - pu  (ji-1,jj   ) ) * r1_e1e2t(ji,jj) & 
    620762                     &                                                                                        * pamsk           & 
    621                      &                             ) * pdt ) * tmask(ji,jj,1)  
     763                     &                             ) * pdt ) * tmask(ji,jj,1) 
     764                  !!clem test 
     765                  !!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) ) 
     766                  !!clem test 
    622767               END DO 
    623768            END DO 
     
    627772         !                                                        !--  ultimate interpolation of pt at v-point  --! 
    628773         IF( ll_hoxy ) THEN 
    629             CALL ultimate_y( kn_umx, pdt, zpt, pv, zt_v, pfv_ho ) 
     774            CALL ultimate_y( pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho ) 
    630775         ELSE 
    631             CALL ultimate_y( kn_umx, pdt, pt , pv, zt_v, pfv_ho ) 
     776            CALL ultimate_y( pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho ) 
    632777         ENDIF 
    633778         !                                                        !--  limiter in y --! 
     
    638783         ! 
    639784         !                                                        !--  ultimate interpolation of pt at v-point  --! 
    640          CALL ultimate_y( kn_umx, pdt, pt, pv, zt_v, pfv_ho ) 
     785         CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho ) 
    641786         !                                                        !--  limiter in y --! 
    642787         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     
    649794                     &                                                                                        * pamsk           & 
    650795                     &                             ) * pdt ) * tmask(ji,jj,1)  
     796                  !!clem test 
     797                  !!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) ) 
     798                  !!clem test 
    651799               END DO 
    652800            END DO 
     
    656804         !                                                        !--  ultimate interpolation of pt at u-point  --! 
    657805         IF( ll_hoxy ) THEN 
    658             CALL ultimate_x( kn_umx, pdt, zpt, pu, zt_u, pfu_ho ) 
     806            CALL ultimate_x( pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho ) 
    659807         ELSE 
    660             CALL ultimate_x( kn_umx, pdt, pt , pu, zt_u, pfu_ho ) 
     808            CALL ultimate_x( pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho ) 
    661809         ENDIF 
    662810         !                                                        !--  limiter in x --! 
     
    670818 
    671819 
    672    SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 
     820   SUBROUTINE ultimate_x( pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 
    673821      !!--------------------------------------------------------------------- 
    674822      !!                    ***  ROUTINE ultimate_x  *** 
     
    680828      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    681829      !!---------------------------------------------------------------------- 
     830      REAL(wp)                        , INTENT(in   ) ::   pamsk     ! advection of concentration (1) or other tracers (0) 
    682831      INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
    683832      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step 
     
    688837      ! 
    689838      INTEGER  ::   ji, jj, jl             ! dummy loop indices 
    690       REAL(wp) ::   zcu, zdx2, zdx4        !   -      - 
     839      REAL(wp) ::   zcu, zdx2, zdx4, zvi_cen2        !   -      - 
    691840      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztu1, ztu2, ztu3, ztu4 
    692841      !!---------------------------------------------------------------------- 
     
    799948      END SELECT 
    800949      ! 
     950      ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge), 
     951      !    interpolated T at u/v points can be non-zero while it should 
     952      !    (because of the high order of the advection scheme). Thus set it to 0 in this case 
     953      IF( ll_icedge ) THEN 
     954         DO jl = 1, jpl 
     955            DO jj = 1, jpjm1 
     956               DO ji = 1, fs_jpim1 
     957                  IF( pt(ji,jj,jl) <= 0._wp .AND. pu(ji,jj) >= 0._wp ) THEN 
     958                     pt_u(ji,jj,jl) = 0._wp 
     959                  ENDIF 
     960               END DO 
     961            END DO 
     962         END DO 
     963      ENDIF 
     964      ! 
    801965      ! if pt at u-point is negative then use the upstream value 
    802966      !    this should not be necessary if a proper sea-ice mask is set in Ultimate 
     
    806970            DO jj = 1, jpjm1 
    807971               DO ji = 1, fs_jpim1 
    808                   IF( pt_u(ji,jj,jl) < 0._wp ) THEN 
     972                  zvi_cen2 = 0.5_wp * ( v_i(ji+1,jj,jl) + v_i(ji,jj,jl) ) 
     973                  IF( pt_u(ji,jj,jl) < 0._wp .OR. ( zvi_cen2 < epsi06 .AND. pamsk == 0._wp ) ) THEN 
    809974                     pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    810975                        &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
     
    826991    
    827992  
    828    SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 
     993   SUBROUTINE ultimate_y( pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 
    829994      !!--------------------------------------------------------------------- 
    830995      !!                    ***  ROUTINE ultimate_y  *** 
     
    8361001      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    8371002      !!---------------------------------------------------------------------- 
     1003      REAL(wp)                        , INTENT(in   ) ::   pamsk     ! advection of concentration (1) or other tracers (0) 
    8381004      INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
    8391005      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step 
     
    8441010      ! 
    8451011      INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    846       REAL(wp) ::   zcv, zdy2, zdy4    !   -      - 
     1012      REAL(wp) ::   zcv, zdy2, zdy4, zvi_cen2    !   -      - 
    8471013      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztv1, ztv2, ztv3, ztv4 
    8481014      !!---------------------------------------------------------------------- 
     
    9521118      END SELECT 
    9531119      ! 
     1120      ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge), 
     1121      !    interpolated T at u/v points can be non-zero while it should 
     1122      !    (because of the high order of the advection scheme). Thus set it to 0 in this case 
     1123      IF( ll_icedge ) THEN 
     1124         DO jl = 1, jpl 
     1125            DO jj = 1, jpjm1 
     1126               DO ji = 1, fs_jpim1 
     1127                  IF( pt(ji,jj,jl) <= 0._wp .AND. pv(ji,jj) >= 0._wp ) THEN 
     1128                     pt_v(ji,jj,jl) = 0._wp 
     1129                  ENDIF 
     1130               END DO 
     1131            END DO 
     1132         END DO 
     1133      ENDIF 
     1134      ! 
    9541135      ! if pt at v-point is negative then use the upstream value 
    9551136      !    this should not be necessary if a proper sea-ice mask is set in Ultimate 
     
    9591140            DO jj = 1, jpjm1 
    9601141               DO ji = 1, fs_jpim1 
    961                   IF( pt_v(ji,jj,jl) < 0._wp ) THEN 
     1142                  zvi_cen2 = 0.5_wp * ( v_i(ji,jj+1,jl) + v_i(ji,jj,jl) ) 
     1143                  IF( pt_v(ji,jj,jl) < 0._wp .OR. ( zvi_cen2 < epsi06 .AND. pamsk == 0._wp ) ) THEN 
    9621144                     pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
    9631145                        &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
     
    11021284               ! 
    11031285               !                                  ! up & down beta terms 
    1104                IF( zpos > 0._wp ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 
    1105                ELSE                    ; zbetup(ji,jj,jl) = 0._wp ! zbig 
     1286               ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!) 
     1287               IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 
     1288               ELSE                     ; zbetup(ji,jj,jl) = 0._wp ! zbig 
    11061289               ENDIF 
    11071290               ! 
    1108                IF( zneg > 0._wp ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 
    1109                ELSE                    ; zbetdo(ji,jj,jl) = 0._wp ! zbig 
     1291               IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 
     1292               ELSE                     ; zbetdo(ji,jj,jl) = 0._wp ! zbig 
    11101293               ENDIF 
    11111294               ! 
     
    11481331            END DO 
    11491332         END DO 
    1150  
    1151          ! clem test 
    1152 !!         DO jj = 2, jpjm1 
    1153 !!            DO ji = 2, fs_jpim1   ! vector opt. 
    1154 !!               zzt = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
    1155 !!                  &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
    1156 !!                  &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
    1157 !!                  &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
    1158 !!                  &         ) * tmask(ji,jj,1) 
    1159 !!               IF( zzt < -epsi20 ) THEN 
    1160 !!                  WRITE(numout,*) 'T<0 nonosc_ice',zzt 
    1161 !!               ENDIF 
    1162 !!            END DO 
    1163 !!         END DO 
    11641333 
    11651334      END DO 
     
    13581527   END SUBROUTINE limiter_y 
    13591528 
     1529 
     1530   SUBROUTINE Hbig( phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     1531      !!------------------------------------------------------------------- 
     1532      !!                  ***  ROUTINE Hbig  *** 
     1533      !! 
     1534      !! ** Purpose : Thickness correction in case advection scheme creates 
     1535      !!              abnormally tick ice or snow 
     1536      !! 
     1537      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points 
     1538      !!                 (before advection) and reduce it by adapting ice concentration 
     1539      !!              2- check whether snow thickness is larger than the surrounding 9-points 
     1540      !!                 (before advection) and reduce it by sending the excess in the ocean 
     1541      !!              3- check whether snow load deplets the snow-ice interface below sea level$ 
     1542      !!                 and reduce it by sending the excess in the ocean 
     1543      !!              4- correct pond fraction to avoid a_ip > a_i 
     1544      !! 
     1545      !! ** input   : Max thickness of the surrounding 9-points 
     1546      !!------------------------------------------------------------------- 
     1547      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
     1548      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip 
     1549      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
     1550      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i 
     1551      ! 
     1552      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
     1553      REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra 
     1554      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch 
     1555      !!------------------------------------------------------------------- 
     1556      ! 
     1557      ! 
     1558      DO jl = 1, jpl 
     1559 
     1560         DO jj = 1, jpj 
     1561            DO ji = 1, jpi 
     1562               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     1563                  ! 
     1564                  !                               ! -- check h_ip -- ! 
     1565                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
     1566                  IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     1567                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
     1568                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
     1569                        pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
     1570                     ENDIF 
     1571                  ENDIF 
     1572                  ! 
     1573                  !                               ! -- check h_i -- ! 
     1574                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
     1575                  zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 
     1576                  IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1577                     pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
     1578                  ENDIF 
     1579                  ! 
     1580                  !                               ! -- check h_s -- ! 
     1581                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
     1582                  zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 
     1583                  IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1584                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 
     1585                     ! 
     1586                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice 
     1587                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
     1588                     ! 
     1589                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     1590                     pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 
     1591                  ENDIF            
     1592                  ! 
     1593                  !                               ! -- check snow load -- ! 
     1594                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 
     1595                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 
     1596                  !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 
     1597                  zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
     1598                  IF( zvs_excess > 0._wp ) THEN 
     1599                     zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 
     1600                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice 
     1601                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
     1602                     ! 
     1603                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     1604                     pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess 
     1605                  ENDIF 
     1606                   
     1607               ENDIF 
     1608            END DO 
     1609         END DO 
     1610      END DO  
     1611      !                                           !-- correct pond fraction to avoid a_ip > a_i 
     1612      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:) 
     1613      ! 
     1614      ! 
     1615   END SUBROUTINE Hbig 
     1616    
    13601617#else 
    13611618   !!---------------------------------------------------------------------- 
  • NEMO/releases/release-4.0/src/ICE/icedyn_rhg.F90

    r10413 r10910  
    5858      !!-------------------------------------------------------------------- 
    5959      INTEGER, INTENT(in) ::   kt     ! ice time step 
     60      ! 
     61      INTEGER  ::   jl   ! dummy loop indices 
    6062      !!-------------------------------------------------------------------- 
    6163      ! controls 
     
    6870         WRITE(numout,*)'~~~~~~~~~~~' 
    6971      ENDIF 
    70  
    71       ! -------- 
    72       ! Rheology 
    73       ! --------    
     72      ! 
     73      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization 
     74         tau_icebfr(:,:) = 0._wp 
     75         DO jl = 1, jpl 
     76            WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
     77         END DO 
     78      ENDIF 
     79      ! 
     80      !--------------! 
     81      !== Rheology ==! 
     82      !--------------!    
    7483      SELECT CASE( nice_rhg ) 
    7584      !                                !------------------------! 
  • NEMO/releases/release-4.0/src/ICE/icewri.F90

    r10785 r10910  
    145145 
    146146        ! record presence of fast ice 
    147          WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk00(:,:) == 1._wp ) ; zfast(:,:) = 1._wp 
     147         WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk15(:,:) == 1._wp ) ; zfast(:,:) = 1._wp 
    148148         ELSEWHERE                                                ; zfast(:,:) = 0._wp 
    149149         END WHERE 
Note: See TracChangeset for help on using the changeset viewer.