Ignore:
Timestamp:
2017-09-21T17:25:29+02:00 (3 years ago)
Author:
clem
Message:

put ridge/raft in 1D

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90

    r8534 r8554  
    2121   USE ice1D          ! sea-ice: thermodynamics 
    2222   USE ice            ! sea-ice: variables 
     23   USE icetab         ! sea-ice: 1D <==> 2D transformation 
    2324   USE icevar         ! sea-ice: operations 
    2425   USE icectl         ! sea-ice: control prints 
     
    3940 
    4041   ! Variables shared among ridging subroutines 
    41    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area 
    42    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged 
    43    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apartf   ! participation function; fraction of ridging/closing associated w/ category n 
    44    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness 
    45    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness 
    46    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice 
    47    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hi_hrdg  ! thickness of ridging ice / mean ridge thickness 
    48    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging 
    49    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting 
     42   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   closing_net     ! net rate at which area is removed    (1/s) 
     43      !                                                 ! (ridging ice area - area of new ridges) / dt 
     44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   opning         ! rate of opening due to divergence/shear 
     45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   closing_gross  ! rate at which area removed, not counting area of new ridges 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   apartf   ! participation function; fraction of ridging/closing associated w/ category n 
     47   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hrmin    ! minimum ridge thickness 
     48   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hrmax    ! maximum ridge thickness 
     49   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hraft    ! thickness of rafted ice 
     50   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hi_hrdg  ! thickness of ridging ice / mean ridge thickness 
     51   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   aridge   ! participating ice ridging 
     52   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   araft    ! participating ice rafting 
     53   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_i_2d 
     54   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_s_2d 
    5055   ! 
    5156   REAL(wp), PARAMETER ::   hrdg_hi_min = 1.1_wp    ! min ridge thickness multiplier: min(hrdg/hi) 
    5257   REAL(wp), PARAMETER ::   hi_hrft     = 0.5_wp    ! rafting multipliyer: (hi/hraft) 
    53    REAL(wp)            ::   zdrho                   !  
    5458   ! 
    5559   ! ** namelist (namdyn_rdgrft) ** 
     
    8387 
    8488   INTEGER FUNCTION ice_dyn_rdgrft_alloc() 
    85       !!---------------------------------------------------------------------! 
     89      !!------------------------------------------------------------------- 
    8690      !!                ***  ROUTINE ice_dyn_rdgrft_alloc *** 
    87       !!---------------------------------------------------------------------! 
    88       ALLOCATE( asum (jpi,jpj)     , apartf (jpi,jpj,0:jpl) , aksum (jpi,jpj)     ,     & 
    89          &      hrmin(jpi,jpj,jpl) , hraft  (jpi,jpj,jpl)   , aridge(jpi,jpj,jpl) ,     & 
    90          &      hrmax(jpi,jpj,jpl) , hi_hrdg(jpi,jpj,jpl)   , araft (jpi,jpj,jpl) , STAT=ice_dyn_rdgrft_alloc ) 
     91      !!------------------------------------------------------------------- 
     92      ALLOCATE( closing_net(jpij), opning(jpij)   , closing_gross(jpij),   & 
     93         &      apartf(jpij,0:jpl), hrmin(jpij,jpl), hraft(jpij,jpl)    , aridge(jpij,jpl),   & 
     94         &      hrmax(jpij,jpl), hi_hrdg(jpij,jpl)  , araft (jpij,jpl),  & 
     95         &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 
    9196 
    9297      IF( lk_mpp                    )   CALL mpp_sum ( ice_dyn_rdgrft_alloc ) 
     
    97102 
    98103   SUBROUTINE ice_dyn_rdgrft( kt ) 
    99       !!---------------------------------------------------------------------! 
     104      !!------------------------------------------------------------------- 
    100105      !!                ***  ROUTINE ice_dyn_rdgrft *** 
    101106      !! 
     
    121126      !!     This routine is based on CICE code and authors William H. Lipscomb, 
    122127      !!  and Elizabeth C. Hunke, LANL are gratefully acknowledged 
    123       !!--------------------------------------------------------------------! 
     128      !!------------------------------------------------------------------- 
    124129      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    125130      !! 
    126       INTEGER  ::   ji, jj, jk, jl     ! dummy loop index 
    127       INTEGER  ::   niter              ! local integer  
    128       INTEGER  ::   iterate_ridging    ! if =1, repeat the ridging 
    129       REAL(wp) ::   zfac               ! local scalar 
    130       REAL(wp), DIMENSION(jpi,jpj) ::   closing_net     ! net rate at which area is removed    (1/s) 
    131       !                                                 ! (ridging ice area - area of new ridges) / dt 
    132       REAL(wp), DIMENSION(jpi,jpj) ::   divu_adv        ! divu as implied by transport scheme  (1/s) 
    133       REAL(wp), DIMENSION(jpi,jpj) ::   opning          ! rate of opening due to divergence/shear 
    134       REAL(wp), DIMENSION(jpi,jpj) ::   closing_gross   ! rate at which area removed, not counting area of new ridges 
     131      INTEGER  ::   ji, jj, jk, jl             ! dummy loop index 
     132      INTEGER  ::   niter, iterate_ridging     ! local integer  
     133      INTEGER  ::   nidx2                      ! local integer 
     134      REAL(wp) ::   zfac                       ! local scalar 
     135      INTEGER , DIMENSION(jpij) ::   idxice2       ! compute ridge/raft or not 
     136      REAL(wp), DIMENSION(jpij) ::   zdivu_adv     ! divu as implied by transport scheme  (1/s) 
     137      REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i 
    135138      ! 
    136139      INTEGER, PARAMETER ::   nitermax = 20     
    137       !!----------------------------------------------------------------------------- 
     140      !!------------------------------------------------------------------- 
    138141      ! controls 
    139142      IF( nn_timing == 1 )   CALL timing_start('icedyn_rdgrft')                                                             ! timing 
     
    143146         IF(lwp) WRITE(numout,*) 
    144147         IF(lwp) WRITE(numout,*)'ice_dyn_rdgrft: ice ridging and rafting' 
    145          IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~' 
    146          ! 
    147          zdrho = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE 
    148          ! 
     148         IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~~' 
    149149      ENDIF       
    150150 
    151       !-----------------------------------------------------------------------------! 
    152       ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 
    153       !-----------------------------------------------------------------------------! 
    154       ! 
    155       CALL rdgrft_prep                             ! prepare ridging 
    156       ! 
    157       DO jj = 1, jpj                                        ! Initialize arrays. 
     151      CALL ice_var_zapsmall   ! Zero out categories with very small areas 
     152 
     153      !-------------------------------- 
     154      ! 0) Identify grid cells with ice 
     155      !-------------------------------- 
     156      nidx = 0   ;   idxice(:) = 0 
     157      DO jj = 1, jpj 
    158158         DO ji = 1, jpi 
     159            IF ( SUM(a_i(ji,jj,:)) > 0._wp ) THEN 
     160               nidx           = nidx + 1 
     161               idxice( nidx ) = (jj - 1) * jpi + ji 
     162            ENDIF 
     163         END DO 
     164      END DO 
     165 
     166      IF( nidx > 0 ) THEN 
     167         
     168         ! just needed here 
     169         CALL tab_2d_1d( nidx, idxice(1:nidx), zdivu(1:nidx), divu_i(:,:) ) 
     170         CALL tab_2d_1d( nidx, idxice(1:nidx), zdelt(1:nidx), delta_i(:,:) ) 
     171         ! needed here and in the iteration loop 
     172         CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i(:,:,:) ) 
     173         CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i(:,:,:) ) 
     174         CALL tab_2d_1d( nidx, idxice(1:nidx), ato_i_1d(1:nidx)      , ato_i(:,:) ) 
     175 
     176         DO ji = 1, nidx 
    159177            !-----------------------------------------------------------------------------! 
    160178            ! 2) Dynamical inputs (closing rate, divu_adv, opning) 
     
    176194            !  (thick, newly ridged ice). 
    177195 
    178             closing_net(ji,jj) = rn_csrdg * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 
     196            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 
    179197 
    180198            ! 2.2 divu_adv 
     
    186204            ! If divu_adv < 0, make sure the closing rate is large enough 
    187205            ! to give asum = 1.0 after ridging. 
    188              
    189             divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep 
    190  
    191             IF( divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) 
     206            zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 
     207 
     208            IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) ) 
    192209 
    193210            ! 2.3 opning 
    194211            !------------ 
    195212            ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging. 
    196             opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj) 
    197          END DO 
    198       END DO 
    199  
    200       !-----------------------------------------------------------------------------! 
    201       ! 3) Ridging iteration 
    202       !-----------------------------------------------------------------------------! 
    203       niter           = 1                 ! iteration counter 
    204       iterate_ridging = 1 
    205  
    206       DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax ) 
    207  
    208          ! 3.1 closing_gross 
     213            opning(ji) = closing_net(ji) + zdivu_adv(ji) 
     214         END DO 
     215 
     216         ! 
     217         !------------------------------------------------ 
     218         ! 2.1) Identify grid cells with nonzero ridging 
     219         !------------------------------------------------ 
     220         CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d ) 
     221         nidx2 = 0 ; idxice2(:) = 0 
     222         DO ji = 1, nidx 
     223            IF( SUM( ABS(apartf(ji,1:jpl)) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 
     224               nidx2 = nidx2 + 1 
     225               idxice2(nidx2) = idxice(ji) 
     226               a_i_2d(nidx2,:) = a_i_2d(ji,:)  ! adjust to new indices 
     227               v_i_2d(nidx2,:) = v_i_2d(ji,:)  ! adjust to new indices 
     228               ato_i_1d(nidx2) = ato_i_1d(ji)  ! adjust to new indices 
     229            ENDIF 
     230         END DO 
     231         idxice(:) = idxice2(:) 
     232         nidx      = nidx2 
     233 
     234      ENDIF 
     235 
     236      !-------------- 
     237      ! Start ridging 
     238      !-------------- 
     239      IF( nidx > 0 ) THEN 
     240          
     241         !----------- 
     242         ! 2D <==> 1D 
     243         !----------- 
     244         ! fields used but not modified 
     245         CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d(1:nidx), sss_m(:,:) ) 
     246         CALL tab_2d_1d( nidx, idxice(1:nidx), sst_1d(1:nidx), sst_m(:,:) ) 
     247         ! the following fields are modified in this routine 
     248         !!CALL tab_2d_1d( nidx, idxice(1:nidx), ato_i_1d(1:nidx), ato_i(:,:) ) 
     249         !!CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d(1:nidx,1:jpl), a_i(:,:,:) ) 
     250         !!CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i  (:,:,:) ) 
     251         CALL tab_3d_2d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s  (:,:,:) ) 
     252         CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) ) 
     253         CALL tab_3d_2d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i (:,:,:) ) 
     254         IF ( nn_pnd_scheme > 0 ) THEN 
     255            CALL tab_3d_2d( nidx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip(:,:,:) ) 
     256            CALL tab_3d_2d( nidx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip(:,:,:) ) 
     257         ENDIF 
     258         DO jl = 1, jpl 
     259            DO jk = 1, nlay_s 
     260               CALL tab_2d_1d( nidx, idxice(1:nidx), ze_s_2d(1:nidx,jk,jl), e_s(:,:,jk,jl) ) 
     261            END DO 
     262            DO jk = 1, nlay_i 
     263               CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) ) 
     264            END DO 
     265         END DO 
     266         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_dyn_1d    (1:nidx), sfx_dyn    (:,:) ) 
     267         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_bri_1d    (1:nidx), sfx_bri    (:,:) ) 
     268         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_dyn_1d    (1:nidx), wfx_dyn    (:,:) ) 
     269         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_dyn_1d    (1:nidx), hfx_dyn    (:,:) ) 
     270         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_snw_dyn_1d(1:nidx), wfx_snw_dyn(:,:) ) 
     271         IF ( nn_pnd_scheme > 0 ) THEN 
     272            CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_pnd_1d(1:nidx), wfx_pnd(:,:) ) 
     273         ENDIF 
     274          
    209275         !-----------------------------------------------------------------------------! 
    210          ! Based on the ITD of ridging and ridged ice, convert the net 
    211          !  closing rate to a gross closing rate.   
    212          ! NOTE: 0 < aksum <= 1 
    213          closing_gross(:,:) = closing_net(:,:) / aksum(:,:) 
    214  
    215          ! correction to closing rate and opening if closing rate is excessive 
    216          !--------------------------------------------------------------------- 
    217          ! Reduce the closing rate if more than 100% of the open water  
    218          ! would be removed.  Reduce the opening rate proportionately. 
    219          DO jj = 1, jpj 
    220             DO ji = 1, jpi 
    221                zfac   = ( opning(ji,jj) - apartf(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice 
    222                IF    ( zfac < 0._wp .AND. zfac > - ato_i(ji,jj) ) THEN                  ! would lead to negative ato_i 
    223                   opning(ji,jj) = apartf(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice  
    224                ELSEIF( zfac > 0._wp .AND. zfac > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN  ! would lead to ato_i > asum 
    225                   opning(ji,jj) = apartf(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice  
     276         ! 3) Ridging iteration 
     277         !-----------------------------------------------------------------------------! 
     278         niter           = 1 
     279         iterate_ridging = 1       
     280         DO WHILE( iterate_ridging > 0 .AND. niter < nitermax ) 
     281 
     282            CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d ) 
     283 
     284            ! 3.2 Redistribute area, volume, and energy. 
     285            !-----------------------------------------------------------------------------! 
     286            CALL rdgrft_shift 
     287 
     288            ! 3.4 Do we keep on iterating? 
     289            !-----------------------------------------------------------------------------! 
     290            ! Check whether asum = 1.  If not (because the closing and opening 
     291            ! rates were reduced above), ridge again with new rates. 
     292            iterate_ridging = 0 
     293            DO ji = 1, nidx 
     294               zfac = 1._wp - ( ato_i_1d(ji) + SUM( a_i_2d(ji,:) ) ) 
     295               IF( ABS( zfac ) < epsi10 ) THEN 
     296                  closing_net(ji) = 0._wp 
     297                  opning     (ji) = 0._wp 
     298                  ato_i_1d   (ji) = MAX( 0._wp, 1._wp - SUM( a_i_2d(ji,:) ) ) 
     299               ELSE 
     300                  iterate_ridging  = 1 
     301                  zdivu_adv  (ji) = zfac * r1_rdtice 
     302                  closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) ) 
     303                  opning     (ji) = MAX( 0._wp,  zdivu_adv(ji) ) 
    226304               ENDIF 
    227305            END DO 
    228          END DO 
    229  
    230          ! correction to closing rate / opening if excessive ice removal 
    231          !--------------------------------------------------------------- 
    232          ! Reduce the closing rate if more than 100% of any ice category  
    233          ! would be removed.  Reduce the opening rate proportionately. 
     306            ! 
     307            niter = niter + 1 
     308            IF( niter  >  nitermax )    CALL ctl_warn( 'icedyn_rdgrft: non-converging ridging scheme' ) 
     309            ! 
     310         END DO 
     311 
     312         !----------- 
     313         ! 1D <==> 2D 
     314         !----------- 
     315         CALL tab_1d_2d( nidx, idxice(1:nidx), ato_i_1d(1:nidx), ato_i(:,:) ) 
     316         CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d(1:nidx,1:jpl), a_i(:,:,:) ) 
     317         CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i  (:,:,:) ) 
     318         CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s  (:,:,:) ) 
     319         CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) ) 
     320         CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i (:,:,:) ) 
     321         IF ( nn_pnd_scheme > 0 ) THEN 
     322            CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip(:,:,:) ) 
     323            CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip(:,:,:) ) 
     324         ENDIF 
    234325         DO jl = 1, jpl 
    235             DO jj = 1, jpj 
    236                DO ji = 1, jpi 
    237                   zfac = apartf(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 
    238                   IF( zfac  >  a_i(ji,jj,jl) ) THEN 
    239                      closing_gross(ji,jj) = closing_gross(ji,jj) * a_i(ji,jj,jl) / zfac 
    240                   ENDIF 
    241                END DO 
    242             END DO 
    243          END DO 
    244  
    245          ! 3.2 Redistribute area, volume, and energy. 
    246          !-----------------------------------------------------------------------------! 
    247          CALL rdgrft_shift( opning, closing_gross ) 
    248           
    249          ! 3.3 Compute total area of ice plus open water after ridging. 
    250          !-----------------------------------------------------------------------------! 
    251          ! This is in general not equal to one because of divergence during transport 
    252          asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 
    253  
    254          ! 3.4 Do we keep on iterating? 
    255          !-----------------------------------------------------------------------------! 
    256          ! Check whether asum = 1.  If not (because the closing and opening 
    257          ! rates were reduced above), ridge again with new rates. 
    258          iterate_ridging = 0 
    259          DO jj = 1, jpj 
    260             DO ji = 1, jpi 
    261                IF( ABS( asum(ji,jj) - 1._wp ) < epsi10 ) THEN 
    262                   closing_net(ji,jj) = 0._wp 
    263                   opning     (ji,jj) = 0._wp 
    264                   ato_i      (ji,jj) = MAX( 0._wp, 1._wp - SUM( a_i(ji,jj,:) ) ) 
    265                ELSE 
    266                   iterate_ridging    = 1 
    267                   divu_adv   (ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice 
    268                   closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) ) 
    269                   opning     (ji,jj) = MAX( 0._wp,  divu_adv(ji,jj) ) 
    270                ENDIF 
    271             END DO 
    272          END DO 
    273          IF( lk_mpp )   CALL mpp_max( iterate_ridging ) 
    274  
    275          ! Repeat if necessary. 
    276          ! NOTE: If strength smoothing is turned on, the ridging must be 
    277          ! iterated globally because of the boundary update in the smoothing. 
    278          niter = niter + 1 
    279          ! 
    280          IF( iterate_ridging == 1 ) THEN 
    281             CALL rdgrft_prep 
    282             IF( niter  >  nitermax ) THEN 
    283                WRITE(numout,*) ' ALERTE : non-converging ridging scheme ' 
    284                WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging 
    285             ENDIF 
     326            DO jk = 1, nlay_s 
     327               CALL tab_1d_2d( nidx, idxice(1:nidx), ze_s_2d(1:nidx,jk,jl), e_s(:,:,jk,jl) ) 
     328            END DO 
     329            DO jk = 1, nlay_i 
     330               CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) ) 
     331            END DO 
     332         END DO 
     333         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_dyn_1d    (1:nidx), sfx_dyn    (:,:) ) 
     334         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_bri_1d    (1:nidx), sfx_bri    (:,:) ) 
     335         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_dyn_1d    (1:nidx), wfx_dyn    (:,:) ) 
     336         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_dyn_1d    (1:nidx), hfx_dyn    (:,:) ) 
     337         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_snw_dyn_1d(1:nidx), wfx_snw_dyn(:,:) ) 
     338         IF ( nn_pnd_scheme > 0 ) THEN 
     339            CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_pnd_1d(1:nidx), wfx_pnd(:,:) ) 
    286340         ENDIF 
    287          ! 
    288       END DO !! on the do while over iter 
    289  
     341 
     342      ENDIF ! nidx > 0 
     343    
    290344      CALL ice_var_agg( 1 )  
    291345 
     
    298352 
    299353 
    300    SUBROUTINE rdgrft_prep 
    301       !!---------------------------------------------------------------------! 
     354   SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i ) 
     355      !!------------------------------------------------------------------- 
    302356      !!                ***  ROUTINE rdgrft_prep *** 
    303357      !! 
    304       !! ** Purpose :   preparation for ridging and strength calculations 
     358      !! ** Purpose :   preparation for ridging calculations 
    305359      !! 
    306360      !! ** Method  :   Compute the thickness distribution of the ice and open water  
    307361      !!              participating in ridging and of the resulting ridges. 
    308       !!---------------------------------------------------------------------! 
    309       INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    310       REAL(wp) ::   z1_gstar, z1_astar, zhmean, zdummy   ! local scalar 
    311       REAL(wp), DIMENSION(jpi,jpj,-1:jpl) ::   zGsum     ! zGsum(n) = sum of areas in categories 0 to n 
    312       !------------------------------------------------------------------------------! 
    313  
    314       z1_gstar      = 1._wp / rn_gstar     
    315       z1_astar      = 1._wp / rn_astar     
    316  
    317       CALL ice_var_zapsmall   ! Zero out categories with very small areas 
     362      !!------------------------------------------------------------------- 
     363      REAL(wp), DIMENSION(:)  , INTENT(in) ::   pato_i  
     364      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pa_i, pv_i  
     365      !! 
     366      INTEGER  ::   ji, jl                     ! dummy loop indices 
     367      REAL(wp) ::   z1_gstar, z1_astar, zhmean, zfac   ! local scalar 
     368      REAL(wp), DIMENSION(jpij)        ::   zasum, z1_asum, zaksum   ! sum of a_i+ato_i and reverse  
     369      REAL(wp), DIMENSION(jpij,jpl)    ::   zhi                      ! ice thickness 
     370      REAL(wp), DIMENSION(jpij,-1:jpl) ::   zGsum                    ! zGsum(n) = sum of areas in categories 0 to n 
     371      !-------------------------------------------------------------------- 
     372 
     373      z1_gstar = 1._wp / rn_gstar 
     374      z1_astar = 1._wp / rn_astar 
    318375 
    319376      !                       ! Ice thickness needed for rafting 
    320       WHERE( a_i(:,:,:) >= epsi20 )   ;   ht_i(:,:,:) = v_i (:,:,:) / a_i(:,:,:) 
    321       ELSEWHERE                       ;   ht_i(:,:,:) = 0._wp 
     377      WHERE( pa_i(1:nidx,:) > 0._wp )   ;   zhi(1:nidx,:) = pv_i(1:nidx,:) / pa_i(1:nidx,:) 
     378      ELSEWHERE                         ;   zhi(1:nidx,:) = 0._wp 
    322379      END WHERE 
    323380 
     
    338395      ! Compute total area of ice plus open water. 
    339396      ! This is in general not equal to one because of divergence during transport 
    340       asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 
    341       ! 
     397      zasum(1:nidx) = pato_i(1:nidx) + SUM( pa_i(1:nidx,:), dim=2 ) 
     398      ! 
     399      WHERE( zasum(1:nidx) > 0._wp )   ;   z1_asum(1:nidx) = 1._wp / zasum(1:nidx) 
     400      ELSEWHERE                        ;   z1_asum(1:nidx) = 0._wp 
     401      END WHERE 
    342402      ! Compute cumulative thickness distribution function 
    343403      ! Compute the cumulative thickness distribution function zGsum, 
    344404      ! where zGsum(n) is the fractional area in categories 0 to n. 
    345405      ! initial value (in h = 0) equals open water area 
    346       zGsum(:,:,-1) = 0._wp 
    347       zGsum(:,:,0 ) = ato_i(:,:) / asum(:,:) 
     406      zGsum(1:nidx,-1) = 0._wp 
     407      zGsum(1:nidx,0 ) = pato_i(1:nidx) * z1_asum(1:nidx) 
    348408      DO jl = 1, jpl 
    349          zGsum(:,:,jl) = ( ato_i(:,:) + SUM( a_i(:,:,1:jl), dim=3 ) ) / asum(:,:) ! sum(1:jl) is ok (and not jpl) 
     409         zGsum(1:nidx,jl) = ( pato_i(1:nidx) + SUM( pa_i(1:nidx,1:jl), dim=2 ) ) * z1_asum(1:nidx) ! sum(1:jl) is ok (and not jpl) 
    350410      END DO 
    351  
    352411      ! 
    353412      IF( ln_partf_lin ) THEN          !--- Linear formulation (Thorndike et al., 1975) 
    354413         DO jl = 0, jpl     
    355             DO jj = 1, jpj  
    356                DO ji = 1, jpi 
    357                   IF    ( zGsum(ji,jj,jl)   < rn_gstar ) THEN 
    358                      apartf(ji,jj,jl) = z1_gstar * ( zGsum(ji,jj,jl) - zGsum(ji,jj,jl-1) ) * & 
    359                         &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + zGsum(ji,jj,jl) ) * z1_gstar ) 
    360                   ELSEIF( zGsum(ji,jj,jl-1) < rn_gstar ) THEN 
    361                      apartf(ji,jj,jl) = z1_gstar * ( rn_gstar        - zGsum(ji,jj,jl-1) ) *  & 
    362                         &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + rn_gstar        ) * z1_gstar ) 
    363                   ELSE 
    364                      apartf(ji,jj,jl) = 0._wp 
    365                   ENDIF 
    366                END DO 
     414            DO ji = 1, nidx 
     415               IF    ( zGsum(ji,jl)   < rn_gstar ) THEN 
     416                  apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * & 
     417                     &                       ( 2._wp - ( zGsum(ji,jl-1) + zGsum(ji,jl) ) * z1_gstar ) 
     418               ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN 
     419                  apartf(ji,jl) = z1_gstar * ( rn_gstar     - zGsum(ji,jl-1) ) *  & 
     420                     &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar        ) * z1_gstar ) 
     421               ELSE 
     422                  apartf(ji,jl) = 0._wp 
     423               ENDIF 
    367424            END DO 
    368425         END DO 
     
    370427      ELSEIF( ln_partf_exp ) THEN      !--- Exponential, more stable formulation (Lipscomb et al, 2007) 
    371428         !                         
    372          zdummy = 1._wp / ( 1._wp - EXP(-z1_astar) )        ! precompute exponential terms using zGsum as a work array 
     429         zfac = 1._wp / ( 1._wp - EXP(-z1_astar) ) 
    373430         DO jl = -1, jpl 
    374             zGsum(:,:,jl) = EXP( -zGsum(:,:,jl) * z1_astar ) * zdummy 
     431            DO ji = 1, nidx 
     432               zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac 
     433            END DO 
    375434         END DO 
    376435         DO jl = 0, jpl 
    377             apartf(:,:,jl) = zGsum(:,:,jl-1) - zGsum(:,:,jl) 
     436            DO ji = 1, nidx 
     437               apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl) 
     438            END DO 
    378439         END DO 
    379440         ! 
     
    383444      IF( ln_rafting .AND. ln_ridging ) THEN             !- ridging & rafting 
    384445         DO jl = 1, jpl 
    385             DO jj = 1, jpj  
    386                DO ji = 1, jpi 
    387                   aridge(ji,jj,jl) = ( 1._wp + TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jj,jl) 
    388                   araft (ji,jj,jl) = apartf(ji,jj,jl) - aridge(ji,jj,jl) 
    389                END DO 
     446            DO ji = 1, nidx 
     447               aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl) 
     448               araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl) 
    390449            END DO 
    391450         END DO 
    392451      ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone 
    393          aridge(:,:,:) = apartf(:,:,1:jpl) 
    394          araft (:,:,:) = 0._wp 
     452         DO jl = 1, jpl 
     453            DO ji = 1, nidx 
     454               aridge(ji,jl) = apartf(ji,jl) 
     455               araft (ji,jl) = 0._wp 
     456            END DO 
     457         END DO 
    395458      ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone    
    396         aridge(:,:,:) = 0._wp 
    397          araft (:,:,:) = apartf(:,:,1:jpl) 
     459         DO jl = 1, jpl 
     460            DO ji = 1, nidx 
     461               aridge(ji,jl) = 0._wp 
     462               araft (ji,jl) = apartf(ji,jl) 
     463            END DO 
     464         END DO 
    398465      ELSE                                               !- no ridging & no rafting 
    399          aridge(:,:,:) = 0._wp 
    400          araft (:,:,:) = 0._wp          
     466         DO jl = 1, jpl 
     467            DO ji = 1, nidx 
     468               aridge(ji,jl) = 0._wp 
     469               araft (ji,jl) = 0._wp          
     470            END DO 
     471         END DO 
    401472      ENDIF 
    402473 
     
    422493      ! ridging. 
    423494      ! 
    424       ! aksum = net area removed/ total area removed 
     495      ! zaksum = net area removed/ total area removed 
    425496      ! where total area removed = area of ice that ridges 
    426497      !         net area removed = total area removed - area of new ridges 
    427498      !----------------------------------------------------------------- 
    428499 
    429       aksum(:,:) = apartf(:,:,0) 
    430       zdummy = 1._wp / hi_hrft 
     500      zaksum(:) = apartf(:,0) 
     501      zfac = 1._wp / hi_hrft 
    431502      ! Transfer function 
    432503      DO jl = 1, jpl !all categories have a specific transfer function 
    433          DO jj = 1, jpj 
    434             DO ji = 1, jpi 
    435                IF ( apartf(ji,jj,jl) > 0._wp ) THEN 
    436                   zhmean            = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * hrdg_hi_min ) 
    437                   hrmin  (ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( zhmean + ht_i(ji,jj,jl) ) ) 
    438                   hrmax  (ji,jj,jl) = 2._wp * zhmean - hrmin(ji,jj,jl) 
    439                   hraft  (ji,jj,jl) = ht_i(ji,jj,jl) * zdummy 
    440                   hi_hrdg(ji,jj,jl) = ht_i(ji,jj,jl) / MAX( zhmean, epsi20 ) 
    441                   ! 
    442                   ! Normalization factor : aksum, ensures mass conservation 
    443                   aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - hi_hrdg(ji,jj,jl) )    & 
    444                      &                        + araft (ji,jj,jl) * ( 1._wp - hi_hrft ) 
    445                ELSE 
    446                   hrmin  (ji,jj,jl) = 0._wp  
    447                   hrmax  (ji,jj,jl) = 0._wp  
    448                   hraft  (ji,jj,jl) = 0._wp  
    449                   hi_hrdg(ji,jj,jl) = 1._wp 
    450                ENDIF 
    451             END DO 
     504         DO ji = 1, nidx 
     505            IF ( apartf(ji,jl) > 0._wp ) THEN 
     506               zhmean         = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min ) 
     507               hrmin  (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) ) 
     508               hrmax  (ji,jl) = 2._wp * zhmean - hrmin(ji,jl) 
     509               hraft  (ji,jl) = zhi(ji,jl) * zfac 
     510               hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 ) 
     511               ! 
     512               ! Normalization factor : zaksum, ensures mass conservation 
     513               zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) )    & 
     514                  &                    + araft (ji,jl) * ( 1._wp - hi_hrft ) 
     515            ELSE 
     516               hrmin  (ji,jl) = 0._wp  
     517               hrmax  (ji,jl) = 0._wp  
     518               hraft  (ji,jl) = 0._wp  
     519               hi_hrdg(ji,jl) = 1._wp 
     520            ENDIF 
    452521         END DO 
    453522      END DO 
    454523      ! 
     524      ! 3.1 closing_gross 
     525      !-----------------------------------------------------------------------------! 
     526      ! Based on the ITD of ridging and ridged ice, convert the net 
     527      !  closing rate to a gross closing rate.   
     528      ! NOTE: 0 < aksum <= 1 
     529      WHERE( zaksum(1:nidx) > 0._wp )   ;   closing_gross(1:nidx) = closing_net(1:nidx) / zaksum(1:nidx) 
     530      ELSEWHERE                         ;   closing_gross(1:nidx) = 0._wp 
     531      END WHERE 
     532       
     533      DO ji = 1, nidx 
     534          
     535         ! correction to closing rate and opening if closing rate is excessive 
     536         !--------------------------------------------------------------------- 
     537         ! Reduce the closing rate if more than 100% of the open water  
     538         ! would be removed.  Reduce the opening rate proportionately. 
     539         zfac = ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice 
     540         IF    ( zfac < 0._wp .AND. zfac > - pato_i(ji) ) THEN                  ! would lead to negative ato_i 
     541            opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_rdtice  
     542         ELSEIF( zfac > 0._wp .AND. zfac > ( zasum(ji) - pato_i(ji) ) ) THEN  ! would lead to ato_i > asum 
     543            opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_rdtice  
     544         ENDIF 
     545      END DO 
     546       
     547      ! correction to closing rate / opening if excessive ice removal 
     548      !--------------------------------------------------------------- 
     549      ! Reduce the closing rate if more than 100% of any ice category  
     550      ! would be removed.  Reduce the opening rate proportionately. 
     551      DO jl = 1, jpl 
     552         DO ji = 1, nidx 
     553            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 
     554!!            IF( zfac > pa_i(ji,jl) .AND. zfac > 0._wp ) THEN 
     555            IF( zfac > pa_i(ji,jl) ) THEN 
     556               closing_gross(ji) = closing_gross(ji) * pa_i(ji,jl) / zfac 
     557            ENDIF 
     558         END DO 
     559      END DO       
     560      ! 
    455561   END SUBROUTINE rdgrft_prep 
    456562 
    457563 
    458    SUBROUTINE rdgrft_shift( opning, closing_gross ) 
    459       !!---------------------------------------------------------------------- 
     564   SUBROUTINE rdgrft_shift 
     565      !!------------------------------------------------------------------- 
    460566      !!                ***  ROUTINE rdgrft_shift *** 
    461567      !! 
     
    464570      !! ** Method  :   Remove area, volume, and energy from each ridging category 
    465571      !!              and add to thicker ice categories. 
    466       !!---------------------------------------------------------------------- 
    467       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   opning         ! rate of opening due to divergence/shear 
    468       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   closing_gross  ! rate at which area retreats, excluding area of new ridges 
     572      !!------------------------------------------------------------------- 
    469573      ! 
    470574      INTEGER  ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices 
    471       INTEGER  ::   ij                         ! horizontal index, combines i and j loops 
    472       INTEGER  ::   icells                     ! number of cells with a_i > puny 
    473575      REAL(wp) ::   hL, hR, farea              ! left and right limits of integration and new area going to jl2 
    474  
    475       INTEGER , DIMENSION(jpij) ::   indxi, indxj     ! compressed indices 
     576      REAL(wp) ::   vsw                        ! vol of water trapped into ridges 
     577      REAL(wp) ::   afrdg, afrft               ! fraction of category area ridged/rafted  
     578      REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1 
     579      REAL(wp)                  ::   airft1, oirft1, aprft1 
     580      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, esrdg  ! area etc of new ridges 
     581      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, esrft  ! area etc of rafted ice 
     582      ! 
     583      REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges 
    476584      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2 
    477  
    478       REAL(wp), DIMENSION(jpij) ::   afrac            ! fraction of category area ridged  
    479       REAL(wp), DIMENSION(jpij) ::   ardg1 , ardg2    ! area of ice ridged & new ridges 
    480       REAL(wp), DIMENSION(jpij) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice 
    481       ! MV MP 2016 
    482       REAL(wp), DIMENSION(jpij) ::   vprdg            ! pond volume of ridging ice 
    483       REAL(wp), DIMENSION(jpij) ::   aprdg1           ! pond area of ridging ice 
    484       REAL(wp), DIMENSION(jpij) ::   aprdg2           ! pond area of ridging ice 
    485       ! END MV MP 2016 
    486       REAL(wp), DIMENSION(jpij) ::   dhr, dhr2        ! hrmax - hrmin  &  hrmax^2 - hrmin^2 
    487  
    488       REAL(wp), DIMENSION(jpij) ::   vrdg1            ! volume of ice ridged 
    489       REAL(wp), DIMENSION(jpij) ::   vrdg2            ! volume of new ridges 
    490       REAL(wp), DIMENSION(jpij) ::   vsw              ! volume of seawater trapped into ridges 
    491       REAL(wp), DIMENSION(jpij) ::   srdg1            ! sal*volume of ice ridged 
    492       REAL(wp), DIMENSION(jpij) ::   srdg2            ! sal*volume of new ridges 
    493       REAL(wp), DIMENSION(jpij) ::   smsw             ! sal*volume of water trapped into ridges 
    494       REAL(wp), DIMENSION(jpij) ::   oirdg1, oirdg2   ! ice age of ice ridged 
    495  
    496       REAL(wp), DIMENSION(jpij) ::   afrft            ! fraction of category area rafted 
    497       REAL(wp), DIMENSION(jpij) ::   arft1 , arft2    ! area of ice rafted and new rafted zone 
    498       REAL(wp), DIMENSION(jpij) ::   virft , vsrft    ! ice & snow volume of rafting ice 
    499       ! MV MP 2016 
    500       REAL(wp), DIMENSION(jpij) ::   vprft            ! pond volume of rafting ice 
    501       REAL(wp), DIMENSION(jpij) ::   aprft1           ! pond area of rafted ice 
    502       REAL(wp), DIMENSION(jpij) ::   aprft2           ! pond area of new rafted ice 
    503       ! END MV MP 2016 
    504       REAL(wp), DIMENSION(jpij) ::   esrft , smrft    ! snow energy & salinity of rafting ice 
    505       REAL(wp), DIMENSION(jpij) ::   oirft1, oirft2   ! ice age of ice rafted 
    506  
     585      ! 
     586      REAL(wp), DIMENSION(jpij,jpl) ::   z1_ai        ! 1 / a 
     587      ! 
    507588      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft     ! ice energy of rafting ice 
    508       REAL(wp), DIMENSION(jpij,nlay_i) ::   erdg1     ! enth*volume of ice ridged 
    509       REAL(wp), DIMENSION(jpij,nlay_i) ::   erdg2     ! enth*volume of new ridges 
    510       REAL(wp), DIMENSION(jpij,nlay_i) ::   ersw      ! enth of water trapped into ridges 
    511       !!---------------------------------------------------------------------- 
     589      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirdg     ! enth*volume of new ridges 
     590      !!------------------------------------------------------------------- 
    512591 
    513592      !------------------------------------------------------------------------------- 
    514593      ! 1) Compute change in open water area due to closing and opening. 
    515594      !------------------------------------------------------------------------------- 
    516       DO jj = 1, jpj 
    517          DO ji = 1, jpi 
    518             ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) +  & 
    519                &                     ( opning(ji,jj) - apartf(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice ) 
    520          END DO 
     595      DO ji = 1, nidx 
     596         ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice ) 
    521597      END DO 
    522  
     598       
    523599      !----------------------------------------------------------------- 
    524       ! 2) Pump everything from ice which is being ridged / rafted 
     600      ! 2) compute categories participating to ridging/rafting (jl1)  
    525601      !----------------------------------------------------------------- 
    526       ! Compute the area, volume, and energy of ice ridging in each 
    527       ! category, along with the area of the resulting ridge. 
    528  
    529       DO jl1 = 1, jpl !jl1 describes the ridging category 
    530  
    531          !------------------------------------------------ 
    532          ! 2.1) Identify grid cells with nonzero ridging 
    533          !------------------------------------------------ 
    534          icells = 0 
    535          DO jj = 1, jpj 
    536             DO ji = 1, jpi 
    537                IF( apartf(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN 
    538                   icells = icells + 1 
    539                   indxi(icells) = ji 
    540                   indxj(icells) = jj 
     602      DO jl1 = 1, jpl 
     603 
     604         CALL tab_2d_1d( nidx, idxice(1:nidx), sm_i_1d(1:nidx), sm_i(:,:,jl1) ) 
     605 
     606         DO ji = 1, nidx 
     607 
     608            !------------------------------------------------ 
     609            ! 2.1) Identify grid cells with nonzero ridging 
     610            !------------------------------------------------ 
     611            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 
     612 
     613               z1_ai(ji,jl1) = 1._wp / a_i_2d(ji,jl1) 
     614 
     615               !-------------------------------------------------------------------- 
     616               ! 2.2) Compute area of ridging ice (airdg1) and of new ridge (airdg2) 
     617               !-------------------------------------------------------------------- 
     618               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice 
     619               airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice 
     620 
     621               airdg2(ji) = airdg1 * hi_hrdg(ji,jl1) 
     622               airft2(ji) = airft1 * hi_hrft 
     623 
     624               !--------------------------------------------------------------- 
     625               ! 2.3) Compute ridging /rafting fractions, make sure afrdg <=1  
     626               !--------------------------------------------------------------- 
     627               afrdg = airdg1 * z1_ai(ji,jl1) 
     628               afrft = airft1 * z1_ai(ji,jl1) 
     629 
     630               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges 
     631               vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg 
     632               ersw(ji) = -rhoic * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?) 
     633 
     634               !--------------------------------------------------------------- 
     635               ! 2.4) Compute ridging ice and new ridges (vi, vs, sm, oi, es, ei) 
     636               !--------------------------------------------------------------- 
     637               virdg1     = v_i_2d  (ji,jl1)   * afrdg 
     638               virdg2(ji) = v_i_2d  (ji,jl1)   * afrdg * ( 1. + rn_porordg ) 
     639               vsrdg(ji)  = v_s_2d  (ji,jl1)   * afrdg 
     640               sirdg1     = smv_i_2d(ji,jl1)   * afrdg 
     641               sirdg2(ji) = smv_i_2d(ji,jl1)   * afrdg + vsw * sss_1d(ji) 
     642               oirdg1     = oa_i_2d (ji,jl1)   * afrdg 
     643               oirdg2(ji) = oa_i_2d (ji,jl1)   * afrdg * hi_hrdg(ji,jl1)  
     644               esrdg(ji)  = ze_s_2d (ji,1,jl1) * afrdg 
     645 
     646               virft(ji)  = v_i_2d  (ji,jl1)   * afrft 
     647               vsrft(ji)  = v_s_2d  (ji,jl1)   * afrft 
     648               sirft(ji)  = smv_i_2d(ji,jl1)   * afrft  
     649               oirft1     = oa_i_2d (ji,jl1)   * afrft  
     650               oirft2(ji) = oa_i_2d (ji,jl1)   * afrft * hi_hrft  
     651               esrft(ji)  = ze_s_2d (ji,1,jl1) * afrft 
     652 
     653               !MV MP 2016 
     654               IF ( nn_pnd_scheme > 0 ) THEN 
     655                  aprdg1     = a_ip_2d(ji,jl1) * afrdg 
     656                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) 
     657                  vprdg (ji) = v_ip_2d(ji,jl1) * afrdg 
     658                  aprft1     = a_ip_2d(ji,jl1) * afrft 
     659                  aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft 
     660                  vprft (ji) = v_ip_2d(ji,jl1) * afrft 
    541661               ENDIF 
    542             END DO 
    543          END DO 
    544  
    545          DO ij = 1, icells 
    546             ji = indxi(ij) ; jj = indxj(ij) 
    547  
    548             !-------------------------------------------------------------------- 
    549             ! 2.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2) 
    550             !-------------------------------------------------------------------- 
    551             ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice 
    552             arft1(ij) = araft (ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice 
    553  
    554             !--------------------------------------------------------------- 
    555             ! 2.3) Compute ridging /rafting fractions, make sure afrac <=1  
    556             !--------------------------------------------------------------- 
    557             afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging 
    558             afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting 
    559             ardg2(ij) = ardg1(ij) * hi_hrdg(ji,jj,jl1) 
    560             arft2(ij) = arft1(ij) * hi_hrft 
    561  
    562             !-------------------------------------------------------------------------- 
    563             ! 2.4) Substract area, volume, and energy from ridging  
    564             !     / rafting category n1. 
    565             !-------------------------------------------------------------------------- 
    566             vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij) 
    567             vrdg2(ij) = vrdg1(ij) * ( 1. + rn_porordg ) 
    568             vsw  (ij) = vrdg1(ij) * rn_porordg 
    569  
    570             vsrdg (ij) = v_s  (ji,jj,  jl1) * afrac(ij) 
    571             esrdg (ij) = e_s  (ji,jj,1,jl1) * afrac(ij) 
    572             !MV MP 2016 
    573             IF ( nn_pnd_scheme > 0 ) THEN 
    574                aprdg1(ij) = a_ip(ji,jj, jl1) * afrac(ij) 
    575                aprdg2(ij) = a_ip(ji,jj, jl1) * afrac(ij) * hi_hrdg(ji,jj,jl1) 
    576                vprdg(ij)  = v_ip(ji,jj, jl1) * afrac(ij) 
     662               ! END MV MP 2016 
     663 
     664               !----------------------------------------------------------------- 
     665               ! 2.5) Ice-ocean exchanges 
     666               !----------------------------------------------------------------- 
     667               wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids 
     668               sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoic * r1_rdtice 
     669               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_rdtice          ! > 0 [W.m-2]  
     670 
     671               ! Put the snow lost by ridging into the ocean 
     672               !------------------------------------------------ 
     673               !  Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow. 
     674               wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhosn * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean 
     675                  &                                      + rhosn * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice 
     676 
     677               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji) * ( 1._wp - rn_fsnwrdg )   &                 ! heat sink for ocean (<0, W.m-2) 
     678                  &                                - esrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice 
     679 
     680               ! Put the melt pond water into the ocean 
     681               !------------------------------------------             
     682               IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN 
     683                  wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean 
     684                     &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice 
     685               ENDIF 
     686 
     687               ! virtual salt flux to keep salinity constant 
     688               !------------------------------------------             
     689               IF( nn_icesal /= 2 )  THEN 
     690                  sirdg2(ji)      = sirdg2(ji)      - vsw * ( sss_1d(ji) - sm_i_1d(ji) )        ! ridge salinity = sm_i 
     691                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji)  * vsw * rhoic * r1_rdtice  &  ! put back sss_m into the ocean 
     692                     &                            - sm_i_1d(ji) * vsw * rhoic * r1_rdtice     ! and get  sm_i  from the ocean  
     693               ENDIF 
     694 
     695 
     696               !------------------------------------------             
     697               ! 2.6) update jl1 (removing ridged/rafted area) 
     698               !------------------------------------------             
     699               a_i_2d  (ji,jl1) = a_i_2d  (ji,jl1) - airdg1    - airft1 
     700               v_i_2d  (ji,jl1) = v_i_2d  (ji,jl1) - virdg1    - virft(ji) 
     701               v_s_2d  (ji,jl1) = v_s_2d  (ji,jl1) - vsrdg(ji) - vsrft(ji) 
     702               smv_i_2d(ji,jl1) = smv_i_2d(ji,jl1) - sirdg1    - sirft(ji) 
     703               oa_i_2d (ji,jl1) = oa_i_2d (ji,jl1) - oirdg1    - oirft1 
     704               ! MV MP 2016 
     705               IF ( nn_pnd_scheme > 0 ) THEN 
     706                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1 
     707                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) 
     708               ENDIF 
     709               ! END MV MP 2016 
     710               ze_s_2d(ji,1,jl1) = ze_s_2d(ji,1,jl1) - esrdg (ji) - esrft (ji) 
    577711            ENDIF 
    578             ! END MV MP 2016 
    579             srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij) 
    580             oirdg1(ij) = oa_i (ji,jj,  jl1) * afrac(ij) 
    581             oirdg2(ij) = oa_i (ji,jj,  jl1) * afrac(ij) * hi_hrdg(ji,jj,jl1)  
    582  
    583             ! rafting volumes, heat contents ... 
    584             virft (ij) = v_i  (ji,jj,  jl1) * afrft(ij) 
    585             vsrft (ij) = v_s  (ji,jj,  jl1) * afrft(ij) 
    586             !MV MP 2016 
    587             IF ( nn_pnd_scheme > 0 ) THEN 
    588                aprft1(ij) = a_ip (ji,jj,  jl1) * afrft(ij) 
    589                aprft2(ij) = a_ip (ji,jj,  jl1) * afrft(ij) * hi_hrft 
    590                vprft(ij)  = v_ip(ji,jj,jl1)    * afrft(ij) 
    591             ENDIF 
    592             ! END MV MP 2016 
    593             srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij) 
    594             esrft (ij) = e_s  (ji,jj,1,jl1) * afrft(ij) 
    595             smrft (ij) = smv_i(ji,jj,  jl1) * afrft(ij)  
    596             oirft1(ij) = oa_i (ji,jj,  jl1) * afrft(ij)  
    597             oirft2(ij) = oa_i (ji,jj,  jl1) * afrft(ij) * hi_hrft  
    598  
    599             !----------------------------------------------------------------- 
    600             ! 2.5) Compute properties of new ridges 
    601             !----------------------------------------------------------------- 
    602             smsw(ij)  = vsw(ij) * sss_m(ji,jj)                   ! salt content of seawater frozen in voids 
    603             srdg2(ij) = srdg1(ij) + smsw(ij)                     ! salt content of new ridge 
    604              
    605             sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ij) * rhoic * r1_rdtice 
    606             wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids 
    607  
    608             ! virtual salt flux to keep salinity constant 
    609             IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN 
    610                srdg2(ij)      = srdg2(ij) - vsw(ij) * ( sss_m(ji,jj) - sm_i(ji,jj,jl1) )           ! ridge salinity = sm_i 
    611                sfx_bri(ji,jj) = sfx_bri(ji,jj) + sss_m(ji,jj)    * vsw(ij) * rhoic * r1_rdtice  &  ! put back sss_m into the ocean 
    612                   &                            - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice     ! and get  sm_i  from the ocean  
    613             ENDIF 
    614  
    615             !------------------------------------------             
    616             ! 2.6 Put the snow somewhere in the ocean 
    617             !------------------------------------------             
    618             !  Place part of the snow lost by ridging into the ocean.  
    619             !  Note that esrdg > 0; the ocean must cool to melt snow. 
    620             !  If the ocean temp = Tf already, new ice must grow. 
    621             !  During the next time step, thermo_rates will determine whether 
    622             !  the ocean cools or new ice grows. 
    623             wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnwrdg )   &  
    624                &                                      + rhosn * vsrft(ij) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice  ! fresh water source for ocean 
    625  
    626             hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnwrdg )         &  
    627                &                                - esrft(ij) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2) 
    628  
    629             ! MV MP 2016 
    630             !------------------------------------------             
    631             ! 2.7 Put the melt pond water in the ocean 
    632             !------------------------------------------             
    633             !  Place part of the melt pond volume into the ocean.  
    634             IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN 
    635                wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + ( rhofw * vprdg(ij) * ( 1._wp - rn_fpndrdg )   &  
    636                &                                 + rhofw * vprft(ij) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice  ! fresh water source for ocean 
    637             ENDIF 
    638             ! END MV MP 2016 
    639  
    640             !----------------------------------------------------------------- 
    641             ! 2.8 Compute quantities used to apportion ice among categories 
    642             ! in the n2 loop below 
    643             !----------------------------------------------------------------- 
    644             dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1)                    - hrmin(ji,jj,jl1)                    ) 
    645             dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) ) 
    646  
    647  
    648             ! update jl1 (removing ridged/rafted area) 
    649             a_i  (ji,jj,  jl1) = a_i  (ji,jj,  jl1) - ardg1 (ij) - arft1 (ij) 
    650             v_i  (ji,jj,  jl1) = v_i  (ji,jj,  jl1) - vrdg1 (ij) - virft (ij) 
    651             v_s  (ji,jj,  jl1) = v_s  (ji,jj,  jl1) - vsrdg (ij) - vsrft (ij) 
    652             e_s  (ji,jj,1,jl1) = e_s  (ji,jj,1,jl1) - esrdg (ij) - esrft (ij) 
    653             smv_i(ji,jj,  jl1) = smv_i(ji,jj,  jl1) - srdg1 (ij) - smrft (ij) 
    654             oa_i (ji,jj,  jl1) = oa_i (ji,jj,  jl1) - oirdg1(ij) - oirft1(ij) 
    655  
    656             ! MV MP 2016 
    657             IF ( nn_pnd_scheme > 0 ) THEN 
    658                v_ip (ji,jj,jl1) = v_ip (ji,jj,jl1) - vprdg (ij) - vprft (ij) 
    659                a_ip (ji,jj,jl1) = a_ip (ji,jj,jl1) - aprdg1(ij) - aprft1(ij) 
    660             ENDIF 
    661             ! END MV MP 2016 
    662  
    663          END DO 
    664  
    665          !-------------------------------------------------------------------- 
    666          ! 2.9 Compute ridging ice enthalpy, remove it from ridging ice and 
    667          !      compute ridged ice enthalpy  
    668          !-------------------------------------------------------------------- 
     712 
     713         END DO ! ji 
     714 
     715         ! special loop for e_i because of layers jk 
    669716         DO jk = 1, nlay_i 
    670             DO ij = 1, icells 
    671                ji = indxi(ij) ; jj = indxj(ij) 
    672                ! heat content of ridged ice 
    673                erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij)  
    674                eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij)                
    675                 
    676                ! enthalpy of the trapped seawater (J/m2, >0) 
    677                ! clem: if sst>0, then ersw <0 (is that possible?) 
    678                ersw(ij,jk)  = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i 
    679  
    680                ! heat flux to the ocean 
    681                hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux  
    682  
    683                ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 
    684                erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk) 
    685  
    686                ! update jl1 
    687                e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk) 
    688  
    689             END DO 
    690          END DO 
     717            DO ji = 1, nidx 
     718               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 
     719                  ! Compute ridging /rafting fractions 
     720                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji,jl1) 
     721                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji,jl1) 
     722                  ! Compute ridging ice and new ridges for ei 
     723                  eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i 
     724                  eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft 
     725                  ! Update jl1 
     726                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft )  
     727               ENDIF 
     728            END DO 
     729         END DO 
     730          
    691731 
    692732         !------------------------------------------------------------------------------- 
     
    694734         !------------------------------------------------------------------------------- 
    695735         DO jl2  = 1, jpl  
    696             ! over categories to which ridged/rafted ice is transferred 
    697             DO ij = 1, icells 
    698                ji = indxi(ij) ; jj = indxj(ij) 
    699  
    700                ! Compute the fraction of ridged ice area and volume going to thickness category jl2. 
    701                IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN 
    702                   hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) ) 
    703                   hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   ) 
    704                   farea    = ( hR      - hL      ) * dhr(ij)  
    705                   fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij) 
    706                ELSE 
    707                   farea    = 0._wp  
    708                   fvol(ij) = 0._wp                   
     736            ! 
     737            DO ji = 1, nidx 
     738 
     739               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 
     740 
     741                  ! Compute the fraction of ridged ice area and volume going to thickness category jl2. 
     742                  IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN 
     743                     !----------------------------------------------------------------- 
     744                     ! 2.8 Compute quantities used to apportion ice among categories 
     745                     ! in the n2 loop below 
     746                     !----------------------------------------------------------------- 
     747                     hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) ) 
     748                     hR = MIN( hrmax(ji,jl1), hi_max(jl2)   ) 
     749                     farea    = ( hR      - hL      ) / ( hrmax(ji,jl1)                 - hrmin(ji,jl1)                 ) 
     750                     fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) ) 
     751                  ELSE 
     752                     farea    = 0._wp  
     753                     fvol(ji) = 0._wp                   
     754                  ENDIF 
     755 
     756                  ! Compute the fraction of rafted ice area and volume going to thickness category jl2 
     757                  !!gm see above               IF( hraft(ji) <= hi_max(jl2) .AND. hraft(ji) >  hi_max(jl2-1) ) THEN 
     758                  IF( hi_max(jl2-1) < hraft(ji,jl1) .AND. hraft(ji,jl1) <= hi_max(jl2)  ) THEN   ;   zswitch(ji) = 1._wp 
     759                  ELSE                                                                           ;   zswitch(ji) = 0._wp 
     760                  ENDIF 
     761                  ! 
     762                  a_i_2d  (ji,jl2) = a_i_2d  (ji,jl2) + ( airdg2(ji) * farea    + airft2(ji) * zswitch(ji) ) 
     763                  oa_i_2d (ji,jl2) = oa_i_2d (ji,jl2) + ( oirdg2(ji) * farea    + oirft2(ji) * zswitch(ji) ) 
     764                  v_i_2d  (ji,jl2) = v_i_2d  (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) ) 
     765                  smv_i_2d(ji,jl2) = smv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) ) 
     766                  v_s_2d  (ji,jl2) = v_s_2d  (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  & 
     767                     &                                    vsrft (ji) * rn_fsnwrft * zswitch(ji) ) 
     768                  ! MV MP 2016 
     769                  IF ( nn_pnd_scheme > 0 ) THEN 
     770                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   & 
     771                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   ) 
     772                     a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         &  
     773                        &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   ) 
     774                  ENDIF 
     775                  ! END MV MP 2016 
     776                  ze_s_2d(ji,1,jl2) = ze_s_2d(ji,1,jl2) + ( esrdg (ji) * rn_fsnwrdg * fvol(ji)  +  & 
     777                     &                                      esrft (ji) * rn_fsnwrft * zswitch(ji) ) 
     778                   
    709779               ENDIF 
    710780 
    711                ! Compute the fraction of rafted ice area and volume going to thickness category jl2 
    712 !!gm see above               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN 
    713                IF( hi_max(jl2-1) < hraft(ji,jj,jl1) .AND. hraft(ji,jj,jl1) <= hi_max(jl2)  ) THEN   ;   zswitch(ij) = 1._wp 
    714                ELSE                                                                                 ;   zswitch(ij) = 0._wp 
    715                ENDIF 
    716                ! 
    717                a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ( ardg2 (ij) * farea    + arft2 (ij) * zswitch(ij) ) 
    718                oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + ( oirdg2(ij) * farea    + oirft2(ij) * zswitch(ij) ) 
    719                v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) ) 
    720                smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) ) 
    721                v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + ( vsrdg (ij) * rn_fsnwrdg * fvol(ij)  +  & 
    722                   &                                        vsrft (ij) * rn_fsnwrft * zswitch(ij) ) 
    723                e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnwrdg * fvol(ij)  +  & 
    724                   &                                        esrft (ij) * rn_fsnwrft * zswitch(ij) ) 
    725                ! MV MP 2016 
    726                IF ( nn_pnd_scheme > 0 ) THEN 
    727                   v_ip (ji,jj,jl2) = v_ip(ji,jj,jl2) + (   vprdg (ij) * rn_fpndrdg * fvol   (ij)   & 
    728                      &                                   + vprft (ij) * rn_fpndrft * zswitch(ij)   ) 
    729                   a_ip (ji,jj,jl2) = a_ip(ji,jj,jl2) + (   aprdg2(ij) * rn_fpndrdg * farea         &  
    730                      &                                   + aprft2(ij) * rn_fpndrft * zswitch(ji)   ) 
    731                ENDIF 
    732                ! END MV MP 2016 
    733             END DO 
    734  
    735             ! Transfer ice energy to category jl2 by ridging 
     781            END DO 
     782 
    736783            DO jk = 1, nlay_i 
    737                DO ij = 1, icells 
    738                   ji = indxi(ij) ; jj = indxj(ij) 
    739                   e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + erdg2(ij,jk) * fvol(ij) + eirft(ij,jk) * zswitch(ij)                   
     784               DO ji = 1, nidx 
     785                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   & 
     786                     &   ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)                   
    740787               END DO 
    741788            END DO 
     
    743790         END DO ! jl2 
    744791         ! 
    745       END DO ! jl1 (deforming categories) 
     792      END DO ! jl1 
     793      ! 
     794      ! In case ridging/rafting lead to very small negative values (sometimes it happens) 
     795      WHERE( a_i_2d(1:nidx,:) < 0._wp )   a_i_2d(1:nidx,:) = 0._wp 
     796      WHERE( v_i_2d(1:nidx,:) < 0._wp )   v_i_2d(1:nidx,:) = 0._wp 
    746797      ! 
    747798   END SUBROUTINE rdgrft_shift 
     
    760811      !!              ice strength has to be smoothed 
    761812      !!---------------------------------------------------------------------- 
    762       INTEGER             ::   ji,jj, jl   ! dummy loop indices 
     813      INTEGER             ::   ji, jj, jl  ! dummy loop indices 
    763814      INTEGER             ::   ismooth     ! smoothing the resistance to deformation 
    764815      INTEGER             ::   itframe     ! number of time steps for the P smoothing 
     
    767818      REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps 
    768819      !!---------------------------------------------------------------------- 
    769  
    770       !                              !--------------------------------------------------! 
    771       CALL rdgrft_prep               ! Thickness distribution of ridging and ridged ice ! 
    772       !                              !--------------------------------------------------! 
    773  
    774820      !                              !--------------------------------------------------! 
    775821      IF( ln_str_R75 ) THEN          ! Ice strength => Rothrock (1975) method           ! 
    776822      !                              !--------------------------------------------------! 
    777          z1_3 = 1._wp / 3._wp 
    778          DO jl = 1, jpl 
    779             WHERE( apartf(:,:,jl) > 0._wp ) 
    780                strength(:,:) =  -         apartf(:,:,jl) * ht_i   (:,:,jl) * ht_i(:,:,jl)   &  ! PE loss from deforming ice 
    781                   &             + 2._wp * araft (:,:,jl) * ht_i   (:,:,jl) * ht_i(:,:,jl)   &  ! PE gain from rafting ice 
    782                   &             +         aridge(:,:,jl) * hi_hrdg(:,:,jl) * z1_3 *         &  ! PE gain from ridging ice 
    783                   &                      ( hrmax(:,:,jl) * hrmax  (:,:,jl) +                & 
    784                   &                        hrmin(:,:,jl) * hrmin  (:,:,jl) +                & 
    785                   &                        hrmax(:,:,jl) * hrmin  (:,:,jl) ) 
    786             ELSEWHERE 
    787                strength(:,:) = 0._wp 
    788             END WHERE 
    789          END DO 
    790          strength(:,:) = rn_perdg * zdrho * strength(:,:) / aksum(:,:) * tmask(:,:,1) 
    791                          ! where zdrho = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_perdg accounts for frictional dissipation 
    792          ismooth = 1 
    793823      !                              !--------------------------------------------------! 
    794824      ELSEIF( ln_str_H79 ) THEN      ! Ice strength => Hibler (1979) method             ! 
     
    805835         DO jj = 2, jpjm1 
    806836            DO ji = 2, jpim1 
    807                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN  
     837               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN  
    808838                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              & 
    809839                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &   
     
    831861         DO jj = 2, jpjm1 
    832862            DO ji = 2, jpim1 
    833                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN  
     863               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN  
    834864                  itframe = 1 ! number of time steps for the running mean 
    835865                  IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1 
Note: See TracChangeset for help on using the changeset viewer.