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 10994 for NEMO/trunk/src/ICE/iceitd.F90 – NEMO

Ignore:
Timestamp:
2019-05-17T15:08:34+02:00 (5 years 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

File:
1 edited

Legend:

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

    r10069 r10994  
    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 
Note: See TracChangeset for help on using the changeset viewer.