Changeset 10993


Ignore:
Timestamp:
2019-05-17T15:07:59+02:00 (18 months ago)
Author:
clem
Message:

couple of commits to 1) deal with ice concentration reaching 1 (so, no more limitations in the max ice concentration imposed in the namelist), 2) deal with very small negative values that can occur due to roundoff errors

Location:
NEMO/releases/release-4.0/src/ICE
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • NEMO/releases/release-4.0/src/ICE/icecor.F90

    r10425 r10993  
    6666         WRITE(numout,*) '~~~~~~~' 
    6767      ENDIF 
    68       ! 
    6968      !                             !----------------------------------------------------- 
    70       !                             !  ice thickness must exceed himin (for ice diff)    ! 
     69      !                             !  ice thickness must exceed himin (for temp. diff.) ! 
    7170      !                             !----------------------------------------------------- 
    7271      WHERE( a_i(:,:,:) >= epsi20 )   ;   h_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:) 
     
    9796         END DO 
    9897      ENDIF 
    99  
    10098      !                             !----------------------------------------------------- 
    10199      !                             !  Rebin categories with thickness out of bounds     ! 
  • NEMO/releases/release-4.0/src/ICE/icectl.F90

    r10581 r10993  
    6969      !! 
    7070      REAL(wp) ::   zv, zs, zt, zfs, zfv, zft 
    71       REAL(wp) ::   zvmin, zamin, zamax  
     71      REAL(wp) ::   zvmin, zamin, zamax, zeimin, zesmin, zsmin 
    7272      REAL(wp) ::   zvtrp, zetrp 
    7373      REAL(wp) ::   zarea, zv_sill, zs_sill, zt_sill 
     
    141141         zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t  ) * zconv 
    142142 
    143          zvmin = glob_min( 'icectl', v_i ) 
    144          zamax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
    145          zamin = glob_min( 'icectl', a_i ) 
     143         zamax  = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
     144         zvmin  = glob_min( 'icectl', v_i ) 
     145         zamin  = glob_min( 'icectl', a_i ) 
     146         zsmin  = glob_min( 'icectl', sv_i ) 
     147         zeimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 
     148         zesmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) 
    146149 
    147150         ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  
     
    152155 
    153156         IF(lwp) THEN 
    154             IF ( ABS( zv   ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv 
    155             IF ( ABS( zs   ) > zs_sill )   WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs 
    156             IF ( ABS( zt   ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt 
    157             IF ( zvmin < -epsi10 )         WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin 
    158             IF ( zamax > MAX( rn_amax_n, rn_amax_s ) + epsi10   & 
    159                & .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' .AND. cd_routine /= 'Hbig' ) & 
    160                &                           WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
    161             IF ( zamin < -epsi10 )         WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
    162 !clem: the following check fails when using UMx advection scheme (see comments in icedyn_adv.F90) 
     157            ! check conservation issues 
     158            IF ( ABS( zv ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv 
     159            IF ( ABS( zs ) > zs_sill )   WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs 
     160            IF ( ABS( zt ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt 
     161            ! check maximum ice concentration 
     162            IF ( zamax > MAX( rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' )  & 
     163               &                         WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
     164            ! check negative values 
     165            IF ( zvmin  < 0. )           WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin 
     166            IF ( zamin  < 0. )           WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
     167            IF ( zsmin  < 0. )           WRITE(numout,*) 'violation s_i<0               (',cd_routine,') = ',zsmin 
     168            IF ( zeimin < 0. )           WRITE(numout,*) 'violation e_i<0               (',cd_routine,') = ',zeimin 
     169            IF ( zesmin < 0. )           WRITE(numout,*) 'violation e_s<0               (',cd_routine,') = ',zesmin 
     170!clem: the following check fails (I think...) 
    163171!            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'icedyn_adv' ) THEN 
    164172!                                           WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp 
  • NEMO/releases/release-4.0/src/ICE/icedyn.F90

    r10910 r10993  
    118118         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice 
    119119         CALL Hpiling                                                       ! -- simple pile-up (replaces ridging/rafting) 
     120         CALL ice_var_zapsmall                                              ! -- zap small areas 
    120121         ! 
    121122      CASE ( np_dynADV1D )         !==  pure advection ==!   (1D) 
     
    189190      ! controls 
    190191      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    191       ! 
    192       CALL ice_var_zapsmall                       !-- zap small areas 
    193192      ! 
    194193      at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
  • NEMO/releases/release-4.0/src/ICE/icedyn_rdgrft.F90

    r10944 r10993  
    142142      INTEGER, PARAMETER ::   jp_itermax = 20     
    143143      !!------------------------------------------------------------------- 
    144       ! clem: The redistribution of ice between categories can lead to small negative values (as for the remapping in ice_itd_rem) 
    145       !       likely due to truncation error ( i.e. 1. - 1. /= 0 ) 
    146       !       I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90) 
    147        
    148144      ! controls 
    149145      IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing 
     
    156152      ENDIF       
    157153 
    158       CALL ice_var_zapsmall   ! Zero out categories with very small areas 
    159  
    160154      !-------------------------------- 
    161155      ! 0) Identify grid cells with ice 
    162156      !-------------------------------- 
     157      at_i(:,:) = SUM( a_i, dim=3 ) 
     158      ! 
    163159      npti = 0   ;   nptidx(:) = 0 
    164160      ipti = 0   ;   iptidx(:) = 0 
    165161      DO jj = 1, jpj 
    166162         DO ji = 1, jpi 
    167             IF ( at_i(ji,jj) > 0._wp ) THEN 
     163            IF ( at_i(ji,jj) > epsi10 ) THEN 
    168164               npti           = npti + 1 
    169165               nptidx( npti ) = (jj - 1) * jpi + ji 
     
    178174         
    179175         ! just needed here 
    180          CALL tab_2d_1d( npti, nptidx(1:npti), zdivu(1:npti), divu_i(:,:) ) 
    181          CALL tab_2d_1d( npti, nptidx(1:npti), zdelt(1:npti), delta_i(:,:) ) 
     176         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i ) 
     177         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i ) 
    182178         ! needed here and in the iteration loop 
    183          CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i(:,:,:) ) 
    184          CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i(:,:,:) ) 
    185          CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i(:,:) ) 
     179         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i  ) 
     180         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  ) 
     181         CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i ) 
    186182 
    187183         DO ji = 1, npti 
     
    310306 
    311307      !                       ! Ice thickness needed for rafting 
    312       WHERE( pa_i(1:npti,:) > 0._wp )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 
    313       ELSEWHERE                         ;   zhi(1:npti,:) = 0._wp 
     308      WHERE( pa_i(1:npti,:) > epsi20 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 
     309      ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp 
    314310      END WHERE 
    315311 
     
    329325      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 
    330326      ! 
    331       WHERE( zasum(1:npti) > 0._wp )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
    332       ELSEWHERE                        ;   z1_asum(1:npti) = 0._wp 
     327      WHERE( zasum(1:npti) > epsi20 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
     328      ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp 
    333329      END WHERE 
    334330      ! 
     
    455451      ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.   
    456452      ! NOTE: 0 < aksum <= 1 
    457       WHERE( zaksum(1:npti) > 0._wp )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 
    458       ELSEWHERE                         ;   closing_gross(1:npti) = 0._wp 
     453      WHERE( zaksum(1:npti) > epsi20 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 
     454      ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp 
    459455      END WHERE 
    460456       
     
    466462         DO ji = 1, npti 
    467463            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 
    468             IF( zfac > pa_i(ji,jl) ) THEN 
     464            IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN 
    469465               closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice 
    470466            ENDIF 
     
    510506      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2 
    511507      REAL(wp), DIMENSION(jpij) ::   z1_ai            ! 1 / a 
     508      REAL(wp), DIMENSION(jpij) ::   zvti             ! sum(v_i) 
    512509      ! 
    513510      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrft     ! snow energy of rafting ice 
     
    518515      INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation 
    519516      !!------------------------------------------------------------------- 
    520  
     517      ! 
     518      zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 )   ! total ice volume 
     519      ! 
    521520      ! 1) Change in open water area due to closing and opening 
    522521      !-------------------------------------------------------- 
     
    535534            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging 
    536535 
    537                z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
    538  
     536               IF( a_i_2d(ji,jl1) > epsi20 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
     537               ELSE                                 ;   z1_ai(ji) = 0._wp 
     538               ENDIF 
     539                
    539540               ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2) 
    540541               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice 
     
    549550 
    550551               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges 
    551                vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg 
     552               IF    ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg                                           ! v <= 10m then porosity = rn_porordg 
     553               ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp                                                                         ! v >= 20m then porosity = 0 
     554               ELSE                           ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg * MAX( 0._wp, 2._wp - 0.1_wp * zvti(ji) ) ! v > 10m and v < 20m then porosity = linear transition to 0 
     555               ENDIF 
    552556               ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?) 
    553557 
    554558               ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei) 
    555559               virdg1     = v_i_2d (ji,jl1)   * afrdg 
    556                virdg2(ji) = v_i_2d (ji,jl1)   * afrdg * ( 1. + rn_porordg ) 
     560               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg + vsw 
    557561               vsrdg(ji)  = v_s_2d (ji,jl1)   * afrdg 
    558562               sirdg1     = sv_i_2d(ji,jl1)   * afrdg 
     
    726730      END DO ! jl1 
    727731      ! 
     732      ! roundoff errors 
     733      !---------------- 
    728734      ! In case ridging/rafting lead to very small negative values (sometimes it happens) 
    729       WHERE( a_i_2d(1:npti,:) < 0._wp )   a_i_2d(1:npti,:) = 0._wp 
    730       WHERE( v_i_2d(1:npti,:) < 0._wp )   v_i_2d(1:npti,:) = 0._wp 
     735      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 
    731736      ! 
    732737   END SUBROUTINE rdgrft_shift 
     
    854859         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
    855860         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) ) 
    856  
     861         ! 
    857862         !                 !---------------------! 
    858863      CASE( 2 )            !==  from 1D to 2D  ==! 
  • NEMO/releases/release-4.0/src/ICE/iceitd.F90

    r10069 r10993  
    2121   USE ice1D          ! sea-ice: thermodynamic variables 
    2222   USE ice            ! sea-ice: variables 
     23   USE icevar         ! sea-ice: operations 
    2324   USE icectl         ! sea-ice: conservation tests 
    2425   USE icetab         ! sea-ice: convert 1D<=>2D 
     
    9192      !  1) Identify grid cells with ice 
    9293      !----------------------------------------------------------------------------------------------- 
     94      at_i(:,:) = SUM( a_i, dim=3 ) 
     95      ! 
    9396      npti = 0   ;   nptidx(:) = 0 
    9497      DO jj = 1, jpj 
     
    249252            ! --- g(h) for each thickness category --- !   
    250253            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 
     254               &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL    (1:npti,jl), hR    (1:npti,jl)   )   ! out 
    252255            ! 
    253256         END DO 
     
    389392      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdvice   ! ice volume transferred across boundary 
    390393      ! 
    391       INTEGER  ::   ji, jj, jl, jk     ! dummy loop indices 
    392       INTEGER  ::   ii, ij, jl2, jl1   ! local integers 
     394      INTEGER  ::   ji, jl, jk         ! dummy loop indices 
     395      INTEGER  ::   jl2, jl1           ! local integers 
    393396      REAL(wp) ::   ztrans             ! ice/snow transferred 
    394       REAL(wp), DIMENSION(jpij)     ::   zworka, zworkv   ! workspace 
    395       REAL(wp), DIMENSION(jpij,jpl) ::   zaTsfn           !  -    - 
     397      REAL(wp), DIMENSION(jpij)            ::   zworka, zworkv   ! workspace 
     398      REAL(wp), DIMENSION(jpij,jpl)        ::   zaTsfn           !  -    - 
     399      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_2d 
     400      REAL(wp), DIMENSION(jpij,nlay_s,jpl) ::   ze_s_2d 
    396401      !!------------------------------------------------------------------ 
    397402          
     
    405410      CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 
    406411      CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 
     412      DO jl = 1, jpl 
     413         DO jk = 1, nlay_s 
     414            CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 
     415         END DO 
     416         DO jk = 1, nlay_i 
     417            CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 
     418         END DO 
     419      END DO 
     420      ! to correct roundoff errors on a_i 
     421      CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti), rn_amax_2d ) 
    407422 
    408423      !---------------------------------------------------------------------------------------------- 
     
    435450               ELSE                                  ;   zworka(ji) = 0._wp 
    436451               ENDIF 
    437                ! 
    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) 
    441452               ! 
    442453               a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas 
     
    476487         ! 
    477488         DO jk = 1, nlay_s         !--- Snow heat content 
    478             ! 
    479489            DO ji = 1, npti 
    480                ii = MOD( nptidx(ji) - 1, jpi ) + 1 
    481                ij = ( nptidx(ji) - 1 ) / jpi + 1 
    482490               ! 
    483491               jl1 = kdonor(ji,jl) 
     
    487495                  ELSE                ;  jl2 = jl 
    488496                  ENDIF 
    489                   ! 
    490                   ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji) 
    491                   e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans 
    492                   e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans 
     497                  ztrans             = ze_s_2d(ji,jk,jl1) * zworkv(ji) 
     498                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) - ztrans 
     499                  ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ztrans 
    493500               ENDIF 
    494501            END DO 
     
    497504         DO jk = 1, nlay_i         !--- Ice heat content 
    498505            DO ji = 1, npti 
    499                ii = MOD( nptidx(ji) - 1, jpi ) + 1 
    500                ij = ( nptidx(ji) - 1 ) / jpi + 1 
    501506               ! 
    502507               jl1 = kdonor(ji,jl) 
     
    506511                  ELSE                ;  jl2 = jl 
    507512                  ENDIF 
    508                   ! 
    509                   ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji) 
    510                   e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans 
    511                   e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans 
     513                  ztrans             = ze_i_2d(ji,jk,jl1) * zworkv(ji) 
     514                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) - ztrans 
     515                  ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + ztrans 
    512516               ENDIF 
    513517            END DO 
     
    515519         ! 
    516520      END DO                   ! boundaries, 1 to jpl-1 
     521 
     522      !------------------- 
     523      ! 3) roundoff errors 
     524      !------------------- 
     525      ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 
     526      !       because of truncation error ( i.e. 1. - 1. /= 0 ) 
     527      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 
     528 
     529      ! at_i must be <= rn_amax 
     530      zworka(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 ) 
     531      DO jl  = 1, jpl 
     532         WHERE( zworka(1:npti) > rn_amax_1d(1:npti) )   & 
     533            &   a_i_2d(1:npti,jl) = a_i_2d(1:npti,jl) * rn_amax_1d(1:npti) / zworka(1:npti) 
     534      END DO 
    517535       
    518536      !------------------------------------------------------------------------------- 
    519       ! 3) Update ice thickness and temperature 
     537      ! 4) Update ice thickness and temperature 
    520538      !------------------------------------------------------------------------------- 
    521539      WHERE( a_i_2d(1:npti,:) >= epsi20 ) 
     
    536554      CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 
    537555      CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 
     556      DO jl = 1, jpl 
     557         DO jk = 1, nlay_s 
     558            CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 
     559         END DO 
     560         DO jk = 1, nlay_i 
     561            CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 
     562         END DO 
     563      END DO 
    538564      ! 
    539565   END SUBROUTINE itd_shiftice 
     
    558584      ! 
    559585      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution'  
     586      ! 
     587      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    560588      ! 
    561589      jdonor(:,:) = 0 
     
    635663      END DO 
    636664      ! 
     665      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     666      ! 
    637667   END SUBROUTINE ice_itd_reb 
    638668 
  • NEMO/releases/release-4.0/src/ICE/icestp.F90

    r10931 r10993  
    323323         WRITE(numout,*) '         maximum ice concentration for SH                              = ', rn_amax_s 
    324324      ENDIF 
     325      !                                        !--- change max ice concentration for roundoff errors 
     326      rn_amax_n = MIN( rn_amax_n, 1._wp - epsi10 ) 
     327      rn_amax_s = MIN( rn_amax_s, 1._wp - epsi10 ) 
    325328      !                                        !--- check consistency 
    326329      IF ( jpl > 1 .AND. ln_virtual_itd ) THEN 
  • NEMO/releases/release-4.0/src/ICE/icethd.F90

    r10931 r10993  
    102102      ENDIF 
    103103       
    104       CALL ice_var_glo2eqv 
    105  
    106104      !---------------------------------------------! 
    107105      ! computation of friction velocity at T points 
     
    162160            qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr ) 
    163161 
    164             ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting  
    165             IF( zqld > 0._wp ) THEN 
     162            ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting  
     163            ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget 
     164            IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN 
    166165               fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 
    167166               qlead(ji,jj) = 0._wp 
     
    242241         ! 
    243242      END DO 
    244  
     243      ! 
    245244      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    246       ! 
    247                               CALL ice_var_zapsmall                 ! --- remove very small ice concentration (<1e-10) --- ! 
    248       !                                                             !     & make sure at_i=sum(a_i) & ato_i=1 where at_i=0 
    249245      !                    
    250246      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- ! 
  • NEMO/releases/release-4.0/src/ICE/icethd_do.F90

    r10425 r10993  
    114114      IF( ln_icediachk )   CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft ) 
    115115 
    116       CALL ice_var_agg(1) 
    117       CALL ice_var_glo2eqv 
    118  
     116      at_i(:,:) = SUM( a_i, dim=3 ) 
    119117      !------------------------------------------------------------------------------! 
    120118      ! 1) Collection thickness of ice formed in leads and polynyas 
     
    319317 
    320318         ! --- lateral ice growth --- ! 
    321          ! If lateral ice growth gives an ice concentration gt 1, then 
     319         ! If lateral ice growth gives an ice concentration > amax, then 
    322320         ! we keep the excessive volume in memory and attribute it later to bottom accretion 
    323321         DO ji = 1, npti 
    324             IF ( za_newice(ji) >  ( rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN 
    325                zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - at_i_1d(ji) ) 
     322            IF ( za_newice(ji) >  MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN ! max is for roundoff error 
     323               zda_res(ji)   = za_newice(ji) - MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) 
    326324               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji)  
    327                za_newice(ji) = za_newice(ji) - zda_res  (ji) 
    328                zv_newice(ji) = zv_newice(ji) - zdv_res  (ji) 
     325               za_newice(ji) = MAX( 0._wp, za_newice(ji) - zda_res  (ji) ) 
     326               zv_newice(ji) = MAX( 0._wp, zv_newice(ji) - zdv_res  (ji) ) 
    329327            ELSE 
    330328               zda_res(ji) = 0._wp 
  • NEMO/releases/release-4.0/src/ICE/icevar.F90

    r10944 r10993  
    4444   !!   ice_var_salprof1d : salinity profile in the ice 1D 
    4545   !!   ice_var_zapsmall  : remove very small area and volume 
    46    !!   ice_var_zapneg    : remove negative ice fields (to debug the advection scheme UM3-5) 
     46   !!   ice_var_zapneg    : remove negative ice fields 
     47   !!   ice_var_roundoff  : remove negative values arising from roundoff erros 
    4748   !!   ice_var_itd       : convert 1-cat to jpl-cat 
    4849   !!   ice_var_itd2      : convert N-cat to jpl-cat 
     
    7172   PUBLIC   ice_var_zapsmall 
    7273   PUBLIC   ice_var_zapneg 
     74   PUBLIC   ice_var_roundoff 
    7375   PUBLIC   ice_var_itd 
    7476   PUBLIC   ice_var_itd2 
     
    622624   END SUBROUTINE ice_var_zapneg 
    623625 
     626 
     627   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     628      !!------------------------------------------------------------------- 
     629      !!                   ***  ROUTINE ice_var_roundoff *** 
     630      !! 
     631      !! ** Purpose :   Remove negative sea ice values arising from roundoff errors 
     632      !!------------------------------------------------------------------- 
     633      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_i       ! ice concentration 
     634      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_i       ! ice volume 
     635      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_s       ! snw volume 
     636      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   psv_i      ! salt content 
     637      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   poa_i      ! age content 
     638      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
     639      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     640      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
     641      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     642      !!------------------------------------------------------------------- 
     643      ! 
     644      WHERE( pa_i (1:npti,:)   < 0._wp .AND. pa_i (1:npti,:)   > -epsi10 )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0 
     645      WHERE( pv_i (1:npti,:)   < 0._wp .AND. pv_i (1:npti,:)   > -epsi10 )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0 
     646      WHERE( pv_s (1:npti,:)   < 0._wp .AND. pv_s (1:npti,:)   > -epsi10 )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0 
     647      WHERE( psv_i(1:npti,:)   < 0._wp .AND. psv_i(1:npti,:)   > -epsi10 )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0 
     648      WHERE( poa_i(1:npti,:)   < 0._wp .AND. poa_i(1:npti,:)   > -epsi10 )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0 
     649      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0 
     650      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0 
     651      IF ( ln_pnd_H12 ) THEN 
     652         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
     653         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 
     654      ENDIF 
     655      ! 
     656   END SUBROUTINE ice_var_roundoff 
     657    
    624658    
    625659   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i ) 
     
    656690      INTEGER  :: idim, i_fill, jl0   
    657691      REAL(wp) :: zarg, zV, zconv, zdh, zdv 
    658       REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables 
    659       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i ! output ice/snow variables 
     692      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input  ice/snow variables 
     693      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables 
    660694      INTEGER , DIMENSION(4)                  ::   itest 
    661695      !!------------------------------------------------------------------- 
     
    786820      !! 
    787821      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1 
    788       !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
     822       !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
    789823      !!               
    790824      !!               3) Expand the filling to the empty cat between jlmin and jlmax  
Note: See TracChangeset for help on using the changeset viewer.