Ignore:
Timestamp:
2019-10-18T17:47:20+02:00 (13 months ago)
Author:
clem
Message:

last commit: avoid division by 0 occuring on some rare occasions (CMMC and ECMWF issues). It might be machine dependant.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/releases/release-4.0.1/src/ICE/icedyn_rdgrft.F90

    r11560 r11730  
    8686      !!                ***  ROUTINE ice_dyn_rdgrft_alloc *** 
    8787      !!------------------------------------------------------------------- 
    88       ALLOCATE( closing_net(jpij), opning(jpij)   , closing_gross(jpij),   & 
    89          &      apartf(jpij,0:jpl), hrmin(jpij,jpl), hraft(jpij,jpl)    , aridge(jpij,jpl),  & 
    90          &      hrmax(jpij,jpl), hi_hrdg(jpij,jpl)  , araft (jpij,jpl),  & 
     88      ALLOCATE( closing_net(jpij)  , opning(jpij)      , closing_gross(jpij) ,               & 
     89         &      apartf(jpij,0:jpl) , hrmin  (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), & 
     90         &      hrmax (jpij,jpl)   , hi_hrdg(jpij,jpl) , araft(jpij,jpl) ,                   & 
    9191         &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 
    9292 
     
    137137      REAL(wp) ::   zfac                       ! local scalar 
    138138      INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not 
    139       REAL(wp), DIMENSION(jpij) ::   zdivu_adv     ! divu as implied by transport scheme  (1/s) 
    140139      REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i 
    141140      ! 
     
    175174         
    176175         ! just needed here 
    177          CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i  ) 
    178176         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i ) 
    179177         ! needed here and in the iteration loop 
     178         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i) ! zdivu is used as a work array here (no change in divu_i) 
    180179         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i   ) 
    181180         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i   ) 
     
    187186            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 
    188187            ! 
    189             ! divergence given by the advection scheme 
    190             !   (which may not be equal to divu as computed from the velocity field) 
    191             IF    ( ln_adv_Pra ) THEN 
    192                zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 
    193             ELSEIF( ln_adv_UMx ) THEN 
    194                zdivu_adv(ji) = zdivu(ji) 
    195             ENDIF 
    196             ! 
    197             IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )   ! make sure the closing rate is large enough 
    198             !                                                                                        ! to give asum = 1.0 after ridging 
     188            IF( zdivu(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) )   ! make sure the closing rate is large enough 
     189            !                                                                                ! to give asum = 1.0 after ridging 
    199190            ! Opening rate (non-negative) that will give asum = 1.0 after ridging. 
    200             opning(ji) = closing_net(ji) + zdivu_adv(ji) 
     191            opning(ji) = closing_net(ji) + zdivu(ji) 
    201192         END DO 
    202193         ! 
     
    215206               ato_i_1d   (ipti)   = ato_i_1d   (ji) 
    216207               closing_net(ipti)   = closing_net(ji) 
    217                zdivu_adv  (ipti)   = zdivu_adv  (ji) 
     208               zdivu      (ipti)   = zdivu      (ji) 
    218209               opning     (ipti)   = opning     (ji) 
    219210            ENDIF 
     
    259250               ELSE 
    260251                  iterate_ridging  = 1 
    261                   zdivu_adv  (ji) = zfac * r1_rdtice 
    262                   closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) ) 
    263                   opning     (ji) = MAX( 0._wp,  zdivu_adv(ji) ) 
     252                  zdivu      (ji) = zfac * r1_rdtice 
     253                  closing_net(ji) = MAX( 0._wp, -zdivu(ji) ) 
     254                  opning     (ji) = MAX( 0._wp,  zdivu(ji) ) 
    264255               ENDIF 
    265256            END DO 
     
    309300 
    310301      !                       ! Ice thickness needed for rafting 
    311       WHERE( pa_i(1:npti,:) > epsi20 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 
     302      WHERE( pa_i(1:npti,:) > epsi10 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 
    312303      ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp 
    313304      END WHERE 
     
    328319      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 
    329320      ! 
    330       WHERE( zasum(1:npti) > epsi20 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
     321      WHERE( zasum(1:npti) > epsi10 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
    331322      ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp 
    332323      END WHERE 
     
    454445      ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.   
    455446      ! NOTE: 0 < aksum <= 1 
    456       WHERE( zaksum(1:npti) > epsi20 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 
     447      WHERE( zaksum(1:npti) > epsi10 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 
    457448      ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp 
    458449      END WHERE 
     
    537528            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging 
    538529 
    539                IF( a_i_2d(ji,jl1) > epsi20 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
     530               IF( a_i_2d(ji,jl1) > epsi10 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
    540531               ELSE                                 ;   z1_ai(ji) = 0._wp 
    541532               ENDIF 
     
    595586               ! virtual salt flux to keep salinity constant 
    596587               IF( nn_icesal /= 2 )  THEN 
    597                   sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i 
     588                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )       ! ridge salinity = s_i 
    598589                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_rdtice  &  ! put back sss_m into the ocean 
    599590                     &                            - s_i_1d(ji) * vsw * rhoi * r1_rdtice     ! and get  s_i  from the ocean  
Note: See TracChangeset for help on using the changeset viewer.