Changeset 9880


Ignore:
Timestamp:
2018-07-05T17:43:44+02:00 (2 years ago)
Author:
clem
Message:

sea-ice: remove negative values mainly created by the advection schemes (UMx and Prather) + adapt BDY to ensure that all the fields are coherent with the rest of the code (max concentration etc)

Location:
NEMO/trunk/src
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icedyn_adv.F90

    r9604 r9880  
    9595      END SELECT 
    9696 
     97      !---------------------------- 
     98      ! Debug the advection schemes 
     99      !---------------------------- 
     100      ! clem: The 2 advection schemes above are not strictly positive. 
     101      !       In Prather, advected fields are bounded by 0 in the routine with a MAX(0,field) ==> likely conservation issues 
     102      !       In UMx    , advected fields are not bounded and negative values can appear. 
     103      !                   These values are usually very small but in some occasions they can also be non-negligible 
     104      !                   Therefore one needs to bound the advected fields by 0 (though this is not a clean fix) 
     105      ! ==> 1) remove negative ice areas and volumes (conservation is ensure) 
     106      CALL ice_var_zapsmall  
     107      ! ==> 2) remove remaining negative advected fields (conservation is not preserved) 
     108      WHERE( v_s (:,:,:)   < 0._wp )   v_s (:,:,:)   = 0._wp 
     109      WHERE( sv_i(:,:,:)   < 0._wp )   sv_i(:,:,:)   = 0._wp 
     110      WHERE( e_i (:,:,:,:) < 0._wp )   e_i (:,:,:,:) = 0._wp 
     111      WHERE( e_s (:,:,:,:) < 0._wp )   e_s (:,:,:,:) = 0._wp 
     112 
    97113      !------------ 
    98114      ! diagnostics 
  • NEMO/trunk/src/ICE/icedyn_rdgrft.F90

    r9604 r9880  
    143143      INTEGER, PARAMETER ::   jp_itermax = 20     
    144144      !!------------------------------------------------------------------- 
     145      ! clem: The redistribution of ice between categories can lead to small negative values (as for the remapping in ice_itd_rem) 
     146      !       likely due to truncation error ( i.e. 1. - 1. /= 0 ) 
     147      !       I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90) 
     148       
    145149      ! controls 
    146150      IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing 
  • NEMO/trunk/src/ICE/iceitd.F90

    r9604 r9880  
    9595         DO ji = 1, jpi 
    9696            IF ( at_i(ji,jj) > epsi10 ) THEN 
    97                npti         = npti + 1 
     97               npti = npti + 1 
    9898               nptidx( npti ) = (jj - 1) * jpi + ji 
    9999            ENDIF 
     
    111111         CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i   ) 
    112112         CALL tab_3d_2d( npti, nptidx(1:npti), h_ib_2d(1:npti,1:jpl), h_i_b ) 
    113          CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i    ) 
    114          CALL tab_3d_2d( npti, nptidx(1:npti), a_ib_2d (1:npti,1:jpl), a_i_b ) 
     113         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i   ) 
     114         CALL tab_3d_2d( npti, nptidx(1:npti), a_ib_2d(1:npti,1:jpl), a_i_b ) 
    115115         ! 
    116116         DO jl = 1, jpl 
    117117            ! Compute thickness change in each ice category 
    118118            DO ji = 1, npti 
    119                IF( a_i_2d(ji,jl) > epsi10 )   zdhice(ji,jl) = h_i_2d(ji,jl) - h_ib_2d(ji,jl)   ! clem: if statement is necessary for repro 
     119               IF( a_i_2d(ji,jl) > epsi10 )   zdhice(ji,jl) = h_i_2d(ji,jl) - h_ib_2d(ji,jl) 
    120120            END DO 
    121121         END DO 
     
    130130               ! 
    131131               IF    ( a_ib_2d(ji,jl) >  epsi10 .AND. a_ib_2d(ji,jl+1) >  epsi10 ) THEN   ! a(jl+1) & a(jl) /= 0 
    132                   zslope           = ( zdhice(ji,jl+1) - zdhice(ji,jl) ) / ( h_ib_2d(ji,jl+1) - h_ib_2d(ji,jl) ) 
     132                  zslope        = ( zdhice(ji,jl+1) - zdhice(ji,jl) ) / ( h_ib_2d(ji,jl+1) - h_ib_2d(ji,jl) ) 
    133133                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) + zslope * ( hi_max(jl) - h_ib_2d(ji,jl) ) 
    134134               ELSEIF( a_ib_2d(ji,jl) >  epsi10 .AND. a_ib_2d(ji,jl+1) <= epsi10 ) THEN   ! a(jl+1)=0 => Hn* = Hn + fn*dt 
     
    198198            ! 
    199199            CALL tab_2d_1d( npti, nptidx(1:npti), h_ib_1d(1:npti), h_i_b(:,:,jl) ) 
    200             CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl)  ) 
    201             CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl)    ) 
    202             CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl)    ) 
     200            CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i  (:,:,jl) ) 
     201            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i  (:,:,jl) ) 
     202            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i  (:,:,jl) ) 
    203203            ! 
    204204            IF( jl == 1 ) THEN 
     
    206206               ! --- g(h) for category 1 --- ! 
    207207               CALL itd_glinear( zhb0(1:npti)  , zhb1(1:npti)  , h_ib_1d(1:npti)  , a_i_1d(1:npti)  ,  &   ! in 
    208                   &              g0  (1:npti,1), g1  (1:npti,1), hL      (1:npti,1), hR    (1:npti,1)   )   ! out 
     208                  &              g0  (1:npti,1), g1  (1:npti,1), hL     (1:npti,1), hR    (1:npti,1)   )   ! out 
    209209                  ! 
    210210               ! Area lost due to melting of thin ice 
     
    228228                           ! Remove area, conserving volume 
    229229                           h_i_1d(ji) = h_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 ) 
    230                            a_i_1d(ji)  = a_i_1d(ji) - zda0 
    231                            v_i_1d(ji)  = a_i_1d(ji) * h_i_1d(ji) ! useless ? 
     230                           a_i_1d(ji) = a_i_1d(ji) - zda0 
     231                           v_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) ! useless ? 
    232232                        ENDIF 
    233233                        ! 
     
    241241               END DO 
    242242               ! 
    243                CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl)  ) 
    244                CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl)    ) 
    245                CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl)    ) 
     243               CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 
     244               CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 
     245               CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) ) 
    246246               ! 
    247247            ENDIF ! jl=1 
     
    249249            ! --- g(h) for each thickness category --- !   
    250250            CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti)   , a_i_1d(1:npti)   ,  &   ! in 
    251                &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL     (1:npti,jl), hR    (1:npti,jl)   )   ! out 
     251               &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL     (1:npti,jl), hR   (1:npti,jl)   )   ! out 
    252252            ! 
    253253         END DO 
     
    300300         DO ji = 1, npti 
    301301            IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN 
    302                a_i_1d (ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin  
     302               a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin  
    303303               IF( ln_pnd_H12 )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 
    304304               h_i_1d(ji) = rn_himin 
     
    436436               ENDIF 
    437437               ! 
     438               ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 
     439               !       because of truncation error ( i.e. 1. - 1. /= 0 ) 
     440               !       I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90) 
     441               ! 
    438442               a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas 
    439443               a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl) 
     
    445449               v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans 
    446450               v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans  
    447                !           
    448                !                                                     ! Ice age  
    449                ztrans          = oa_i_2d(ji,jl1) * zworka(ji) 
     451               ! 
     452               ztrans          = oa_i_2d(ji,jl1) * zworka(ji)        ! Ice age 
    450453               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - ztrans 
    451454               oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ztrans 
     
    455458               sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ztrans 
    456459               ! 
    457                !                                                     ! Surface temperature 
    458                ztrans          = zaTsfn(ji,jl1) * zworka(ji) 
     460               ztrans          = zaTsfn(ji,jl1) * zworka(ji)         ! Surface temperature 
    459461               zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans 
    460462               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans 
    461463               !   
    462464               IF ( ln_pnd_H12 ) THEN 
    463                   !                                                  ! Pond fraction 
    464                   ztrans          = a_ip_2d(ji,jl1) * zworka(ji) 
     465                  ztrans          = a_ip_2d(ji,jl1) * zworka(ji)     ! Pond fraction 
    465466                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans 
    466467                  a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans 
    467                   !                                                  ! Pond volume (also proportional to da/a) 
    468                   ztrans          = v_ip_2d(ji,jl1) * zworka(ji) 
     468                  !                                               
     469                  ztrans          = v_ip_2d(ji,jl1) * zworka(ji)     ! Pond volume (also proportional to da/a) 
    469470                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 
    470471                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 
     
    519520      !------------------------------------------------------------------------------- 
    520521      WHERE( a_i_2d(1:npti,:) >= epsi20 ) 
    521          h_i_2d(1:npti,:)  =  v_i_2d(1:npti,:) / a_i_2d(1:npti,:)  
     522         h_i_2d (1:npti,:)  =  v_i_2d(1:npti,:) / a_i_2d(1:npti,:)  
    522523         t_su_2d(1:npti,:)  =  zaTsfn(1:npti,:) / a_i_2d(1:npti,:)  
    523524      ELSEWHERE 
    524          h_i_2d(1:npti,:)  = 0._wp 
     525         h_i_2d (1:npti,:)  = 0._wp 
    525526         t_su_2d(1:npti,:)  = rt0 
    526527      END WHERE 
     
    568569         DO jj = 1, jpj 
    569570            DO ji = 1, jpi 
    570                IF( a_i(ji,jj,jl) > epsi10 .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN 
     571               IF( a_i(ji,jj,jl) > 0._wp .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN 
    571572                  npti = npti + 1 
    572573                  nptidx( npti ) = (jj - 1) * jpi + ji                   
     
    575576         END DO 
    576577         ! 
    577 !!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl)  ) 
    578          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl)    ) 
    579          CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl)    ) 
     578!!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 
     579         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 
     580         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) ) 
    580581         ! 
    581582         DO ji = 1, npti 
    582583            jdonor(ji,jl)  = jl  
    583584            ! how much of a_i you send in cat sup is somewhat arbitrary 
    584 !!clem: these do not work properly after a restart (I do not know why) 
     585!!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    585586!!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji)   
    586587!!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 
    587 !!clem: these do not work properly after a restart (I do not know why) 
    588 !!            zdaice(ji,jl)  = a_i_1d(ji) 
    589 !!            zdvice(ji,jl)  = v_i_1d(ji) 
     588!!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
     589!!          zdaice(ji,jl)  = a_i_1d(ji) 
     590!!          zdvice(ji,jl)  = v_i_1d(ji) 
    590591!!clem: these are from UCL and work ok 
    591592            zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp 
     
    609610         DO jj = 1, jpj 
    610611            DO ji = 1, jpi 
    611                IF( a_i(ji,jj,jl+1) > epsi10 .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN 
     612               IF( a_i(ji,jj,jl+1) > 0._wp .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN 
    612613                  npti = npti + 1 
    613614                  nptidx( npti ) = (jj - 1) * jpi + ji                   
     
    616617         END DO 
    617618         ! 
    618          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl+1)    ) ! jl+1 is ok 
    619          CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl+1)    ) ! jl+1 is ok 
     619         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok 
     620         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok 
    620621         DO ji = 1, npti 
    621622            jdonor(ji,jl) = jl + 1 
  • NEMO/trunk/src/OCE/BDY/bdyice.F90

    r9657 r9880  
    1818   USE ice             ! sea-ice: variables 
    1919   USE icevar          ! sea-ice: operations 
    20    USE iceitd          ! sea-ice: rebining 
     20   USE icecor          ! sea-ice: corrections 
    2121   USE icectl          ! sea-ice: control prints 
    2222   USE phycst          ! physical constant 
     
    7373      END DO 
    7474      ! 
    75       CALL ice_var_zapsmall 
     75      CALL ice_cor( kt , 0 )      ! -- In case categories are out of bounds, do a remapping 
     76      !                           !    i.e. inputs have not the same ice thickness distribution  
     77      !                           !    (set by rn_himean) than the regional simulation 
    7678      CALL ice_var_agg(1) 
    7779      ! 
     
    250252      END DO !jl 
    251253      ! 
    252       ! --- In case categories are out of bounds, do a remapping --- ! 
    253       !     i.e. inputs have not the same ice thickness distribution  
    254       !          (set by rn_himean) than the regional simulation 
    255       IF( jpl > 1 )   CALL ice_itd_reb( kt ) 
     254!!      ! --- In case categories are out of bounds, do a remapping --- ! 
     255!!      !     i.e. inputs have not the same ice thickness distribution  
     256!!      !          (set by rn_himean) than the regional simulation 
     257!!      IF( jpl > 1 )   CALL ice_itd_reb( kt ) 
    256258      !       
    257259   END SUBROUTINE bdy_ice_frs 
Note: See TracChangeset for help on using the changeset viewer.