Ignore:
Timestamp:
2015-06-19T17:18:00+02:00 (5 years ago)
Author:
davestorkey
Message:

Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r5234 r5443  
    1313   !!   'key_lim3' :                                   LIM3 sea-ice model 
    1414   !!---------------------------------------------------------------------- 
    15    !!   lim_itd_th       : thermodynamics of ice thickness distribution 
    1615   !!   lim_itd_th_rem   : 
    1716   !!   lim_itd_th_reb   : 
     
    2524   USE thd_ice          ! LIM-3 thermodynamic variables 
    2625   USE ice              ! LIM-3 variables 
    27    USE par_ice          ! LIM-3 parameters 
    28    USE limthd_lac       ! LIM-3 lateral accretion 
    2926   USE limvar           ! LIM-3 variables 
    30    USE limcons          ! LIM-3 conservation 
    3127   USE prtctl           ! Print control 
    3228   USE in_out_manager   ! I/O manager 
     
    3430   USE wrk_nemo         ! work arrays 
    3531   USE lib_fortran      ! to use key_nosignedzero 
    36    USE timing          ! Timing 
    37    USE limcons        ! conservation tests 
     32   USE limcons          ! conservation tests 
    3833 
    3934   IMPLICIT NONE 
    4035   PRIVATE 
    4136 
    42    PUBLIC   lim_itd_th         ! called by ice_stp 
    4337   PUBLIC   lim_itd_th_rem 
    4438   PUBLIC   lim_itd_th_reb 
    45    PUBLIC   lim_itd_fitline 
    46    PUBLIC   lim_itd_shiftice 
    4739 
    4840   !!---------------------------------------------------------------------- 
     
    5345CONTAINS 
    5446 
    55    SUBROUTINE lim_itd_th( kt ) 
    56       !!------------------------------------------------------------------ 
    57       !!                ***  ROUTINE lim_itd_th *** 
    58       !! 
    59       !! ** Purpose :   computes the thermodynamics of ice thickness distribution 
    60       !! 
    61       !! ** Method  : 
    62       !!------------------------------------------------------------------ 
    63       INTEGER, INTENT(in) ::   kt   ! time step index 
    64       ! 
    65       INTEGER ::   ji, jj, jk, jl   ! dummy loop index          
    66       ! 
    67       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    68       !!------------------------------------------------------------------ 
    69       IF( nn_timing == 1 )  CALL timing_start('limitd_th') 
    70  
    71       ! conservation test 
    72       IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    73  
    74       IF( kt == nit000 .AND. lwp ) THEN 
    75          WRITE(numout,*) 
    76          WRITE(numout,*) 'lim_itd_th  : Thermodynamics of the ice thickness distribution' 
    77          WRITE(numout,*) '~~~~~~~~~~~' 
    78       ENDIF 
    79  
    80       !------------------------------------------------------------------------------| 
    81       !  1) Transport of ice between thickness categories.                           | 
    82       !------------------------------------------------------------------------------| 
    83       ! Given thermodynamic growth rates, transport ice between 
    84       ! thickness categories. 
    85       IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt ) 
    86       ! 
    87       CALL lim_var_glo2eqv    ! only for info 
    88       CALL lim_var_agg(1) 
    89  
    90       !------------------------------------------------------------------------------| 
    91       !  3) Add frazil ice growing in leads. 
    92       !------------------------------------------------------------------------------| 
    93       CALL lim_thd_lac 
    94       CALL lim_var_glo2eqv    ! only for info 
    95       
    96       IF(ln_ctl) THEN   ! Control print 
    97          CALL prt_ctl_info(' ') 
    98          CALL prt_ctl_info(' - Cell values : ') 
    99          CALL prt_ctl_info('   ~~~~~~~~~~~~~ ') 
    100          CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_th  : cell area :') 
    101          CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th  : at_i      :') 
    102          CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th  : vt_i      :') 
    103          CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th  : vt_s      :') 
    104          DO jl = 1, jpl 
    105             CALL prt_ctl_info(' ') 
    106             CALL prt_ctl_info(' - Category : ', ivar1=jl) 
    107             CALL prt_ctl_info('   ~~~~~~~~~~') 
    108             CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : a_i      : ') 
    109             CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_i     : ') 
    110             CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_s     : ') 
    111             CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_i      : ') 
    112             CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_s      : ') 
    113             CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : e_s      : ') 
    114             CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_itd_th  : t_su     : ') 
    115             CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : t_snow   : ') 
    116             CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : sm_i     : ') 
    117             CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_th  : smv_i    : ') 
    118             DO jk = 1, nlay_i 
    119                CALL prt_ctl_info(' ') 
    120                CALL prt_ctl_info(' - Layer : ', ivar1=jk) 
    121                CALL prt_ctl_info('   ~~~~~~~') 
    122                CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : t_i      : ') 
    123                CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : e_i      : ') 
    124             END DO 
    125          END DO 
    126       ENDIF 
    127       ! 
    128       ! conservation test 
    129       IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    130       ! 
    131      IF( nn_timing == 1 )  CALL timing_stop('limitd_th') 
    132    END SUBROUTINE lim_itd_th 
    133    ! 
    134  
    13547   SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) 
    13648      !!------------------------------------------------------------------ 
     
    15365      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars 
    15466      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      - 
    155       REAL(wp) ::   zx3,             zareamin          !   -      - 
     67      REAL(wp) ::   zx3         
    15668      CHARACTER (len = 15) :: fieldid 
    15769 
     
    17991      !!------------------------------------------------------------------ 
    18092 
    181       CALL wrk_alloc( jpi,jpj, zremap_flag )    ! integer 
    182       CALL wrk_alloc( jpi,jpj,jpl-1, zdonor )   ! integer 
     93      CALL wrk_alloc( jpi,jpj, zremap_flag ) 
     94      CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) 
    18395      CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 
    18496      CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    18597      CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
    18698      CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )    
    187       CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j )   ! integer  
     99      CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j )  
    188100      CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 
    189  
    190       zareamin = epsi10   !minimum area in thickness categories tolerated by the conceptors of the model 
    191101 
    192102      !!---------------------------------------------------------------------------------------------- 
     
    216126         DO jj = 1, jpj 
    217127            DO ji = 1, jpi 
    218                rswitch             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) )     !0 if no ice and 1 if yes 
     128               rswitch           = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) )     !0 if no ice and 1 if yes 
    219129               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 
    220                rswitch             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 
     130               rswitch           = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) ) 
    221131               zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 
    222                IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)  
     132               IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) ! clem: useless IF statement?  
    223133            END DO 
    224134         END DO 
     
    239149      DO jj = 1, jpj 
    240150         DO ji = 1, jpi 
    241             IF ( at_i(ji,jj) .gt. zareamin ) THEN 
     151            IF ( at_i(ji,jj) > epsi10 ) THEN 
    242152               nbrem         = nbrem + 1 
    243153               nind_i(nbrem) = ji 
     
    247157               zremap_flag(ji,jj) = 0 
    248158            ENDIF 
    249          END DO !ji 
    250       END DO !jj 
     159         END DO 
     160      END DO 
    251161 
    252162      !----------------------------------------------------------------------------------------------- 
     
    254164      !----------------------------------------------------------------------------------------------- 
    255165      !- 4.1 Compute category boundaries 
    256       ! Tricky trick see limitd_me.F90 
    257       ! will be soon removed, CT 
    258       ! hi_max(kubnd) = 99. 
    259166      zhbnew(:,:,:) = 0._wp 
    260167 
     
    265172            ! 
    266173            zhbnew(ii,ij,jl) = hi_max(jl) 
    267             IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 
     174            IF    ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 
    268175               !interpolate between adjacent category growth rates 
    269176               zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 
    270177               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 
    271             ELSEIF ( a_i_b(ii,ij,jl) > epsi10) THEN 
     178            ELSEIF( a_i_b(ii,ij,jl) > epsi10) THEN 
    272179               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 
    273             ELSEIF ( a_i_b(ii,ij,jl+1) > epsi10) THEN 
     180            ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN 
    274181               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 
    275182            ENDIF 
     
    280187            ii = nind_i(ji) 
    281188            ij = nind_j(ji) 
    282             IF( a_i(ii,ij,jl) > epsi10 .AND. ht_i(ii,ij,jl) >= zhbnew(ii,ij,jl) ) THEN 
     189 
     190            ! clem: we do not want ht_i to be too close to either HR or HL otherwise a division by nearly 0 is possible  
     191            ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
     192            IF    ( a_i(ii,ij,jl  ) > epsi10 .AND. ht_i(ii,ij,jl  ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) THEN 
    283193               zremap_flag(ii,ij) = 0 
    284             ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) <= zhbnew(ii,ij,jl) ) THEN 
     194            ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) THEN 
    285195               zremap_flag(ii,ij) = 0 
    286196            ENDIF 
    287197 
    288198            !- 4.3 Check that each zhbnew does not exceed maximal values hi_max   
     199            IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 
    289200            IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 
    290             IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 
    291          END DO 
    292  
    293       END DO !jl 
     201            ! clem bug: why is not the following instead? 
     202            !!IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 
     203            !!IF( zhbnew(ii,ij,jl) > hi_max(jl  ) ) zremap_flag(ii,ij) = 0 
     204  
     205         END DO 
     206 
     207      END DO 
    294208 
    295209      !----------------------------------------------------------------------------------------------- 
     
    312226      DO jj = 1, jpj 
    313227         DO ji = 1, jpi 
    314             zhb0(ji,jj) = hi_max(0) ! 0eme 
    315             zhb1(ji,jj) = hi_max(1) ! 1er 
    316  
    317             zhbnew(ji,jj,klbnd-1) = 0._wp 
     228            zhb0(ji,jj) = hi_max(0) 
     229            zhb1(ji,jj) = hi_max(1) 
    318230 
    319231            IF( a_i(ji,jj,kubnd) > epsi10 ) THEN 
    320                zhbnew(ji,jj,kubnd) = 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) 
     232               zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) ) 
    321233            ELSE 
    322                zhbnew(ji,jj,kubnd) = hi_max(kubnd)   
    323                !!? clem bug: since hi_max(jpl)=99, this limit is very high  
    324                !!? but I think it is erased in fitline subroutine  
    325             ENDIF 
    326  
    327             IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 
    328  
    329          END DO !jj 
    330       END DO !jj 
     234!clem bug               zhbnew(ji,jj,kubnd) = hi_max(kubnd)   
     235               zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) ! not used anyway 
     236            ENDIF 
     237 
     238            ! clem: we do not want ht_i_b to be too close to either HR or HL otherwise a division by nearly 0 is possible  
     239            ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
     240            IF    ( zht_i_b(ji,jj,klbnd) < ( zhb0(ji,jj) + epsi10 ) )  THEN 
     241               zremap_flag(ji,jj) = 0 
     242            ELSEIF( zht_i_b(ji,jj,klbnd) > ( zhb1(ji,jj) - epsi10 ) )  THEN 
     243               zremap_flag(ji,jj) = 0 
     244            ENDIF 
     245 
     246         END DO 
     247      END DO 
    331248 
    332249      !----------------------------------------------------------------------------------------------- 
     
    334251      !----------------------------------------------------------------------------------------------- 
    335252      !- 7.1 g(h) for category 1 at start of time step 
    336       CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd),         & 
    337          &                  g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),   & 
     253      CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),   & 
    338254         &                  hR(:,:,klbnd), zremap_flag ) 
    339255 
     
    343259         ij = nind_j(ji)  
    344260 
    345          !ji 
    346          IF (a_i(ii,ij,klbnd) .gt. epsi10) THEN 
     261         IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 
     262 
    347263            zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 
    348             ! ji, a_i > epsi10 
    349             IF (zdh0 .lt. 0.0) THEN !remove area from category 1 
    350                ! ji, a_i > epsi10; zdh0 < 0 
    351                zdh0 = MIN(-zdh0,hi_max(klbnd)) 
    352  
     264 
     265            IF( zdh0 < 0.0 ) THEN      !remove area from category 1 
     266               zdh0 = MIN( -zdh0, hi_max(klbnd) ) 
    353267               !Integrate g(1) from 0 to dh0 to estimate area melted 
    354                zetamax = MIN(zdh0,hR(ii,ij,klbnd)) - hL(ii,ij,klbnd) 
    355                IF (zetamax.gt.0.0) THEN 
    356                   zx1  = zetamax 
    357                   zx2  = 0.5 * zetamax*zetamax  
    358                   zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 
    359                   ! Constrain new thickness <= ht_i 
    360                   zdamax = a_i(ii,ij,klbnd) * &  
    361                      (1.0 - ht_i(ii,ij,klbnd)/zht_i_b(ii,ij,klbnd)) ! zdamax > 0 
    362                   !ice area lost due to melting of thin ice 
    363                   zda0   = MIN(zda0, zdamax) 
    364  
     268               zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd) 
     269 
     270               IF( zetamax > 0.0 ) THEN 
     271                  zx1    = zetamax 
     272                  zx2    = 0.5 * zetamax * zetamax  
     273                  zda0   = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1                        ! ice area removed 
     274                  zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! Constrain new thickness <= ht_i                 
     275                  zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting  
     276                                                                                                !     of thin ice (zdamax > 0) 
    365277                  ! Remove area, conserving volume 
    366                   ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) &  
    367                      * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 
     278                  ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 
    368279                  a_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) - zda0 
    369                   v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem-useless ? 
    370                ENDIF     ! zetamax > 0 
    371                ! ji, a_i > epsi10 
    372  
    373             ELSE ! if ice accretion 
    374                ! ji, a_i > epsi10; zdh0 > 0 
    375                zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd))  
    376                ! zhbnew was 0, and is shifted to the right to account for thin ice 
    377                ! growth in openwater (F0 = f1) 
    378             ENDIF ! zdh0  
    379  
    380             ! a_i > epsi10 
    381          ENDIF ! a_i > epsi10 
    382  
    383       END DO ! ji 
     280                  v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ? 
     281               ENDIF 
     282 
     283            ELSE ! if ice accretion zdh0 > 0 
     284               ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 
     285               zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) )  
     286            ENDIF 
     287 
     288         ENDIF 
     289 
     290      END DO 
    384291 
    385292      !- 7.3 g(h) for each thickness category   
    386293      DO jl = klbnd, kubnd 
    387          CALL lim_itd_fitline(jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 
    388             g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag) 
     294         CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 
     295            &                  g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) 
    389296      END DO 
    390297 
     
    406313            ij = nind_j(ji) 
    407314 
    408             IF (zhbnew(ii,ij,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1 
    409  
     315            IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1 
    410316               ! left and right integration limits in eta space 
    411                zvetamin(ji) = MAX(hi_max(jl), hL(ii,ij,jl)) - hL(ii,ij,jl) 
    412                zvetamax(ji) = MIN(zhbnew(ii,ij,jl), hR(ii,ij,jl)) - hL(ii,ij,jl) 
     317               zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl) 
     318               zvetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) 
    413319               zdonor(ii,ij,jl) = jl 
    414320 
    415             ELSE  ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 
    416  
     321            ELSE                                    ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 
    417322               ! left and right integration limits in eta space 
    418323               zvetamin(ji) = 0.0 
    419                zvetamax(ji) = MIN(hi_max(jl), hR(ii,ij,jl+1)) - hL(ii,ij,jl+1) 
     324               zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1) 
    420325               zdonor(ii,ij,jl) = jl + 1 
    421326 
    422             ENDIF  ! zhbnew(jl) > hi_max(jl) 
    423  
    424             zetamax = MAX(zvetamax(ji), zvetamin(ji)) ! no transfer if etamax < etamin 
     327            ENDIF 
     328 
     329            zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin 
    425330            zetamin = zvetamin(ji) 
    426331 
    427332            zx1  = zetamax - zetamin 
    428             zwk1 = zetamin*zetamin 
    429             zwk2 = zetamax*zetamax 
    430             zx2  = 0.5 * (zwk2 - zwk1) 
     333            zwk1 = zetamin * zetamin 
     334            zwk2 = zetamax * zetamax 
     335            zx2  = 0.5 * ( zwk2 - zwk1 ) 
    431336            zwk1 = zwk1 * zetamin 
    432337            zwk2 = zwk2 * zetamax 
    433             zx3  = 1.0/3.0 * (zwk2 - zwk1) 
     338            zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 ) 
    434339            nd   = zdonor(ii,ij,jl) 
    435340            zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 
    436341            zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 
    437342 
    438          END DO ! ji 
    439       END DO ! jl klbnd -> kubnd - 1 
     343         END DO 
     344      END DO 
    440345 
    441346      !!---------------------------------------------------------------------------------------------- 
     
    451356         ii = nind_i(ji) 
    452357         ij = nind_j(ji) 
    453          IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < hiclim ) THEN 
    454             a_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim  
    455             ht_i(ii,ij,1) = hiclim 
     358         IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN 
     359            a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin  
     360            ht_i(ii,ij,1) = rn_himin 
    456361         ENDIF 
    457       END DO !ji 
     362      END DO 
    458363 
    459364      !!---------------------------------------------------------------------------------------------- 
     
    479384      ENDIF 
    480385 
    481       CALL wrk_dealloc( jpi,jpj, zremap_flag )    ! integer 
    482       CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor )   ! integer 
     386      CALL wrk_dealloc( jpi,jpj, zremap_flag ) 
     387      CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) 
    483388      CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 
    484389      CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    485390      CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
    486391      CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )    
    487       CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j )   ! integer  
     392      CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j )  
    488393      CALL wrk_dealloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 
    489394 
     
    491396 
    492397 
    493    SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice,   & 
    494       &                        g0, g1, hL, hR, zremap_flag ) 
     398   SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 
    495399      !!------------------------------------------------------------------ 
    496400      !!                ***  ROUTINE lim_itd_fitline *** 
     
    511415      INTEGER , DIMENSION(jpi,jpj), INTENT(in   ) ::   zremap_flag  ! 
    512416      ! 
    513       INTEGER ::   ji,jj           ! horizontal indices 
     417      INTEGER  ::   ji,jj        ! horizontal indices 
    514418      REAL(wp) ::   zh13         ! HbL + 1/3 * (HbR - HbL) 
    515419      REAL(wp) ::   zh23         ! HbL + 2/3 * (HbR - HbL) 
     
    518422      !!------------------------------------------------------------------ 
    519423      ! 
    520       ! 
    521424      DO jj = 1, jpj 
    522425         DO ji = 1, jpi 
    523426            ! 
    524427            IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10   & 
    525                &                        .AND. hice(ji,jj)        > 0._wp     ) THEN 
     428               &                        .AND. hice(ji,jj)        > 0._wp ) THEN 
    526429 
    527430               ! Initialize hL and hR 
    528  
    529431               hL(ji,jj) = HbL(ji,jj) 
    530432               hR(ji,jj) = HbR(ji,jj) 
    531433 
    532434               ! Change hL or hR if hice falls outside central third of range 
    533  
    534                zh13 = 1.0/3.0 * (2.0*hL(ji,jj) + hR(ji,jj)) 
    535                zh23 = 1.0/3.0 * (hL(ji,jj) + 2.0*hR(ji,jj)) 
     435               zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) ) 
     436               zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) ) 
    536437 
    537438               IF    ( hice(ji,jj) < zh13 ) THEN   ;   hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) 
     
    540441 
    541442               ! Compute coefficients of g(eta) = g0 + g1*eta 
    542  
    543443               zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj)) 
    544444               zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 
    545445               zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr 
    546                g0(ji,jj) = zwk1 * ( 2._wp/3._wp - zwk2 ) 
    547                g1(ji,jj) = 2._wp * zdhr * zwk1 * (zwk2 - 0.5) 
     446               g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 ) 
     447               g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) 
    548448               ! 
    549             ELSE                   ! remap_flag = .false. or a_i < epsi10  
     449            ELSE  ! remap_flag = .false. or a_i < epsi10  
    550450               hL(ji,jj) = 0._wp 
    551451               hR(ji,jj) = 0._wp 
    552452               g0(ji,jj) = 0._wp 
    553453               g1(ji,jj) = 0._wp 
    554             ENDIF                  ! a_i > epsi10 
     454            ENDIF 
    555455            ! 
    556456         END DO 
     
    576476 
    577477      INTEGER ::   ji, jj, jl, jl2, jl1, jk   ! dummy loop indices 
    578       INTEGER ::   ii, ij          ! indices when changing from 2D-1D is done 
     478      INTEGER ::   ii, ij                     ! indices when changing from 2D-1D is done 
    579479 
    580480      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zaTsfn 
     
    589489      INTEGER, POINTER, DIMENSION(:) ::   nind_i, nind_j   ! compressed indices for i/j directions 
    590490 
    591       INTEGER ::   nbrem             ! number of cells with ice to transfer 
    592  
    593       LOGICAL ::   zdaice_negative         ! true if daice < -puny 
    594       LOGICAL ::   zdvice_negative         ! true if dvice < -puny 
    595       LOGICAL ::   zdaice_greater_aicen    ! true if daice > aicen 
    596       LOGICAL ::   zdvice_greater_vicen    ! true if dvice > vicen 
     491      INTEGER  ::   nbrem             ! number of cells with ice to transfer 
    597492      !!------------------------------------------------------------------ 
    598493 
    599494      CALL wrk_alloc( jpi,jpj,jpl, zaTsfn ) 
    600495      CALL wrk_alloc( jpi,jpj, zworka ) 
    601       CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j )   ! integer 
     496      CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 
    602497 
    603498      !---------------------------------------------------------------------------------------------- 
     
    606501 
    607502      DO jl = klbnd, kubnd 
    608          zaTsfn(:,:,jl) = a_i(:,:,jl)*t_su(:,:,jl) 
    609       END DO 
    610  
    611       !---------------------------------------------------------------------------------------------- 
    612       ! 2) Check for daice or dvice out of range, allowing for roundoff error 
    613       !---------------------------------------------------------------------------------------------- 
    614       ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl 
    615       ! has a small area, with h(n) very close to a boundary.  Then 
    616       ! the coefficients of g(h) are large, and the computed daice and 
    617       ! dvice can be in error. If this happens, it is best to transfer 
    618       ! either the entire category or nothing at all, depending on which 
    619       ! side of the boundary hice(n) lies. 
    620       !----------------------------------------------------------------- 
    621       DO jl = klbnd, kubnd-1 
    622  
    623          zdaice_negative = .false. 
    624          zdvice_negative = .false. 
    625          zdaice_greater_aicen = .false. 
    626          zdvice_greater_vicen = .false. 
    627  
    628          DO jj = 1, jpj 
    629             DO ji = 1, jpi 
    630  
    631                IF (zdonor(ji,jj,jl) .GT. 0) THEN 
    632                   jl1 = zdonor(ji,jj,jl) 
    633  
    634                   IF (zdaice(ji,jj,jl) .LT. 0.0) THEN 
    635                      IF (zdaice(ji,jj,jl) .GT. -epsi10) THEN 
    636                         IF ( ( jl1.EQ.jl   .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) )           & 
    637                            .OR.                                      & 
    638                            ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) )           &   
    639                            ) THEN                                                              
    640                            zdaice(ji,jj,jl) = a_i(ji,jj,jl1)  ! shift entire category 
    641                            zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
    642                         ELSE 
    643                            zdaice(ji,jj,jl) = 0.0 ! shift no ice 
    644                            zdvice(ji,jj,jl) = 0.0 
    645                         ENDIF 
    646                      ELSE 
    647                         zdaice_negative = .true. 
    648                      ENDIF 
    649                   ENDIF 
    650  
    651                   IF (zdvice(ji,jj,jl) .LT. 0.0) THEN 
    652                      IF (zdvice(ji,jj,jl) .GT. -epsi10 ) THEN 
    653                         IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) )     & 
    654                            .OR.                                     & 
    655                            ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 
    656                            ) THEN 
    657                            zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 
    658                            zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
    659                         ELSE 
    660                            zdaice(ji,jj,jl) = 0.0    ! shift no ice 
    661                            zdvice(ji,jj,jl) = 0.0 
    662                         ENDIF 
    663                      ELSE 
    664                         zdvice_negative = .true. 
    665                      ENDIF 
    666                   ENDIF 
    667  
    668                   ! If daice is close to aicen, set daice = aicen. 
    669                   IF (zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - epsi10 ) THEN 
    670                      IF (zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1)+epsi10) THEN 
    671                         zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 
    672                         zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
    673                      ELSE 
    674                         zdaice_greater_aicen = .true. 
    675                      ENDIF 
    676                   ENDIF 
    677  
    678                   IF (zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1)-epsi10) THEN 
    679                      IF (zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1)+epsi10) THEN 
    680                         zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 
    681                         zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
    682                      ELSE 
    683                         zdvice_greater_vicen = .true. 
    684                      ENDIF 
    685                   ENDIF 
    686  
    687                ENDIF               ! donor > 0 
    688             END DO                   ! i 
    689          END DO                 ! j 
    690  
    691       END DO !jl 
     503         zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) 
     504      END DO 
    692505 
    693506      !------------------------------------------------------------------------------- 
    694       ! 3) Transfer volume and energy between categories 
     507      ! 2) Transfer volume and energy between categories 
    695508      !------------------------------------------------------------------------------- 
    696509 
     
    699512         DO jj = 1, jpj 
    700513            DO ji = 1, jpi 
    701                IF (zdaice(ji,jj,jl) .GT. 0.0 ) THEN ! daice(n) can be < puny 
     514               IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny 
    702515                  nbrem = nbrem + 1 
    703516                  nind_i(nbrem) = ji 
    704517                  nind_j(nbrem) = jj 
    705                ENDIF ! tmask 
     518               ENDIF 
    706519            END DO 
    707520         END DO 
     
    712525 
    713526            jl1 = zdonor(ii,ij,jl) 
    714             rswitch             = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) ) 
    715             zworka(ii,ij)   = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * rswitch 
     527            rswitch       = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi10 ) ) 
     528            zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch 
    716529            IF( jl1 == jl) THEN   ;   jl2 = jl1+1 
    717             ELSE                    ;   jl2 = jl  
     530            ELSE                  ;   jl2 = jl  
    718531            ENDIF 
    719532 
     
    721534            ! Ice areas 
    722535            !-------------- 
    723  
    724536            a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 
    725537            a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) 
     
    728540            ! Ice volumes 
    729541            !-------------- 
    730  
    731542            v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl)  
    732543            v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) 
     
    735546            ! Snow volumes 
    736547            !-------------- 
    737  
    738548            zdvsnow        = v_s(ii,ij,jl1) * zworka(ii,ij) 
    739549            v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 
     
    743553            ! Snow heat content   
    744554            !-------------------- 
    745  
    746555            zdesnow            = e_s(ii,ij,1,jl1) * zworka(ii,ij) 
    747556            e_s(ii,ij,1,jl1)   = e_s(ii,ij,1,jl1) - zdesnow 
     
    751560            ! Ice age  
    752561            !-------------- 
    753  
    754562            zdo_aice           = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 
    755563            oa_i(ii,ij,jl1)    = oa_i(ii,ij,jl1) - zdo_aice 
     
    759567            ! Ice salinity 
    760568            !-------------- 
    761  
    762569            zdsm_vice          = smv_i(ii,ij,jl1) * zworka(ii,ij) 
    763570            smv_i(ii,ij,jl1)   = smv_i(ii,ij,jl1) - zdsm_vice 
     
    767574            ! Surface temperature 
    768575            !--------------------- 
    769  
    770576            zdaTsf             = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 
    771577            zaTsfn(ii,ij,jl1)  = zaTsfn(ii,ij,jl1) - zdaTsf 
    772578            zaTsfn(ii,ij,jl2)  = zaTsfn(ii,ij,jl2) + zdaTsf  
    773579 
    774          END DO                 ! ji 
     580         END DO 
    775581 
    776582         !------------------ 
     
    779585 
    780586         DO jk = 1, nlay_i 
    781 !CDIR NODEP 
    782587            DO ji = 1, nbrem 
    783588               ii = nind_i(ji) 
     
    785590 
    786591               jl1 = zdonor(ii,ij,jl) 
    787                IF (jl1 .EQ. jl) THEN 
     592               IF (jl1 == jl) THEN 
    788593                  jl2 = jl+1 
    789594               ELSE             ! n1 = n+1 
     
    794599               e_i(ii,ij,jk,jl1) =  e_i(ii,ij,jk,jl1) - zdeice 
    795600               e_i(ii,ij,jk,jl2) =  e_i(ii,ij,jk,jl2) + zdeice  
    796             END DO              ! ji 
    797          END DO                 ! jk 
     601            END DO 
     602         END DO 
    798603 
    799604      END DO                   ! boundaries, 1 to ncat-1 
     
    809614                  ht_i(ji,jj,jl)  =  v_i   (ji,jj,jl) / a_i(ji,jj,jl)  
    810615                  t_su(ji,jj,jl)  =  zaTsfn(ji,jj,jl) / a_i(ji,jj,jl)  
    811                   rswitch         =  1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 
    812616               ELSE 
    813617                  ht_i(ji,jj,jl)  = 0._wp 
    814                   t_su(ji,jj,jl)  = rtt 
     618                  t_su(ji,jj,jl)  = rt0 
    815619               ENDIF 
    816             END DO                 ! ji 
    817          END DO                 ! jj 
    818       END DO                    ! jl 
     620            END DO 
     621         END DO 
     622      END DO 
    819623      ! 
    820624      CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn ) 
    821625      CALL wrk_dealloc( jpi,jpj, zworka ) 
    822       CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j )   ! integer 
     626      CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 
    823627      ! 
    824628   END SUBROUTINE lim_itd_shiftice 
     
    846650      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   ! snow volume summed over categories 
    847651      !!------------------------------------------------------------------ 
    848       !! clem 2014/04: be carefull, rebining does not conserve salt(maybe?) => the difference is taken into account in limupdate 
    849652       
    850653      CALL wrk_alloc( jpi,jpj,jpl, zdonor )   ! interger 
     
    864667         DO jj = 1, jpj 
    865668            DO ji = 1, jpi  
    866                IF( a_i(ji,jj,jl) > epsi10 ) THEN  
    867                   ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    868                ELSE 
    869                   ht_i(ji,jj,jl) = 0._wp 
    870                ENDIF 
     669               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
     670               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 
    871671            END DO 
    872672         END DO 
     
    874674 
    875675      !------------------------------------------------------------------------------ 
    876       ! 2) Make sure thickness of cat klbnd is at least hi_max(klbnd) 
    877       !------------------------------------------------------------------------------ 
    878       DO jj = 1, jpj  
    879          DO ji = 1, jpi  
    880             IF( a_i(ji,jj,klbnd) > epsi10 ) THEN 
    881                IF( ht_i(ji,jj,klbnd) <= hi_max(0) .AND. hi_max(0) > 0._wp ) THEN 
    882                   a_i(ji,jj,klbnd)  = v_i(ji,jj,klbnd) / hi_max(0)  
    883                   ht_i(ji,jj,klbnd) = hi_max(0) 
    884                ENDIF 
    885             ENDIF 
    886          END DO 
    887       END DO 
    888  
    889       !------------------------------------------------------------------------------ 
    890       ! 3) If a category thickness is not in bounds, shift the 
     676      ! 2) If a category thickness is not in bounds, shift the 
    891677      ! entire area, volume, and energy to the neighboring category 
    892678      !------------------------------------------------------------------------------ 
     
    917703                  zdonor(ji,jj,jl)  = jl  
    918704                  ! begin TECLIM change 
    919                   !zdaice(ji,jj,jl)  = a_i(ji,jj,jl) 
    920                   !zdvice(ji,jj,jl)  = v_i(ji,jj,jl) 
    921705                  !zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * 0.5_wp 
    922706                  !zdvice(ji,jj,jl)  = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1)) * 0.5_wp 
    923707                  ! end TECLIM change  
    924708                  ! clem: how much of a_i you send in cat sup is somewhat arbitrary 
    925                   zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi10 ) / ht_i(ji,jj,jl)   
    926                   zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi10 ) 
     709                  zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl)   
     710                  zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi20 ) 
    927711               ENDIF 
    928             END DO                 ! ji 
    929          END DO                 ! jj 
     712            END DO 
     713         END DO 
    930714         IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
    931715 
     
    938722         ENDIF 
    939723         ! 
    940       END DO                    ! jl 
     724      END DO 
    941725 
    942726      !---------------------------- 
     
    951735         zshiftflag = 0 
    952736 
    953 !clem-change 
    954737         DO jj = 1, jpj 
    955738            DO ji = 1, jpi 
     
    961744                  zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 
    962745               ENDIF 
    963             END DO                 ! ji 
    964          END DO                 ! jj 
     746            END DO 
     747         END DO 
    965748 
    966749         IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
     
    973756            zdvice(:,:,jl) = 0._wp 
    974757         ENDIF 
    975 !clem-change 
    976  
    977 !         ! clem-change begin: why not doing that? 
    978 !         DO jj = 1, jpj 
    979 !            DO ji = 1, jpi 
    980 !               IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 
    981 !                  ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10 
    982 !                  a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1)  
    983 !               ENDIF 
    984 !            END DO                 ! ji 
    985 !         END DO                 ! jj 
    986          ! clem-change end 
    987  
    988       END DO                    ! jl 
     758 
     759      END DO 
    989760 
    990761      !------------------------------------------------------------------------------ 
    991       ! 4) Conservation check 
     762      ! 3) Conservation check 
    992763      !------------------------------------------------------------------------------ 
    993764 
     
    1002773      ENDIF 
    1003774      ! 
    1004       CALL wrk_dealloc( jpi,jpj,jpl, zdonor )   ! interger 
     775      CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) 
    1005776      CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice ) 
    1006777      CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final ) 
     
    1013784   !!---------------------------------------------------------------------- 
    1014785CONTAINS 
    1015    SUBROUTINE lim_itd_th           ! Empty routines 
    1016    END SUBROUTINE lim_itd_th 
    1017    SUBROUTINE lim_itd_th_ini 
    1018    END SUBROUTINE lim_itd_th_ini 
    1019786   SUBROUTINE lim_itd_th_rem 
    1020787   END SUBROUTINE lim_itd_th_rem 
Note: See TracChangeset for help on using the changeset viewer.