Changeset 8369


Ignore:
Timestamp:
2017-07-25T16:38:38+02:00 (3 years ago)
Author:
clem
Message:

STEP4 (4): put all thermodynamics in 1D (limitd_th OK)

Location:
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
Files:
8 edited

Legend:

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

    r8361 r8369  
    2525   USE limvar           ! LIM-3 variables 
    2626   USE limcons          ! conservation tests 
     27   USE limtab 
    2728   ! 
    2829   USE prtctl           ! Print control 
     
    3536   PRIVATE 
    3637 
    37    PUBLIC   lim_itd_th_rem 
    38    PUBLIC   lim_itd_th_reb 
     38   PUBLIC   lim_itd_th_rem   ! called in limthd 
     39   PUBLIC   lim_itd_th_reb   ! called in limupdate 
    3940 
    4041   !!---------------------------------------------------------------------- 
     
    5859      INTEGER , INTENT (in) ::   kt      ! Ocean time step  
    5960      ! 
    60       INTEGER  ::   ji, jj, jl     ! dummy loop index 
    61       INTEGER  ::   ii, ij         ! 2D corresponding indices to ji 
    62       INTEGER  ::   nd             ! local integer 
     61      INTEGER  ::   ji, jj, jl, jcat     ! dummy loop index 
     62      INTEGER  ::   nidx2                ! local integer 
    6363      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars 
    6464      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      - 
    6565      REAL(wp) ::   zx3         
    66  
    67       INTEGER , POINTER, DIMENSION(:,:,:) ::   zdonor   ! donor category index 
    68  
    69       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdhice      ! ice thickness increment 
    70       REAL(wp), POINTER, DIMENSION(:,:,:) ::   g0          ! coefficients for fitting the line of the ITD 
    71       REAL(wp), POINTER, DIMENSION(:,:,:) ::   g1          ! coefficients for fitting the line of the ITD 
    72       REAL(wp), POINTER, DIMENSION(:,:,:) ::   hL          ! left boundary for the ITD for each thickness 
    73       REAL(wp), POINTER, DIMENSION(:,:,:) ::   hR          ! left boundary for the ITD for each thickness 
    74       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice          ! local increment of ice area and volume 
    75       INTEGER , DIMENSION(jpij)  ::   nind_i, nind_j          ! compressed indices for i/j directions 
    76       REAL(wp)                            ::   zslope                  ! used to compute local thermodynamic "speeds" 
    77       REAL(wp), POINTER, DIMENSION(:,:)   ::   zhb0, zhb1              ! category boundaries for thinnes categories 
    78       INTEGER , POINTER, DIMENSION(:,:)   ::   zremap_flag      ! compute remapping or not ???? 
    79       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhbnew           ! new boundaries of ice categories 
    80       !!------------------------------------------------------------------ 
    81  
    82       CALL wrk_alloc( jpi,jpj, zremap_flag ) 
    83       CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) 
    84       CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR ) 
    85       CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    86       CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
    87       CALL wrk_alloc( jpi,jpj, zhb0, zhb1 ) 
     66      REAL(wp) ::   zslope          ! used to compute local thermodynamic "speeds" 
     67 
     68      INTEGER , DIMENSION(jpij)       ::   idxice2         ! compute remapping or not 
     69      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor          ! donor category index 
     70      REAL(wp), DIMENSION(jpij,jpl)   ::   zdhice          ! ice thickness increment 
     71      REAL(wp), DIMENSION(jpij,jpl)   ::   g0, g1          ! coefficients for fitting the line of the ITD 
     72      REAL(wp), DIMENSION(jpij,jpl)   ::   hL, hR          ! left and right boundary for the ITD for each thickness 
     73      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice  ! local increment of ice area and volume 
     74      REAL(wp), DIMENSION(jpij)       ::   zhb0, zhb1      ! category boundaries for thinnes categories 
     75      REAL(wp), DIMENSION(jpij,0:jpl) ::   zhbnew          ! new boundaries of ice categories 
     76      !!------------------------------------------------------------------ 
    8877 
    8978      IF( kt == nit000 .AND. lwp) THEN 
     
    9483 
    9584      !----------------------------------------------------------------------------------------------- 
    96       !  3) Identify grid cells with ice 
     85      !  1) Identify grid cells with ice 
    9786      !----------------------------------------------------------------------------------------------- 
    9887      nidx = 0 ; idxice(:) = 0 
     
    10190            IF ( at_i(ji,jj) > epsi10 ) THEN 
    10291               nidx         = nidx + 1 
    103                nind_i(nidx) = ji 
    104                nind_j(nidx) = jj 
    10592               idxice( nidx ) = (jj - 1) * jpi + ji 
    106                zremap_flag(ji,jj) = 1 
    107             ELSE 
    108                zremap_flag(ji,jj) = 0 
    10993            ENDIF 
    11094         END DO 
    11195      END DO 
    112  
     96       
    11397      !----------------------------------------------------------------------------------------------- 
    114       !  4) Compute new category boundaries: 1:jpl-1 
     98      !  2) Compute new category boundaries 
    11599      !----------------------------------------------------------------------------------------------- 
    116       zdhice(:,:,:) = 0._wp 
    117       zhbnew(:,:,:) = 0._wp 
    118       zhb0(:,:) = hi_max(0) 
    119       zhb1(:,:) = hi_max(1) 
    120  
    121100      IF( nidx > 0 ) THEN 
    122101 
    123          ! Compute thickness change in each ice category 
     102         zdhice(:,:) = 0._wp 
     103         zhbnew(:,:) = 0._wp 
     104 
     105         CALL tab_3d_2d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i   ) 
     106         CALL tab_3d_2d( nidx, idxice(1:nidx), ht_ib_2d(1:nidx,1:jpl), ht_i_b ) 
     107         CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i    ) 
     108         CALL tab_3d_2d( nidx, idxice(1:nidx), a_ib_2d (1:nidx,1:jpl), a_i_b  ) 
     109 
    124110         DO jl = 1, jpl 
     111            ! Compute thickness change in each ice category 
    125112            DO ji = 1, nidx 
    126                ii = nind_i(ji) 
    127                ij = nind_j(ji) 
    128                zdhice(ii,ij,jl) = ht_i(ii,ij,jl) - ht_i_b(ii,ij,jl) 
     113               zdhice(ji,jl) = ht_i_2d(ji,jl) - ht_ib_2d(ji,jl) 
    129114            END DO 
    130115         END DO 
    131116          
     117         ! --- New boundaries for category 1:jpl-1 --- ! 
    132118         DO jl = 1, jpl - 1 
     119 
    133120            DO ji = 1, nidx 
    134                ii = nind_i(ji) 
    135                ij = nind_j(ji) 
    136121               ! 
    137                ! --- New boundary categories Hn* = Hn + Fn*dt --- ! 
     122               ! --- New boundary: Hn* = Hn + Fn*dt --- ! 
    138123               !     Fn*dt = ( fn + (fn+1 - fn)/(hn+1 - hn) * (Hn - hn) ) * dt = zdhice + zslope * (Hmax - ht_i_b) 
    139124               ! 
    140                IF    ( a_i_b(ii,ij,jl) >  epsi10 .AND. a_i_b(ii,ij,jl+1) >  epsi10 ) THEN   ! a_i_b(jl+1) & a_i_b(jl) /= 0 
    141                   zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( ht_i_b(ii,ij,jl+1) - ht_i_b(ii,ij,jl) ) 
    142                   zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - ht_i_b(ii,ij,jl) ) 
    143                ELSEIF( a_i_b(ii,ij,jl) >  epsi10 .AND. a_i_b(ii,ij,jl+1) <= epsi10 ) THEN   ! a_i_b(jl+1)=0 => Hn* = Hn + fn*dt 
    144                   zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 
    145                ELSEIF( a_i_b(ii,ij,jl) <= epsi10 .AND. a_i_b(ii,ij,jl+1) >  epsi10 ) THEN   ! a_i_b(jl)=0 => Hn* = Hn + fn+1*dt 
    146                   zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 
    147                ELSE                                                                         ! a_i_b(jl+1) & a_i_b(jl) = 0  
    148                   zhbnew(ii,ij,jl) = hi_max(jl) 
     125               IF    ( a_ib_2d(ji,jl) >  epsi10 .AND. a_ib_2d(ji,jl+1) >  epsi10 ) THEN   ! a(jl+1) & a(jl) /= 0 
     126                  zslope           = ( zdhice(ji,jl+1) - zdhice(ji,jl) ) / ( ht_ib_2d(ji,jl+1) - ht_ib_2d(ji,jl) ) 
     127                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) + zslope * ( hi_max(jl) - ht_ib_2d(ji,jl) ) 
     128               ELSEIF( a_ib_2d(ji,jl) >  epsi10 .AND. a_ib_2d(ji,jl+1) <= epsi10 ) THEN   ! a(jl+1)=0 => Hn* = Hn + fn*dt 
     129                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) 
     130               ELSEIF( a_ib_2d(ji,jl) <= epsi10 .AND. a_ib_2d(ji,jl+1) >  epsi10 ) THEN   ! a(jl)=0 => Hn* = Hn + fn+1*dt 
     131                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl+1) 
     132               ELSE                                                                       ! a(jl+1) & a(jl) = 0  
     133                  zhbnew(ji,jl) = hi_max(jl) 
    149134               ENDIF 
    150135             
     
    153138               !    Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible  
    154139               !          in lim_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
    155                IF( a_i(ii,ij,jl  ) > epsi10 .AND. ht_i(ii,ij,jl  ) > ( zhbnew(ii,ij,jl) - epsi10 ) )   zremap_flag(ii,ij) = 0 
    156                IF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) )   zremap_flag(ii,ij) = 0 
     140               IF( a_i_2d(ji,jl  ) > epsi10 .AND. ht_i_2d(ji,jl  ) > ( zhbnew(ji,jl) - epsi10 ) )   idxice(ji) = 0 
     141               IF( a_i_2d(ji,jl+1) > epsi10 .AND. ht_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) )   idxice(ji) = 0 
    157142                
    158143               ! 2) Hn-1 < Hn* < Hn+1   
    159                IF( zhbnew(ii,ij,jl) < hi_max(jl-1) )   zremap_flag(ii,ij) = 0 
    160                IF( zhbnew(ii,ij,jl) > hi_max(jl+1) )   zremap_flag(ii,ij) = 0 
     144               IF( zhbnew(ji,jl) < hi_max(jl-1) )   idxice(ji) = 0 
     145               IF( zhbnew(ji,jl) > hi_max(jl+1) )   idxice(ji) = 0 
    161146                
    162147            END DO 
    163148         END DO 
    164149          
    165          ! --- New boundary for category jpl --- ! 
    166          DO ji = 1, nidx 
    167             ii = nind_i(ji) 
    168             ij = nind_j(ji) 
    169             IF( a_i(ii,ij,jpl) > epsi10 ) THEN 
    170                zhbnew(ii,ij,jpl) = MAX( hi_max(jpl-1), 3._wp * ht_i(ii,ij,jpl) - 2._wp * zhbnew(ii,ij,jpl-1) ) 
     150         ! --- New boundaries for category jpl --- ! 
     151         DO ji = 1, nidx 
     152            IF( a_i_2d(ji,jpl) > epsi10 ) THEN 
     153               zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * ht_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) ) 
    171154            ELSE 
    172                zhbnew(ii,ij,jpl) = hi_max(jpl)   
    173             ENDIF 
    174              
    175             ! 1 condition for remapping for the 1st category 
     155               zhbnew(ji,jpl) = hi_max(jpl)   
     156            ENDIF 
     157             
     158            ! --- 1 additional condition for remapping (1st category) --- ! 
    176159            ! H0+epsi < h1(t) < H1-epsi  
    177160            !    h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible  
    178161            !    in lim_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
    179             IF( ht_i_b(ii,ij,1) < ( hi_max(0) + epsi10 ) )   zremap_flag(ii,ij) = 0 
    180             IF( ht_i_b(ii,ij,1) > ( hi_max(1) - epsi10 ) )   zremap_flag(ii,ij) = 0 
    181          END DO 
     162            IF( ht_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) )   idxice(ji) = 0 
     163            IF( ht_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) )   idxice(ji) = 0 
     164         END DO 
     165 
     166         !----------------------------------------------------------------------------------------------- 
     167         !  3) Identify cells where remapping 
     168         !----------------------------------------------------------------------------------------------- 
     169         nidx2 = 0 ; idxice2(:) = 0 
     170         DO ji = 1, nidx 
     171            IF( idxice(ji) /= 0 ) THEN 
     172               nidx2 = nidx2 + 1 
     173               idxice2(nidx2) = idxice(ji) 
     174               zhbnew(nidx2,:) = zhbnew(ji,:)  ! adjust zhbnew to new indices 
     175            ENDIF 
     176         END DO 
     177         idxice(:) = idxice2(:) 
     178         nidx      = nidx2 
    182179          
    183180      ENDIF 
    184181    
     182 
    185183      !----------------------------------------------------------------------------------------------- 
    186       !  5) Identify cells where ITD is to be remapped 
    187       !----------------------------------------------------------------------------------------------- 
    188       nidx = 0 
    189       DO jj = 1, jpj 
    190          DO ji = 1, jpi 
    191             IF( zremap_flag(ji,jj) == 1 ) THEN 
    192                nidx         = nidx + 1 
    193                nind_i(nidx) = ji 
    194                nind_j(nidx) = jj 
    195             ENDIF 
    196          END DO  
    197       END DO  
    198  
    199       !----------------------------------------------------------------------------------------------- 
    200       !  7) Compute g(h)  
     184      !  4) Compute g(h)  
    201185      !----------------------------------------------------------------------------------------------- 
    202186      IF( nidx > 0 ) THEN 
    203187 
    204          !- 7.1 g(h) for category 1 at start of time step 
    205          CALL lim_itd_glinear( 1, zhb0, zhb1, ht_i_b(:,:,1), g0(:,:,1), g1(:,:,1), hL(:,:,1),   & 
    206             &                  hR(:,:,1), zremap_flag ) 
    207           
    208          !- 7.2 Area lost due to melting of thin ice (first category,  1) 
    209          DO ji = 1, nidx 
    210             ii = nind_i(ji)  
    211             ij = nind_j(ji)  
    212              
    213             IF( a_i(ii,ij,1) > epsi10 ) THEN 
    214                 
    215                zdh0 = zdhice(ii,ij,1) !decrease of ice thickness in the lower category 
    216                 
    217                IF( zdh0 < 0.0 ) THEN      !remove area from category 1 
    218                   zdh0 = MIN( -zdh0, hi_max(1) ) 
    219                   !Integrate g(1) from 0 to dh0 to estimate area melted 
    220                   zetamax = MIN( zdh0, hR(ii,ij,1) ) - hL(ii,ij,1) 
     188         zhb0(:) = hi_max(0) ; zhb1(:) = hi_max(1) 
     189         g0(:,:) = 0._wp     ; g1(:,:) = 0._wp  
     190         hL(:,:) = 0._wp     ; hR(:,:) = 0._wp  
     191          
     192         DO jl = 1, jpl 
     193 
     194            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_ib_1d(1:nidx), ht_i_b(:,:,jl) ) 
     195            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   ) 
     196            CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    ) 
     197            CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    ) 
     198 
     199            IF( jl == 1 ) THEN 
     200                
     201               ! --- g(h) for category 1 --- ! 
     202               CALL lim_itd_glinear( zhb0(1:nidx), zhb1(1:nidx), ht_ib_1d(1:nidx), a_i_1d(1:nidx),  &   ! in 
     203                  &                  g0(1:nidx,1), g1(1:nidx,1), hL(1:nidx,1)    , hR(1:nidx,1) )       ! out 
     204                
     205               ! Area lost due to melting of thin ice 
     206               DO ji = 1, nidx 
    221207                   
    222                   IF( zetamax > 0.0 ) THEN 
    223                      zx1    = zetamax 
    224                      zx2    = 0.5 * zetamax * zetamax  
    225                      zda0   = g1(ii,ij,1) * zx2 + g0(ii,ij,1) * zx1                        ! ice area removed 
    226                      zdamax = a_i(ii,ij,1) * (1.0 - ht_i(ii,ij,1) / ht_i_b(ii,ij,1) ) ! Constrain new thickness <= ht_i                 
    227                      zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting  
    228                      !     of thin ice (zdamax > 0) 
    229                      ! Remove area, conserving volume 
    230                      ht_i(ii,ij,1) = ht_i(ii,ij,1) * a_i(ii,ij,1) / ( a_i(ii,ij,1) - zda0 ) 
    231                      a_i(ii,ij,1)  = a_i(ii,ij,1) - zda0 
    232                      v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) ! clem-useless ? 
     208                  IF( a_i_1d(ji) > epsi10 ) THEN 
     209                      
     210                     zdh0 =  ht_i_1d(ji) - ht_ib_1d(ji)                 
     211                     IF( zdh0 < 0.0 ) THEN      !remove area from category 1 
     212                        zdh0 = MIN( -zdh0, hi_max(1) ) 
     213                        !Integrate g(1) from 0 to dh0 to estimate area melted 
     214                        zetamax = MIN( zdh0, hR(ji,1) ) - hL(ji,1) 
     215                         
     216                        IF( zetamax > 0.0 ) THEN 
     217                           zx1    = zetamax 
     218                           zx2    = 0.5 * zetamax * zetamax  
     219                           zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                        ! ice area removed 
     220                           zdamax = a_i_1d(ji) * (1.0 - ht_i_1d(ji) / ht_ib_1d(ji) ) ! Constrain new thickness <= ht_i                 
     221                           zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting  
     222                           !     of thin ice (zdamax > 0) 
     223                           ! Remove area, conserving volume 
     224                           ht_i_1d(ji) = ht_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 ) 
     225                           a_i_1d(ji)  = a_i_1d(ji) - zda0 
     226                           v_i_1d(ji)  = a_i_1d(ji) * ht_i_1d(ji) ! clem-useless ? 
     227                        ENDIF 
     228                         
     229                     ELSE ! if ice accretion zdh0 > 0 
     230                        ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 
     231                        zhbnew(ji,0) = MIN( zdh0, hi_max(1) )  
     232                     ENDIF 
     233                      
    233234                  ENDIF 
    234235                   
    235                ELSE ! if ice accretion zdh0 > 0 
    236                   ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 
    237                   zhbnew(ii,ij,0) = MIN( zdh0, hi_max(1) )  
    238                ENDIF 
    239                 
    240             ENDIF 
    241              
    242          END DO 
    243           
    244          !- 7.3 g(h) for each thickness category   
    245          DO jl = 1, jpl 
    246             CALL lim_itd_glinear( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),  & 
    247                &                  g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) 
     236               END DO 
     237 
     238               CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   ) 
     239               CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    ) 
     240               CALL tab_1d_2d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    ) 
     241                
     242            ENDIF ! jl=1 
     243             
     244            ! --- g(h) for each thickness category --- !   
     245            CALL lim_itd_glinear( zhbnew(1:nidx,jl-1), zhbnew(1:nidx,jl), ht_i_1d(1:nidx), a_i_1d(1:nidx),  &  ! in 
     246               &                  g0(1:nidx,jl)      , g1(1:nidx,jl)    , hL(1:nidx,jl)  , hR(1:nidx,jl) )     ! out 
     247 
    248248         END DO 
    249249          
    250250         !----------------------------------------------------------------------------------------------- 
    251          !  8) Compute area and volume to be shifted across each boundary (Eq. 18) 
     251         !  5) Compute area and volume to be shifted across each boundary (Eq. 18) 
    252252         !----------------------------------------------------------------------------------------------- 
    253           
    254253         DO jl = 1, jpl - 1 
    255             DO jj = 1, jpj 
    256                DO ji = 1, jpi 
    257                   zdonor(ji,jj,jl) = 0 
    258                   zdaice(ji,jj,jl) = 0.0 
    259                   zdvice(ji,jj,jl) = 0.0 
    260                END DO 
    261             END DO 
    262254             
    263255            DO ji = 1, nidx 
    264                ii = nind_i(ji) 
    265                ij = nind_j(ji) 
    266256                
    267257               ! left and right integration limits in eta space 
    268                IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1 
    269                   zetamin = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl)       ! hi_max(jl) - hL  
    270                   zetamax = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) ! hR - hL 
    271                   zdonor(ii,ij,jl) = jl 
    272                ELSE                                    ! Hn* <= Hn => transfer from jl+1 to jl 
     258               IF (zhbnew(ji,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1 
     259                  zetamin = MAX( hi_max(jl)   , hL(ji,jl) ) - hL(ji,jl)   ! hi_max(jl) - hL  
     260                  zetamax = MIN( zhbnew(ji,jl), hR(ji,jl) ) - hL(ji,jl)  ! hR - hL 
     261                  jdonor(ji,jl) = jl 
     262               ELSE                                 ! Hn* <= Hn => transfer from jl+1 to jl 
    273263                  zetamin = 0.0 
    274                   zetamax = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1)   ! hi_max(jl) - hL 
    275                   zdonor(ii,ij,jl) = jl + 1 
     264                  zetamax = MIN( hi_max(jl), hR(ji,jl+1) ) - hL(ji,jl+1)  ! hi_max(jl) - hL 
     265                  jdonor(ji,jl) = jl + 1 
    276266               ENDIF 
    277267               zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin 
     
    284274               zwk2 = zwk2 * zetamax 
    285275               zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 ) 
    286                nd   = zdonor(ii,ij,jl) 
    287                zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 
    288                zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 
     276               jcat = jdonor(ji,jl) 
     277               zdaice(ji,jl) = g1(ji,jcat)*zx2 + g0(ji,jcat)*zx1 
     278               zdvice(ji,jl) = g1(ji,jcat)*zx3 + g0(ji,jcat)*zx2 + zdaice(ji,jl)*hL(ji,jcat) 
    289279                
    290280            END DO 
    291281         END DO 
    292282          
    293          !!---------------------------------------------------------------------------------------------- 
    294          !! 9) Shift ice between categories 
    295          !!---------------------------------------------------------------------------------------------- 
    296          CALL lim_itd_shiftice ( zdonor, zdaice, zdvice ) 
    297           
    298          !!---------------------------------------------------------------------------------------------- 
    299          !! 10) Make sure ht_i >= minimum ice thickness hi_min 
    300          !!---------------------------------------------------------------------------------------------- 
    301           
    302          DO ji = 1, nidx 
    303             ii = nind_i(ji) 
    304             ij = nind_j(ji) 
    305             IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN 
    306                a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin  
     283         !---------------------------------------------------------------------------------------------- 
     284         ! 6) Shift ice between categories 
     285         !---------------------------------------------------------------------------------------------- 
     286         CALL lim_itd_shiftice ( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) ) 
     287          
     288         !---------------------------------------------------------------------------------------------- 
     289         ! 7) Make sure ht_i >= minimum ice thickness hi_min 
     290         !---------------------------------------------------------------------------------------------- 
     291         CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1)   ) 
     292         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    ) 
     293         CALL tab_2d_1d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)    ) 
     294          
     295         DO ji = 1, nidx 
     296            IF ( a_i_1d(ji) > epsi10 .AND. ht_i_1d(ji) < rn_himin ) THEN 
     297               a_i_1d (ji) = a_i_1d(ji) * ht_i_1d(ji) / rn_himin  
    307298               ! MV MP 2016 
    308299               IF ( nn_pnd_scheme > 0 ) THEN 
    309                   a_ip(ii,ij,1) = a_ip(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 
     300                  a_ip_1d(ji) = a_ip_1d(ji) * ht_i_1d(ji) / rn_himin 
    310301               ENDIF 
    311302               ! END MV MP 2016 
    312                ht_i(ii,ij,1) = rn_himin 
    313             ENDIF 
    314          END DO 
    315           
     303               ht_i_1d(ji) = rn_himin 
     304            ENDIF 
     305         END DO 
     306 
     307         CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1)   ) 
     308         CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    ) 
     309         CALL tab_1d_2d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)    ) 
     310         
    316311      ENDIF 
    317312 
    318       CALL wrk_dealloc( jpi,jpj, zremap_flag ) 
    319       CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) 
    320       CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR ) 
    321       CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    322       CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
    323       CALL wrk_dealloc( jpi,jpj, zhb0, zhb1 ) 
    324  
    325313   END SUBROUTINE lim_itd_th_rem 
    326314 
    327315 
    328    SUBROUTINE lim_itd_glinear( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 
     316   SUBROUTINE lim_itd_glinear( HbL, Hbr, phice, paice, pg0, pg1, phL, phR ) 
    329317      !!------------------------------------------------------------------ 
    330318      !!                ***  ROUTINE lim_itd_glinear *** 
     
    336324      !!                 
    337325      !!------------------------------------------------------------------ 
    338       INTEGER                     , INTENT(in   ) ::   num_cat      ! category index 
    339       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   HbL, HbR     ! left and right category boundaries 
    340       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   hice         ! ice thickness 
    341       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   g0, g1       ! coefficients in linear equation for g(eta) 
    342       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   hL           ! min value of range over which g(h) > 0 
    343       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   hR           ! max value of range over which g(h) > 0 
    344       INTEGER , DIMENSION(jpi,jpj), INTENT(in   ) ::   zremap_flag  ! 
     326      REAL(wp), DIMENSION(:), INTENT(in   ) ::   HbL, HbR      ! left and right category boundaries 
     327      REAL(wp), DIMENSION(:), INTENT(in   ) ::   phice, paice  ! ice thickness and concentration 
     328      REAL(wp), DIMENSION(:), INTENT(inout) ::   pg0, pg1      ! coefficients in linear equation for g(eta) 
     329      REAL(wp), DIMENSION(:), INTENT(inout) ::   phL, phR      ! min and max value of range over which g(h) > 0 
    345330      ! 
    346       INTEGER  ::   ji,jj        ! horizontal indices 
     331      INTEGER  ::   ji           ! horizontal indices 
    347332      REAL(wp) ::   zh13         ! HbL + 1/3 * (HbR - HbL) 
    348333      REAL(wp) ::   zh23         ! HbL + 2/3 * (HbR - HbL) 
     
    351336      !!------------------------------------------------------------------ 
    352337      ! 
    353       DO jj = 1, jpj 
    354          DO ji = 1, jpi 
     338         DO ji = 1, nidx 
    355339            ! 
    356             IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10   & 
    357                &                        .AND. hice(ji,jj)        > 0._wp )  THEN 
     340            IF( paice(ji) > epsi10  .AND. phice(ji) > 0._wp )  THEN 
    358341 
    359342               ! Initialize hL and hR 
    360                hL(ji,jj) = HbL(ji,jj) 
    361                hR(ji,jj) = HbR(ji,jj) 
     343               phL(ji) = HbL(ji) 
     344               phR(ji) = HbR(ji) 
    362345 
    363346               ! Change hL or hR if hice falls outside central third of range, 
    364347               ! so that hice is in the central third of the range [HL HR] 
    365                zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) +       hR(ji,jj) ) 
    366                zh23 = 1.0 / 3.0 * (       hL(ji,jj) + 2.0 * hR(ji,jj) ) 
    367  
    368                IF    ( hice(ji,jj) < zh13 ) THEN   ;   hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) ! move HR to the left 
    369                ELSEIF( hice(ji,jj) > zh23 ) THEN   ;   hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj) ! move HL to the right 
     348               zh13 = 1.0 / 3.0 * ( 2.0 * phL(ji) +       phR(ji) ) 
     349               zh23 = 1.0 / 3.0 * (       phL(ji) + 2.0 * phR(ji) ) 
     350 
     351               IF    ( phice(ji) < zh13 ) THEN   ;   phR(ji) = 3._wp * phice(ji) - 2._wp * phL(ji) ! move HR to the left 
     352               ELSEIF( phice(ji) > zh23 ) THEN   ;   phL(ji) = 3._wp * phice(ji) - 2._wp * phR(ji) ! move HL to the right 
    370353               ENDIF 
    371354 
    372355               ! Compute coefficients of g(eta) = g0 + g1*eta 
    373                zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj)) 
    374                zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 
    375                zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr 
    376                g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 )      ! Eq. 14 
    377                g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) ! Eq. 14 
     356               zdhr = 1._wp / (phR(ji) - phL(ji)) 
     357               zwk1 = 6._wp * paice(ji) * zdhr 
     358               zwk2 = ( phice(ji) - phL(ji) ) * zdhr 
     359               pg0(ji) = zwk1 * ( 2._wp / 3._wp - zwk2 )      ! Eq. 14 
     360               pg1(ji) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) ! Eq. 14 
    378361               ! 
    379362            ELSE  ! remap_flag = .false. or a_i < epsi10  
    380                hL(ji,jj) = 0._wp 
    381                hR(ji,jj) = 0._wp 
    382                g0(ji,jj) = 0._wp 
    383                g1(ji,jj) = 0._wp 
     363               phL(ji) = 0._wp 
     364               phR(ji) = 0._wp 
     365               pg0(ji) = 0._wp 
     366               pg1(ji) = 0._wp 
    384367            ENDIF 
    385368            ! 
    386369         END DO 
    387       END DO 
    388370      ! 
    389371   END SUBROUTINE lim_itd_glinear 
    390372 
    391373 
    392    SUBROUTINE lim_itd_shiftice( zdonor, zdaice, zdvice ) 
     374   SUBROUTINE lim_itd_shiftice( kdonor, pdaice, pdvice ) 
    393375      !!------------------------------------------------------------------ 
    394376      !!                ***  ROUTINE lim_itd_shiftice *** 
     
    399381      !! ** Method  : 
    400382      !!------------------------------------------------------------------ 
    401       INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(in   ) ::   zdonor   ! donor category index 
    402       REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) ::   zdaice   ! ice area transferred across boundary 
    403       REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) ::   zdvice   ! ice volume transferred across boundary 
    404  
    405       INTEGER ::   ji, jj, jl, jl2, jl1, jk   ! dummy loop indices 
    406  
    407       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zaTsfn 
    408       REAL(wp), DIMENSION(jpi,jpj)     ::   zworka            ! temporary array used here 
     383      INTEGER , DIMENSION(:,:), INTENT(in) ::   kdonor   ! donor category index 
     384      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdaice   ! ice area transferred across boundary 
     385      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdvice   ! ice volume transferred across boundary 
     386 
     387      INTEGER ::   ii,ij, ji, jj, jl, jl2, jl1, jk   ! dummy loop indices 
     388 
     389      REAL(wp), DIMENSION(jpij,jpl) ::   zaTsfn 
     390      REAL(wp), DIMENSION(jpij)     ::   zworka            ! temporary array used here 
     391      REAL(wp), DIMENSION(jpij)     ::   zworkv            ! temporary array used here 
    409392 
    410393      REAL(wp) ::   ztrans             ! ice/snow transferred 
    411394      !!------------------------------------------------------------------ 
     395          
     396      CALL tab_3d_2d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i   ) 
     397      CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i    ) 
     398      CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i    ) 
     399      CALL tab_3d_2d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s    ) 
     400      CALL tab_3d_2d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i   ) 
     401      CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i  ) 
     402      CALL tab_3d_2d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip   ) 
     403      CALL tab_3d_2d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip   ) 
     404      CALL tab_3d_2d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su   ) 
    412405 
    413406      !---------------------------------------------------------------------------------------------- 
     
    415408      !---------------------------------------------------------------------------------------------- 
    416409      DO jl = 1, jpl 
    417          DO jj = 1, jpj 
    418             DO ji = 1, jpi 
    419                zaTsfn(ji,jj,jl) = a_i(ji,jj,jl) * t_su(ji,jj,jl) 
    420             END DO 
     410         DO ji = 1, nidx 
     411            zaTsfn(ji,jl) = a_i_2d(ji,jl) * t_su_2d(ji,jl) 
    421412         END DO 
    422413      END DO 
     
    426417      !------------------------------------------------------------------------------- 
    427418      DO jl = 1, jpl - 1 
    428          DO jj = 1, jpj 
    429             DO ji = 1, jpi 
    430  
    431                jl1 = zdonor(ji,jj,jl) 
    432                rswitch       = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl1) - epsi10 ) ) 
    433                zworka(ji,jj) = zdvice(ji,jj,jl) / MAX( v_i(ji,jj,jl1), epsi10 ) * rswitch 
    434  
    435                IF( jl1 == jl) THEN   ;   jl2 = jl1+1 
    436                ELSE                  ;   jl2 = jl  
    437                ENDIF 
    438                 
    439                ! Ice areas 
    440                a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - zdaice(ji,jj,jl) 
    441                a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + zdaice(ji,jj,jl) 
    442                 
    443                ! Ice volumes 
    444                v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - zdvice(ji,jj,jl)  
    445                v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + zdvice(ji,jj,jl) 
    446                 
    447                ! Snow volumes 
    448                ztrans         = v_s(ji,jj,jl1) * zworka(ji,jj) 
    449                v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - ztrans 
    450                v_s(ji,jj,jl2) = v_s(ji,jj,jl2) + ztrans  
    451                 
    452                ! Snow heat content   
    453                ztrans             = e_s(ji,jj,1,jl1) * zworka(ji,jj) 
    454                e_s(ji,jj,1,jl1)   = e_s(ji,jj,1,jl1) - ztrans 
    455                e_s(ji,jj,1,jl2)   = e_s(ji,jj,1,jl2) + ztrans 
    456                 
    457                ! Ice age  
    458                ztrans             = oa_i(ji,jj,jl1) * zdaice(ji,jj,jl) 
    459                oa_i(ji,jj,jl1)    = oa_i(ji,jj,jl1) - ztrans 
    460                oa_i(ji,jj,jl2)    = oa_i(ji,jj,jl2) + ztrans 
    461                 
    462                ! Ice salinity 
    463                ztrans             = smv_i(ji,jj,jl1) * zworka(ji,jj) 
    464                smv_i(ji,jj,jl1)   = smv_i(ji,jj,jl1) - ztrans 
    465                smv_i(ji,jj,jl2)   = smv_i(ji,jj,jl2) + ztrans 
    466                 
    467                ! Surface temperature 
    468                ztrans             = t_su(ji,jj,jl1) * zdaice(ji,jj,jl) 
    469                zaTsfn(ji,jj,jl1)  = zaTsfn(ji,jj,jl1) - ztrans 
    470                zaTsfn(ji,jj,jl2)  = zaTsfn(ji,jj,jl2) + ztrans 
    471                 
    472                ! MV MP 2016  
    473                IF ( nn_pnd_scheme > 0 ) THEN 
    474                   ! Pond fraction 
    475                   ztrans          = a_ip(ji,jj,jl1) * zdaice(ji,jj,jl) 
    476                   a_ip(ji,jj,jl1) = a_ip(ji,jj,jl1) - ztrans 
    477                   a_ip(ji,jj,jl2) = a_ip(ji,jj,jl2) + ztrans 
    478                    
    479                   ! Pond volume 
    480                   ztrans          = v_ip(ji,jj,jl1) * zdaice(ji,jj,jl) 
    481                   v_ip(ji,jj,jl1) = v_ip(ji,jj,jl1) - ztrans 
    482                   v_ip(ji,jj,jl2) = v_ip(ji,jj,jl2) + ztrans 
    483                ENDIF 
    484                ! END MV MP 2016  
    485                 
     419         DO ji = 1, nidx 
     420             
     421            jl1 = kdonor(ji,jl) 
     422            IF    ( jl1 == jl  ) THEN   ;   jl2 = jl1+1 
     423            ELSE                        ;   jl2 = jl  
     424            ENDIF 
     425 
     426            rswitch    = MAX( 0._wp , SIGN( 1._wp , v_i_2d(ji,jl1) - epsi10 ) ) 
     427            zworkv(ji) = pdvice(ji,jl) / MAX( v_i_2d(ji,jl1), epsi10 ) * rswitch 
     428 
     429            rswitch    = MAX( 0._wp , SIGN( 1._wp , a_i_2d(ji,jl1) - epsi10 ) ) 
     430            zworka(ji) = pdaice(ji,jl) / MAX( a_i_2d(ji,jl1), epsi10 ) * rswitch         
     431                
     432            ! Ice areas 
     433            a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl) 
     434            a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl) 
     435                
     436            ! Ice volumes 
     437            v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl)  
     438            v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl) 
     439             
     440            ! Snow volumes 
     441            ztrans         = v_s_2d(ji,jl1) * zworkv(ji) 
     442            v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans 
     443            v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans  
     444                        
     445            ! Ice age  
     446            ztrans             = oa_i_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work 
     447            oa_i_2d(ji,jl1)    = oa_i_2d(ji,jl1) - ztrans 
     448            oa_i_2d(ji,jl2)    = oa_i_2d(ji,jl2) + ztrans 
     449             
     450            ! Ice salinity 
     451            ztrans             = smv_i_2d(ji,jl1) * zworkv(ji) 
     452            smv_i_2d(ji,jl1)   = smv_i_2d(ji,jl1) - ztrans 
     453            smv_i_2d(ji,jl2)   = smv_i_2d(ji,jl2) + ztrans 
     454             
     455            ! Surface temperature 
     456            ztrans          = t_su_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work 
     457            zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans 
     458            zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans 
     459                
     460            ! MV MP 2016  
     461            IF ( nn_pnd_scheme > 0 ) THEN 
     462               ! Pond fraction 
     463               ztrans          = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work 
     464               a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans 
     465               a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans 
     466                
     467               ! Pond volume (also proportional to da/a) 
     468               ztrans          = v_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work 
     469               v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 
     470               v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 
     471            ENDIF 
     472            ! END MV MP 2016  
     473                
     474         END DO 
     475 
     476         ! Snow heat content 
     477         DO jk = 1, nlay_s 
     478 
     479            DO ji = 1, nidx 
     480               ii = MOD( idxice(ji) - 1, jpi ) + 1 
     481               ij = ( idxice(ji) - 1 ) / jpi + 1 
     482 
     483               jl1 = kdonor(ji,jl) 
     484               IF(jl1 == jl) THEN  ;  jl2 = jl+1 
     485               ELSE                ;  jl2 = jl 
     486               ENDIF 
     487                
     488               ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji) 
     489               e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans 
     490               e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans  
    486491            END DO 
    487492         END DO 
     493 
    488494          
    489495         ! Ice heat content 
    490496         DO jk = 1, nlay_i 
    491             DO jj = 1, jpj 
    492                DO ji = 1, jpi 
    493                    
    494                   jl1 = zdonor(ji,jj,jl) 
    495                   IF(jl1 == jl) THEN  ;  jl2 = jl+1 
    496                   ELSE                ;  jl2 = jl 
    497                   ENDIF 
    498                    
    499                   ztrans            = e_i(ji,jj,jk,jl1) * zworka(ji,jj) 
    500                   e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - ztrans 
    501                   e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + ztrans  
    502                END DO 
     497            DO ji = 1, nidx 
     498               ii = MOD( idxice(ji) - 1, jpi ) + 1 
     499               ij = ( idxice(ji) - 1 ) / jpi + 1 
     500                
     501               jl1 = kdonor(ji,jl) 
     502               IF(jl1 == jl) THEN  ;  jl2 = jl+1 
     503               ELSE                ;  jl2 = jl 
     504               ENDIF 
     505                
     506               ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji) 
     507               e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans 
     508               e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans  
    503509            END DO 
    504510         END DO 
     
    506512      END DO                   ! boundaries, 1 to jpl-1 
    507513       
    508       ! Update ice thickness and temperature 
     514      !------------------------------------------------------------------------------- 
     515      ! 3) Update ice thickness and temperature 
     516      !------------------------------------------------------------------------------- 
    509517      DO jl = 1, jpl 
    510          DO jj = 1, jpj 
    511             DO ji = 1, jpi 
    512                IF ( a_i(ji,jj,jl) > epsi10 ) THEN  
    513                   ht_i(ji,jj,jl)  =  v_i   (ji,jj,jl) / a_i(ji,jj,jl)  
    514                   t_su(ji,jj,jl)  =  zaTsfn(ji,jj,jl) / a_i(ji,jj,jl)  
    515                ELSE 
    516                   ht_i(ji,jj,jl)  = 0._wp 
    517                   t_su(ji,jj,jl)  = rt0 
    518                ENDIF 
    519             END DO 
    520          END DO 
    521       END DO       
     518         DO ji = 1, nidx 
     519            IF ( a_i_2d(ji,jl) > epsi10 ) THEN  
     520               ht_i_2d(ji,jl)  =  v_i_2d   (ji,jl) / a_i_2d(ji,jl)  
     521               t_su_2d(ji,jl)  =  zaTsfn(ji,jl) / a_i_2d(ji,jl)  
     522            ELSE 
     523               ht_i_2d(ji,jl)  = 0._wp 
     524               t_su_2d(ji,jl)  = rt0 
     525            ENDIF 
     526         END DO 
     527      END DO 
     528       
     529      CALL tab_2d_3d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i   ) 
     530      CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i    ) 
     531      CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i    ) 
     532      CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s    ) 
     533      CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i   ) 
     534      CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i  ) 
     535      CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip   ) 
     536      CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip   ) 
     537      CALL tab_2d_3d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su   ) 
     538 
    522539      ! 
    523540   END SUBROUTINE lim_itd_shiftice 
     
    530547      !! ** Purpose : rebin - rebins thicknesses into defined categories 
    531548      !! 
    532       !! ** Method  : 
     549      !! ** Method  : If a category thickness is out of bounds, shift part (for down to top) 
     550      !!              or entire (for top to down) area, volume, and energy 
     551      !!              to the neighboring category 
    533552      !!------------------------------------------------------------------ 
    534553      INTEGER ::   ji, jj, jl   ! dummy loop indices 
    535       INTEGER ::   zshiftflag          ! = .true. if ice must be shifted 
    536  
    537       INTEGER , DIMENSION(jpi,jpj,jpl) ::   zdonor           ! donor category index 
    538       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zdaice, zdvice   ! ice area and volume transferred 
     554      ! 
     555      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor           ! donor category index 
     556      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred 
    539557      !!------------------------------------------------------------------ 
    540558       
    541       !------------------------------------------------------------------------------ 
    542       ! 1) Compute ice thickness. 
    543       !------------------------------------------------------------------------------ 
    544       DO jl = 1, jpl 
     559      DO jl = 1, jpl-1 
     560         jdonor(:,jl) = 0 
     561         zdaice(:,jl) = 0._wp 
     562         zdvice(:,jl) = 0._wp 
     563      END DO 
     564 
     565      !--------------------------------------- 
     566      ! identify thicknesses that are too big 
     567      !--------------------------------------- 
     568      DO jl = 1, jpl - 1 
     569 
     570         nidx = 0 ; idxice(:) = 0 
    545571         DO jj = 1, jpj 
    546             DO ji = 1, jpi  
    547                rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
    548                ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 
    549             END DO 
    550          END DO 
    551       END DO 
    552  
    553       !------------------------------------------------------------------------------ 
    554       ! 2) If a category thickness is not in bounds, shift the 
    555       ! entire area, volume, and energy to the neighboring category 
    556       !------------------------------------------------------------------------------ 
    557       !------------------------- 
    558       ! Initialize shift arrays 
    559       !------------------------- 
    560       DO jl = 1, jpl 
    561          zdonor(:,:,jl) = 0 
    562          zdaice(:,:,jl) = 0._wp 
    563          zdvice(:,:,jl) = 0._wp 
    564       END DO 
    565  
    566       !------------------------- 
    567       ! Move thin categories up 
    568       !------------------------- 
    569  
    570       DO jl = 1, jpl - 1  ! loop over category boundaries 
    571  
    572          !--------------------------------------- 
    573          ! identify thicknesses that are too big 
    574          !--------------------------------------- 
    575          zshiftflag = 0 
    576  
    577          DO jj = 1, jpj  
    578             DO ji = 1, jpi  
    579                IF( a_i(ji,jj,jl) > epsi10 .AND. ht_i(ji,jj,jl) > hi_max(jl) ) THEN  
    580                   zshiftflag        = 1 
    581                   zdonor(ji,jj,jl)  = jl  
    582                   ! clem: how much of a_i you send in cat sup is somewhat arbitrary 
    583                   zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl)   
    584                   zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi20 ) 
    585                ENDIF 
    586             END DO 
    587          END DO 
    588          IF(lk_mpp)   CALL mpp_max( zshiftflag ) ! clem: for reproducibility ??? 
    589  
    590          IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
    591             CALL lim_itd_shiftice( zdonor, zdaice, zdvice ) 
     572            DO ji = 1, jpi 
     573               IF( a_i(ji,jj,jl) > epsi10 .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN 
     574                  nidx = nidx + 1 
     575                  idxice( nidx ) = (jj - 1) * jpi + ji                   
     576               ENDIF 
     577            ENDDO 
     578         ENDDO 
     579          
     580!!clem   CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   ) 
     581         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    ) 
     582         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    ) 
     583 
     584         DO ji = 1, nidx 
     585            jdonor(ji,jl)  = jl  
     586            ! how much of a_i you send in cat sup is somewhat arbitrary 
     587!!clem: these do not work properly after a restart (I do not know why) 
     588!!          zdaice(ji,jl)  = a_i_1d(ji) * ( ht_i_1d(ji) - hi_max(jl) + epsi10 ) / ht_i_1d(ji)   
     589!!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 
     590!!clem: these do not work properly after a restart (I do not know why) 
     591!!            zdaice(ji,jl)  = a_i_1d(ji) 
     592!!            zdvice(ji,jl)  = v_i_1d(ji) 
     593!!clem: these are from UCL and work ok 
     594            zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp 
     595            zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
     596             
     597         END DO 
     598 
     599         IF( nidx > 0 ) THEN 
     600            CALL lim_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl=>jl+1 
    592601            ! Reset shift parameters 
    593             zdonor(:,:,jl) = 0 
    594             zdaice(:,:,jl) = 0._wp 
    595             zdvice(:,:,jl) = 0._wp 
     602            jdonor(1:nidx,jl) = 0 
     603            zdaice(1:nidx,jl) = 0._wp 
     604            zdvice(1:nidx,jl) = 0._wp 
    596605         ENDIF 
    597606         ! 
    598607      END DO 
    599608 
    600       !---------------------------- 
    601       ! Move thick categories down 
    602       !---------------------------- 
    603  
    604       DO jl = jpl - 1, 1, -1       ! loop over category boundaries 
    605  
    606          !----------------------------------------- 
    607          ! Identify thicknesses that are too small 
    608          !----------------------------------------- 
    609          zshiftflag = 0 
    610  
     609      !----------------------------------------- 
     610      ! Identify thicknesses that are too small 
     611      !----------------------------------------- 
     612      DO jl = jpl - 1, 1, -1 
     613 
     614         nidx = 0 ; idxice(:) = 0 
    611615         DO jj = 1, jpj 
    612616            DO ji = 1, jpi 
    613                IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 
    614                   ! 
    615                   zshiftflag = 1 
    616                   zdonor(ji,jj,jl) = jl + 1 
    617                   zdaice(ji,jj,jl) = a_i(ji,jj,jl+1)  
    618                   zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 
    619                ENDIF 
    620             END DO 
    621          END DO 
    622  
    623          IF(lk_mpp)   CALL mpp_max( zshiftflag ) ! clem: for reproducibility ??? 
    624           
    625          IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
    626             CALL lim_itd_shiftice( zdonor, zdaice, zdvice ) 
     617               IF( a_i(ji,jj,jl+1) > epsi10 .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN 
     618                  nidx = nidx + 1 
     619                  idxice( nidx ) = (jj - 1) * jpi + ji                   
     620               ENDIF 
     621            ENDDO 
     622         ENDDO 
     623 
     624         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl+1)    ) ! jl+1 is ok 
     625         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl+1)    ) ! jl+1 is ok 
     626         DO ji = 1, nidx 
     627            jdonor(ji,jl) = jl + 1 
     628            zdaice(ji,jl) = a_i_1d(ji)  
     629            zdvice(ji,jl) = v_i_1d(ji) 
     630         END DO 
     631          
     632         IF( nidx > 0 ) THEN 
     633            CALL lim_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl+1=>jl 
    627634            ! Reset shift parameters 
    628             zdonor(:,:,jl) = 0 
    629             zdaice(:,:,jl) = 0._wp 
    630             zdvice(:,:,jl) = 0._wp 
     635            jdonor(1:nidx,jl) = 0 
     636            zdaice(1:nidx,jl) = 0._wp 
     637            zdvice(1:nidx,jl) = 0._wp 
    631638         ENDIF 
    632639 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90

    r8321 r8369  
    110110 
    111111      IF ( .NOT. ln_pnd ) THEN 
    112          WRITE(numout,*) 
    113          WRITE(numout,*) ' Melt ponds are not activated ' 
    114          WRITE(numout,*) ' ln_pnd_rad and ln_pnd_fw set to .FALSE. ' 
    115          WRITE(numout,*) ' nn_pnd_scheme, rn_apnd, rn_hpnd set to zero ' 
     112         IF(lwp) WRITE(numout,*) 
     113         IF(lwp) WRITE(numout,*) ' Melt ponds are not activated ' 
     114         IF(lwp) WRITE(numout,*) ' ln_pnd_rad and ln_pnd_fw set to .FALSE. ' 
     115         IF(lwp) WRITE(numout,*) ' nn_pnd_scheme, rn_apnd, rn_hpnd set to zero ' 
    116116         ln_pnd_rad    = .FALSE. 
    117117         ln_pnd_fw     = .FALSE. 
     
    130130 
    131131      IF ( ln_pnd .AND. ( nn_pnd_scheme == 0 ) .AND. ( ln_pnd_fw ) ) THEN 
    132          WRITE(numout,*) ' Prescribed melt ponds do not conserve fresh water mass, hence ln_pnd_fw must be set to false '  
     132         IF(lwp) WRITE(numout,*) ' Prescribed melt ponds do not conserve fresh water mass, hence ln_pnd_fw must be set to false '  
    133133         ln_pnd_fw     = .FALSE. 
    134134         IF(lwp) THEN                     ! control print 
     
    138138 
    139139      IF ( ln_pnd .AND. ( nn_pnd_scheme == 2 ) .AND. ( jpl == 1 ) ) THEN 
    140          WRITE(numout,*) ' Topographic melt ponds are incompatible with jpl = 1 ' 
    141          WRITE(numout,*) ' Run aborted ' 
     140         IF(lwp) WRITE(numout,*) ' Topographic melt ponds are incompatible with jpl = 1 ' 
     141         IF(lwp) WRITE(numout,*) ' Run aborted ' 
    142142         CALL ctl_stop( 'STOP', 'lim_mp_init: uncompatible options, reset namelist_ice_ref ' ) 
    143143      ENDIF 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limtab.F90

    r8342 r8369  
    1313   USE par_kind 
    1414   USE par_oce 
     15   USE ice, ONLY : jpl 
    1516    
    1617   IMPLICIT NONE 
    1718   PRIVATE 
    1819 
    19    PUBLIC   tab_2d_1d   ! called by limthd 
    20    PUBLIC   tab_1d_2d   ! called by limthd 
     20   PUBLIC   tab_3d_2d 
     21   PUBLIC   tab_2d_1d 
     22   PUBLIC   tab_2d_3d 
     23   PUBLIC   tab_1d_2d 
    2124 
    2225   !!---------------------------------------------------------------------- 
     
    2629   !!---------------------------------------------------------------------- 
    2730CONTAINS 
     31 
     32 
     33   SUBROUTINE tab_3d_2d( ndim1d, tab_ind, tab1d, tab2d ) 
     34      !!---------------------------------------------------------------------- 
     35      !!                  ***  ROUTINE tab_2d_1d  *** 
     36      !!---------------------------------------------------------------------- 
     37      INTEGER                         , INTENT(in   ) ::   ndim1d   ! 1d size 
     38      INTEGER , DIMENSION(ndim1d)     , INTENT(in   ) ::   tab_ind  ! input index 
     39      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in   ) ::   tab2d    ! input 2D field 
     40      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(  out) ::   tab1d    ! output 1D field 
     41      ! 
     42      INTEGER ::   jl, jn, jid, jjd 
     43      !!---------------------------------------------------------------------- 
     44      DO jl = 1, jpl 
     45         DO jn = 1, ndim1d 
     46            jid          = MOD( tab_ind(jn) - 1 , jpi ) + 1 
     47            jjd          =    ( tab_ind(jn) - 1 ) / jpi + 1 
     48            tab1d(jn,jl) = tab2d(jid,jjd,jl) 
     49         END DO 
     50      END DO 
     51   END SUBROUTINE tab_3d_2d 
    2852 
    2953   SUBROUTINE tab_2d_1d( ndim1d, tab_ind, tab1d, tab2d ) 
     
    4569   END SUBROUTINE tab_2d_1d 
    4670 
     71   SUBROUTINE tab_2d_3d( ndim1d, tab_ind, tab1d, tab2d ) 
     72      !!---------------------------------------------------------------------- 
     73      !!                  ***  ROUTINE tab_2d_1d  *** 
     74      !!---------------------------------------------------------------------- 
     75      INTEGER                         , INTENT(in   ) ::   ndim1d    ! 1D size 
     76      INTEGER , DIMENSION(ndim1d)     , INTENT(in   ) ::   tab_ind   ! input index 
     77      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in   ) ::   tab1d     ! input 1D field 
     78      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   tab2d     ! output 2D field 
     79      ! 
     80      INTEGER ::   jl, jn, jid, jjd 
     81      !!---------------------------------------------------------------------- 
     82      DO jl = 1, jpl 
     83         DO jn = 1, ndim1d 
     84            jid               = MOD( tab_ind(jn) - 1 ,  jpi ) + 1 
     85            jjd               =    ( tab_ind(jn) - 1 ) / jpi  + 1 
     86            tab2d(jid,jjd,jl) = tab1d(jn,jl) 
     87         END DO 
     88      END DO 
     89   END SUBROUTINE tab_2d_3d 
    4790 
    4891   SUBROUTINE tab_1d_2d( ndim1d, tab_ind, tab1d, tab2d ) 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r8355 r8369  
    237237         IF( nidx > 0 ) THEN  ! If there is no ice, do nothing. 
    238238            !                                                                 
    239             s_i_new   (:) = 0._wp ; dh_s_tot (:) = 0._wp                     ! --- some init --- ! 
     239            s_i_new   (:) = 0._wp ; dh_s_tot (:) = 0._wp            ! --- some init --- ! 
    240240            dh_i_surf (:) = 0._wp ; dh_i_bott(:) = 0._wp 
    241241            dh_snowice(:) = 0._wp ; dh_i_sub (:) = 0._wp 
    242242 
    243                               CALL lim_thd_1d2d( jl, 1 )               ! --- Move to 1D arrays --- ! 
    244             ! 
    245             DO jk = 1, nlay_i                                                ! --- Change units from J/m2 to J/m3 --- ! 
    246                WHERE( ht_i_1d(:)>0._wp ) e_i_1d(:,jk) = e_i_1d(:,jk) / (ht_i_1d(:) * a_i_1d(:)) * nlay_i 
     243                              CALL lim_thd_1d2d( jl, 1 )            ! --- Move to 1D arrays --- ! 
     244            ! 
     245            DO jk = 1, nlay_i                                       ! --- Change units from J/m2 to J/m3 --- ! 
     246               WHERE( ht_i_1d(1:nidx)>0._wp ) e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) / (ht_i_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_i 
    247247            ENDDO 
    248248            DO jk = 1, nlay_s 
    249                WHERE( ht_s_1d(:)>0._wp ) e_s_1d(:,jk) = e_s_1d(:,jk) / (ht_s_1d(:) * a_i_1d(:)) * nlay_s 
     249               WHERE( ht_s_1d(1:nidx)>0._wp ) e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) / (ht_s_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_s 
    250250            ENDDO 
    251251            ! 
    252             IF( ln_limdH )    CALL lim_thd_dif                    ! --- Ice/Snow Temperature profile --- ! 
    253             ! 
    254             IF( ln_limdH )    CALL lim_thd_dh                     ! --- Ice/Snow thickness --- !     
     252            IF( ln_limdH )    CALL lim_thd_dif                      ! --- Ice/Snow Temperature profile --- ! 
     253            ! 
     254            IF( ln_limdH )    CALL lim_thd_dh                       ! --- Ice/Snow thickness --- !     
    255255            ! 
    256256            IF( ln_limdH )    CALL lim_thd_ent( e_i_1d(1:nidx,:) )  ! --- Ice enthalpy remapping --- ! 
    257257            ! 
    258                               CALL lim_thd_sal                    ! --- Ice salinity --- !     
    259             ! 
    260                               CALL lim_thd_temp                   ! --- temperature update --- ! 
     258                              CALL lim_thd_sal                      ! --- Ice salinity --- !     
     259            ! 
     260                              CALL lim_thd_temp                     ! --- temperature update --- ! 
    261261            ! 
    262262            IF( ln_limdH ) THEN 
    263263               IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN 
    264                               CALL lim_thd_lam                    ! --- extra lateral melting if monocat --- ! 
     264                              CALL lim_thd_lam                      ! --- extra lateral melting if monocat --- ! 
    265265               END IF 
    266266            END IF 
    267267            ! 
    268             DO jk = 1, nlay_i                                                ! --- Change units from J/m3 to J/m2 --- ! 
    269                e_i_1d(:,jk) = e_i_1d(:,jk) * ht_i_1d(:) * a_i_1d(:) * r1_nlay_i 
     268            DO jk = 1, nlay_i                                       ! --- Change units from J/m3 to J/m2 --- ! 
     269               e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) * ht_i_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_i 
    270270            ENDDO 
    271271            DO jk = 1, nlay_s 
    272                e_s_1d(:,jk) = e_s_1d(:,jk) * ht_s_1d(:) * a_i_1d(:) * r1_nlay_s 
     272               e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) * ht_s_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_s 
    273273            ENDDO 
    274274            ! 
    275                               CALL lim_thd_1d2d( jl, 2 )               ! --- Move to 2D arrays --- ! 
     275                              CALL lim_thd_1d2d( jl, 2 )            ! --- Move to 2D arrays --- ! 
    276276            ! 
    277277            IF( lk_mpp )      CALL mpp_comm_free( ncomm_ice ) !RB necessary ?? 
     
    288288 
    289289      IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limitd_thd_da', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    290       IF( ln_limdA)           CALL lim_thd_da                                ! --- lateral melting --- ! 
     290      IF( ln_limdA)           CALL lim_thd_da                       ! --- lateral melting --- ! 
    291291 
    292292      ! Change thickness to volume 
     
    441441         ! 
    442442         CALL tab_2d_1d( nidx, idxice(1:nidx), qprec_ice_1d(1:nidx), qprec_ice        ) 
    443          CALL tab_2d_1d( nidx, idxice(1:nidx), qevap_ice_1d(1:nidx), qevap_ice(:,:,jl)  ) 
    444443         CALL tab_2d_1d( nidx, idxice(1:nidx), qsr_ice_1d  (1:nidx), qsr_ice(:,:,jl)  ) 
    445444         CALL tab_2d_1d( nidx, idxice(1:nidx), fr1_i0_1d   (1:nidx), fr1_i0           ) 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_da.F90

    r8360 r8369  
    163163             
    164164            DO ji = 1, nidx 
    165                ! decrease of concentration for the category jl 
    166                !    each category contributes to melting in proportion to its concentration 
    167                zda = zda_tot(ji) * a_i_1d(ji) / at_i_1d(ji) 
    168                 
    169                ! Contribution to salt flux 
    170                sfx_lam_1d(ji) = sfx_lam_1d(ji) - rhoic *  ht_i_1d(ji) * zda * sm_i_1d(ji) * r1_rdtice 
    171                 
    172                ! Contribution to heat flux into the ocean [W.m-2], (<0)   
    173                hfx_thd_1d(ji) = hfx_thd_1d(ji) + zda_tot(ji) / at_i_1d(ji) * SUM( e_i_1d(ji,1:nlay_i) + e_s_1d(ji,1) ) * r1_rdtice  
    174                
    175                ! Contribution to mass flux 
    176                wfx_lam_1d(ji) =  wfx_lam_1d(ji) - zda * r1_rdtice * ( rhoic * ht_i_1d(ji) + rhosn * ht_s_1d(ji) ) 
    177             
    178                !! adjust e_i ??? 
    179                e_i_1d(ji,1:nlay_i) = e_i_1d(ji,1:nlay_i) * ( 1._wp + zda_tot(ji) / at_i_1d(ji) ) 
    180                e_s_1d(ji,1)        = e_s_1d(ji,1)        * ( 1._wp + zda_tot(ji) / at_i_1d(ji) ) 
    181   
    182                ! new concentration 
    183                a_i_1d(ji) = a_i_1d(ji) + zda 
    184  
    185                ! ensure that ht_i = 0 where a_i = 0 
    186                IF( a_i_1d(ji) == 0._wp ) THEN 
    187                   ht_i_1d(ji) = 0._wp 
    188                   ht_s_1d(ji) = 0._wp 
    189                ENDIF             
    190             END DO 
    191  
     165               IF( a_i_1d(ji) > 0._wp ) THEN 
     166                  ! decrease of concentration for the category jl 
     167                  !    each category contributes to melting in proportion to its concentration 
     168                  zda = zda_tot(ji) * a_i_1d(ji) / at_i_1d(ji) 
     169                   
     170                  ! Contribution to salt flux 
     171                  sfx_lam_1d(ji) = sfx_lam_1d(ji) - rhoic *  ht_i_1d(ji) * zda * sm_i_1d(ji) * r1_rdtice 
     172                   
     173                  ! Contribution to heat flux into the ocean [W.m-2], (<0)   
     174                  hfx_thd_1d(ji) = hfx_thd_1d(ji) + zda_tot(ji) / at_i_1d(ji) * SUM( e_i_1d(ji,1:nlay_i) + e_s_1d(ji,1) ) * r1_rdtice  
     175                   
     176                  ! Contribution to mass flux 
     177                  wfx_lam_1d(ji) =  wfx_lam_1d(ji) - zda * r1_rdtice * ( rhoic * ht_i_1d(ji) + rhosn * ht_s_1d(ji) ) 
     178                   
     179                  !! adjust e_i ??? 
     180                  e_i_1d(ji,1:nlay_i) = e_i_1d(ji,1:nlay_i) * ( 1._wp + zda_tot(ji) / at_i_1d(ji) ) 
     181                  e_s_1d(ji,1)        = e_s_1d(ji,1)        * ( 1._wp + zda_tot(ji) / at_i_1d(ji) ) 
     182                   
     183                  ! new concentration 
     184                  a_i_1d(ji) = a_i_1d(ji) + zda 
     185                   
     186                  ! ensure that ht_i = 0 where a_i = 0 
     187                  IF( a_i_1d(ji) == 0._wp ) THEN 
     188                     ht_i_1d(ji) = 0._wp 
     189                     ht_s_1d(ji) = 0._wp 
     190                  ENDIF 
     191               ENDIF 
     192            END DO 
     193             
    192194            CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d    (1:nidx), a_i (:,:,jl) ) 
    193195            CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d   (1:nidx), ht_i(:,:,jl) ) 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r8342 r8369  
    631631         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN 
    632632            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m     into the ocean 
    633                &                            - sm_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  rn_icesal from the ocean  
     633               &                            - sm_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  rn_icesal from the ocean  
    634634         ENDIF 
    635635 
     
    665665 
    666666      ! --- ensure that a_i = 0 where ht_i = 0 --- 
    667       WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 
    668        
     667      DO ji = 1, nidx 
     668         IF( ht_i_1d(ji) == 0._wp )   a_i_1d(ji) = 0._wp 
     669      END DO 
     670          
    669671      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema ) 
    670672      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zeh_i ) 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r8342 r8369  
    735735      DO ji = 1, nidx 
    736736         zfac        = 1. / MAX( epsi10 , rn_cdsn * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) ) 
    737          t_si_1d(ji) = ( rn_cdsn        * zh_i(ji) * t_s_1d(ji,1) + & 
    738             &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) * zfac 
    739       END DO 
    740       WHERE( ( zh_s < 1.0e-3 ) ) ; t_si_1d(:) = t_su_1d(:) ; END WHERE 
     737         IF( zh_s(ji) >= 1.e-3 ) THEN 
     738            t_si_1d(ji) = ( rn_cdsn        * zh_i(ji) * t_s_1d(ji,1) + & 
     739               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) * zfac 
     740         ELSE 
     741            t_si_1d(ji) = t_su_1d(ji) 
     742         ENDIF 
     743      END DO 
    741744      ! END MV SIMIP 2016 
    742745 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90

    r8342 r8369  
    1212   USE in_out_manager ! I/O manager 
    1313   USE lib_mpp        ! MPP library 
    14    USE ice, ONLY :   nlay_i, nlay_s 
     14   USE ice, ONLY :   nlay_i, nlay_s, jpl 
    1515 
    1616   IMPLICIT NONE 
     
    1919   PUBLIC thd_ice_alloc ! Routine called by nemogcm.F90 
    2020 
    21    !!----------------------------- 
    22    !! * Share 1D Module variables 
    23    !!----------------------------- 
     21   !!---------------------- 
     22   !! * 1D Module variables 
     23   !!---------------------- 
    2424   !: In ice thermodynamics, to spare memory, the vectors are folded 
    2525   !: from 1D to 2D vectors. The following variables, with ending _1d 
     
    9090   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   evap_ice_1d   !: <==> the 2D  evap_ice 
    9191   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qprec_ice_1d  !: <==> the 2D  qprec_ice 
    92    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qevap_ice_1d  !: <==> the 3D  qevap_ice 
    9392   !                                                     ! to reintegrate longwave flux inside the ice thermodynamics 
    9493   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   i0            !: fraction of radiation transmitted to the ice 
     
    9897   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_si_1d       !: <==> the 2D  t_si 
    9998   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_i_1d        !: <==> the 2D  a_i 
    100    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_i_1d       !: <==> the 2D  ht_s 
    101    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_s_1d       !: <==> the 2D  ht_i 
     99   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ib_1d       !: <==> the 2D  a_i_b 
     100   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_i_1d       !: <==> the 2D  ht_i 
     101   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_ib_1d      !: <==> the 2D  ht_i_b 
     102   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_s_1d       !: <==> the 2D  ht_s 
    102103   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_su         !: Surface Conduction flux  
    103104   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_bo_i       !: Bottom  Conduction flux  
     
    109110   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sm_i_1d       !: Ice bulk salinity [ppt] 
    110111   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_new       !: Salinity of new ice at the bottom 
     112   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_1d        !: 
     113   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_ip_1d        !: 
     114   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_i_1d        !: 
    111115 
    112116   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d   !: corresponding to the 2D var  t_s 
     
    127131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sss_1d 
    128132 
     133   !  
     134   !!---------------------- 
     135   !! * 2D Module variables 
     136   !!---------------------- 
     137   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   a_i_2d  
     138   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_i_2d  
     139   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_s_2d  
     140   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   oa_i_2d  
     141   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   smv_i_2d  
     142   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   a_ip_2d 
     143   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_ip_2d  
     144   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_su_2d  
     145   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_i_2d 
     146 
     147   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   a_ib_2d 
     148   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_ib_2d 
     149    
    129150   INTEGER , PUBLIC ::   jiindex_1d   ! 1D index of debugging point 
    130151 
     
    141162      !!---------------------------------------------------------------------! 
    142163      INTEGER ::   thd_ice_alloc   ! return value 
    143       INTEGER ::   ierr(6), ii 
     164      INTEGER ::   ierr(7), ii 
    144165      !!---------------------------------------------------------------------! 
    145166      ierr(:) = 0 
     
    164185         &      wfx_snw_sub_1d(jpij), wfx_ice_sub_1d(jpij), wfx_err_sub_1d(jpij) ,              & 
    165186         &      wfx_lam_1d(jpij)  , dqns_ice_1d(jpij) , evap_ice_1d (jpij),                     & 
    166          &      qprec_ice_1d(jpij), qevap_ice_1d(jpij), i0         (jpij) ,                     &   
     187         &      qprec_ice_1d(jpij), i0         (jpij) ,                                         &   
    167188         &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  & 
    168189         &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  & 
     
    170191      ! 
    171192      ii = ii + 1 
    172       ALLOCATE( t_su_1d   (jpij) , t_si_1d  (jpij)  , a_i_1d   (jpij) , ht_i_1d  (jpij) ,    & 
    173          &      ht_s_1d   (jpij) , fc_su     (jpij) , fc_bo_i  (jpij) ,                      &     
    174          &      dh_s_tot  (jpij) , dh_i_surf (jpij) , dh_i_sub (jpij) ,                      &     
    175          &      dh_i_bott (jpij) , dh_snowice(jpij) , sm_i_1d  (jpij) , s_i_new  (jpij) ,    & 
    176          &      STAT=ierr(ii) ) 
     193      ALLOCATE( t_su_1d  (jpij) , t_si_1d   (jpij) , a_i_1d  (jpij) , a_ib_1d(jpij) ,                  & 
     194         &      ht_i_1d  (jpij) , ht_ib_1d  (jpij) , ht_s_1d (jpij) , fc_su  (jpij) , fc_bo_i(jpij) ,  &     
     195         &      dh_s_tot (jpij) , dh_i_surf (jpij) , dh_i_sub(jpij) ,                                  &     
     196         &      dh_i_bott(jpij) , dh_snowice(jpij) , sm_i_1d (jpij) , s_i_new(jpij) , & 
     197         &      a_ip_1d  (jpij) , v_ip_1d   (jpij) , v_i_1d  (jpij) , STAT=ierr(ii) ) 
    177198      ! 
    178199      ii = ii + 1 
     
    186207      ii = ii + 1 
    187208      ALLOCATE( sst_1d(jpij) , sss_1d(jpij) , STAT=ierr(ii) ) 
     209      ! 
     210      ii = ii + 1 
     211      ALLOCATE( a_i_2d(jpij,jpl) , a_ib_2d(jpij,jpl) , ht_i_2d(jpij,jpl) , ht_ib_2d(jpij,jpl) , & 
     212         &      v_i_2d(jpij,jpl) ,v_s_2d(jpij,jpl) ,oa_i_2d(jpij,jpl) ,smv_i_2d(jpij,jpl) ,  & 
     213         &      a_ip_2d(jpij,jpl) ,v_ip_2d(jpij,jpl) ,t_su_2d(jpij,jpl) ,  & 
     214         &      STAT=ierr(ii) ) 
    188215 
    189216      thd_ice_alloc = MAXVAL( ierr(:) ) 
Note: See TracChangeset for help on using the changeset viewer.