New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 4921 for branches/2014/dev_r4650_UKMO13_CICE_changes_take2/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 – NEMO

Ignore:
Timestamp:
2014-11-28T14:59:01+01:00 (9 years ago)
Author:
timgraham
Message:

merged with revision 4879 of trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO13_CICE_changes_take2/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r4333 r4921  
    66   !! History :   -   !          (W. H. Lipscomb and E.C. Hunke) CICE (c) original code 
    77   !!            3.0  ! 2005-12  (M. Vancoppenolle) adaptation to LIM-3 
    8    !!             -   ! 2006-06  (M. Vancoppenolle) adaptation to include salt, age and types 
     8   !!             -   ! 2006-06  (M. Vancoppenolle) adaptation to include salt, age 
    99   !!             -   ! 2007-04  (M. Vancoppenolle) Mass conservation checked 
    1010   !!---------------------------------------------------------------------- 
     
    3535   USE lib_fortran      ! to use key_nosignedzero 
    3636   USE timing          ! Timing 
     37   USE limcons        ! conservation tests 
    3738 
    3839   IMPLICIT NONE 
     
    6566      INTEGER, INTENT(in) ::   kt   ! time step index 
    6667      ! 
    67       INTEGER ::   jl, ja, jm, jbnd1, jbnd2   ! ice types    dummy loop index          
    68       REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
    69       REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
     68      INTEGER ::   ji, jj, jk, jl   ! dummy loop index          
     69      ! 
     70      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    7071      !!------------------------------------------------------------------ 
    7172      IF( nn_timing == 1 )  CALL timing_start('limitd_th') 
    7273 
    73       ! ------------------------------- 
    74       !- check conservation (C Rousset) 
    75       IF (ln_limdiahsb) THEN 
    76          zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    77          zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    78          zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 
    79          zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 
    80        ENDIF 
    81       !- check conservation (C Rousset) 
    82       ! ------------------------------- 
     74      ! conservation test 
     75      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    8376 
    8477      IF( kt == nit000 .AND. lwp ) THEN 
     
    9386      ! Given thermodynamic growth rates, transport ice between 
    9487      ! thickness categories. 
    95       DO jm = 1, jpm 
    96          jbnd1 = ice_cat_bounds(jm,1) 
    97          jbnd2 = ice_cat_bounds(jm,2) 
    98          IF( ice_ncat_types(jm) > 1 )   CALL lim_itd_th_rem( jbnd1, jbnd2, jm, kt ) 
    99       END DO 
     88      IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt ) 
    10089      ! 
    10190      CALL lim_var_glo2eqv    ! only for info 
     
    10594      !  3) Add frazil ice growing in leads. 
    10695      !------------------------------------------------------------------------------| 
    107  
    10896      CALL lim_thd_lac 
    10997      CALL lim_var_glo2eqv    ! only for info 
    110  
    111      IF(ln_ctl) THEN   ! Control print 
     98      
     99      IF(ln_ctl) THEN   ! Control print 
    112100         CALL prt_ctl_info(' ') 
    113101         CALL prt_ctl_info(' - Cell values : ') 
     
    131119            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : sm_i     : ') 
    132120            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_th  : smv_i    : ') 
    133             DO ja = 1, nlay_i 
     121            DO jk = 1, nlay_i 
    134122               CALL prt_ctl_info(' ') 
    135                CALL prt_ctl_info(' - Layer : ', ivar1=ja) 
     123               CALL prt_ctl_info(' - Layer : ', ivar1=jk) 
    136124               CALL prt_ctl_info('   ~~~~~~~') 
    137                CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_itd_th  : t_i      : ') 
    138                CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_itd_th  : e_i      : ') 
     125               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : t_i      : ') 
     126               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : e_i      : ') 
    139127            END DO 
    140128         END DO 
    141129      ENDIF 
    142130      ! 
    143       ! ------------------------------- 
    144       !- check conservation (C Rousset) 
    145       IF( ln_limdiahsb ) THEN 
    146          zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
    147          zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 
    148   
    149          zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 
    150          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
    151  
    152          zchk_vmin = glob_min(v_i) 
    153          zchk_amax = glob_max(SUM(a_i,dim=3)) 
    154          zchk_amin = glob_min(a_i) 
    155  
    156          IF(lwp) THEN 
    157             IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_th) = ',(zchk_v_i * rday) 
    158             IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * rday) 
    159             IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_th) = ',(zchk_vmin * 1.e-3) 
    160             IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limitd_th) = ',zchk_amax 
    161             IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_th) = ',zchk_amin 
    162          ENDIF 
    163        ENDIF 
    164       !- check conservation (C Rousset) 
    165       ! ------------------------------- 
     131      ! conservation test 
     132      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    166133      ! 
    167134     IF( nn_timing == 1 )  CALL timing_stop('limitd_th') 
     
    169136   ! 
    170137 
    171    SUBROUTINE lim_itd_th_rem( klbnd, kubnd, ntyp, kt ) 
     138   SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) 
    172139      !!------------------------------------------------------------------ 
    173140      !!                ***  ROUTINE lim_itd_th_rem *** 
     
    182149      INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point 
    183150      INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied 
    184       INTEGER , INTENT (in) ::   ntyp    ! Number of the type used 
    185151      INTEGER , INTENT (in) ::   kt      ! Ocean time step  
    186152      ! 
     
    200166      REAL(wp), POINTER, DIMENSION(:,:,:) ::   hL          ! left boundary for the ITD for each thickness 
    201167      REAL(wp), POINTER, DIMENSION(:,:,:) ::   hR          ! left boundary for the ITD for each thickness 
    202       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zht_i_o     ! old ice thickness 
     168      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zht_i_b     ! old ice thickness 
    203169      REAL(wp), POINTER, DIMENSION(:,:,:) ::   dummy_es 
    204170      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice          ! local increment of ice area and volume 
     
    218184      CALL wrk_alloc( jpi,jpj, zremap_flag )    ! integer 
    219185      CALL wrk_alloc( jpi,jpj,jpl-1, zdonor )   ! integer 
    220       CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_o, dummy_es ) 
     186      CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 
    221187      CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    222188      CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
     
    247213         WRITE(numout,*) ' klbnd :       ', klbnd 
    248214         WRITE(numout,*) ' kubnd :       ', kubnd 
    249          WRITE(numout,*) ' ntyp  :       ', ntyp  
    250215      ENDIF 
    251216 
     
    256221               zindb             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) )     !0 if no ice and 1 if yes 
    257222               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb 
    258                zindb             = 1.0 - MAX( 0.0, SIGN( 1.0, - old_a_i(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 
    259                zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX( old_a_i(ji,jj,jl), epsi10 ) * zindb 
    260                IF( a_i(ji,jj,jl) > epsi06 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl)  
     223               zindb             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 
     224               zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * zindb 
     225               IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)  
    261226            END DO 
    262227         END DO 
     
    302267            ij = nind_j(ji) 
    303268            ! 
    304             IF ( ( zht_i_o(ii,ij,jl) .GT. epsi10 ) .AND. &  
    305                ( zht_i_o(ii,ij,jl+1) .GT. epsi10 ) ) THEN 
     269            zhbnew(ii,ij,jl) = hi_max(jl) 
     270            IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 
    306271               !interpolate between adjacent category growth rates 
    307                zslope = ( zdhice(ii,ij,jl+1)     - zdhice(ii,ij,jl) ) / & 
    308                   ( zht_i_o   (ii,ij,jl+1) - zht_i_o   (ii,ij,jl) ) 
    309                zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + & 
    310                   zslope * ( hi_max(jl) - zht_i_o(ii,ij,jl) ) 
    311             ELSEIF (zht_i_o(ii,ij,jl).gt.epsi10) THEN 
     272               zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 
     273               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 
     274            ELSEIF ( a_i_b(ii,ij,jl) > epsi10) THEN 
    312275               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 
    313             ELSEIF (zht_i_o(ii,ij,jl+1).gt.epsi10) THEN 
     276            ELSEIF ( a_i_b(ii,ij,jl+1) > epsi10) THEN 
    314277               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 
    315             ELSE 
    316                zhbnew(ii,ij,jl) = hi_max(jl) 
    317278            ENDIF 
    318279         END DO 
     
    320281         !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness 
    321282         DO ji = 1, nbrem 
    322             ! jl, ji 
    323283            ii = nind_i(ji) 
    324284            ij = nind_j(ji) 
    325             ! jl, ji 
    326             IF ( ( a_i(ii,ij,jl) .GT.epsi10) .AND. &  
    327                ( ht_i(ii,ij,jl).GE. zhbnew(ii,ij,jl) ) & 
    328                ) THEN 
     285            IF( a_i(ii,ij,jl) > epsi10 .AND. ht_i(ii,ij,jl) >= zhbnew(ii,ij,jl) ) THEN 
    329286               zremap_flag(ii,ij) = 0 
    330             ELSEIF ( ( a_i(ii,ij,jl+1) .GT. epsi10 ) .AND. & 
    331                ( ht_i(ii,ij,jl+1).LE. zhbnew(ii,ij,jl) ) & 
    332                ) THEN 
     287            ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) <= zhbnew(ii,ij,jl) ) THEN 
    333288               zremap_flag(ii,ij) = 0 
    334289            ENDIF 
    335290 
    336291            !- 4.3 Check that each zhbnew does not exceed maximal values hi_max   
    337             ! jl, ji 
    338             IF (zhbnew(ii,ij,jl).gt.hi_max(jl+1)) THEN 
    339                zremap_flag(ii,ij) = 0 
    340             ENDIF 
    341             ! jl, ji 
    342             IF (zhbnew(ii,ij,jl).lt.hi_max(jl-1)) THEN 
    343                zremap_flag(ii,ij) = 0 
    344             ENDIF 
    345             ! jl, ji 
    346          END DO !ji 
    347          ! ji 
     292            IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 
     293            IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 
     294         END DO 
     295 
    348296      END DO !jl 
    349297 
     
    354302      DO jj = 1, jpj 
    355303         DO ji = 1, jpi 
    356             IF ( zremap_flag(ji,jj) == 1 ) THEN 
     304            IF( zremap_flag(ji,jj) == 1 ) THEN 
    357305               nbrem         = nbrem + 1 
    358306               nind_i(nbrem) = ji 
    359307               nind_j(nbrem) = jj 
    360308            ENDIF 
    361          END DO !ji 
    362       END DO !jj 
     309         END DO  
     310      END DO  
    363311 
    364312      !----------------------------------------------------------------------------------------------- 
     
    367315      DO jj = 1, jpj 
    368316         DO ji = 1, jpi 
    369             zhb0(ji,jj) = hi_max_typ(0,ntyp) ! 0eme 
    370             zhb1(ji,jj) = hi_max_typ(1,ntyp) ! 1er 
     317            zhb0(ji,jj) = hi_max(0) ! 0eme 
     318            zhb1(ji,jj) = hi_max(1) ! 1er 
    371319 
    372320            zhbnew(ji,jj,klbnd-1) = 0._wp 
     
    380328            ENDIF 
    381329 
    382             IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) )   zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 
     330            IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 
    383331 
    384332         END DO !jj 
     
    389337      !----------------------------------------------------------------------------------------------- 
    390338      !- 7.1 g(h) for category 1 at start of time step 
    391       CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_o(:,:,klbnd),         & 
     339      CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd),         & 
    392340         &                  g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),   & 
    393341         &                  hR(:,:,klbnd), zremap_flag ) 
     
    414362                  ! Constrain new thickness <= ht_i 
    415363                  zdamax = a_i(ii,ij,klbnd) * &  
    416                      (1.0 - ht_i(ii,ij,klbnd)/zht_i_o(ii,ij,klbnd)) ! zdamax > 0 
     364                     (1.0 - ht_i(ii,ij,klbnd)/zht_i_b(ii,ij,klbnd)) ! zdamax > 0 
    417365                  !ice area lost due to melting of thin ice 
    418366                  zda0   = MIN(zda0, zdamax) 
     
    428376            ELSE ! if ice accretion 
    429377               ! ji, a_i > epsi10; zdh0 > 0 
    430                IF ( ntyp .EQ. 1 ) zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd))  
     378               zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd))  
    431379               ! zhbnew was 0, and is shifted to the right to account for thin ice 
    432380               ! growth in openwater (F0 = f1) 
    433                IF ( ntyp .NE. 1 ) zhbnew(ii,ij,0) = 0  
    434                ! in other types there is 
    435                ! no open water growth (F0 = 0) 
    436381            ENDIF ! zdh0  
    437382 
     
    444389      DO jl = klbnd, kubnd 
    445390         CALL lim_itd_fitline(jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 
    446             g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl),     & 
    447             zremap_flag) 
     391            g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag) 
    448392      END DO 
    449393 
     
    493437            nd   = zdonor(ii,ij,jl) 
    494438            zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 
    495             zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + & 
    496                zdaice(ii,ij,jl)*hL(ii,ij,nd) 
     439            zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 
    497440 
    498441         END DO ! ji 
     
    511454         ii = nind_i(ji) 
    512455         ij = nind_j(ji) 
    513          IF ( ( a_i(ii,ij,1) > epsi10 ) .AND. ( ht_i(ii,ij,1) < hiclim ) ) THEN 
     456         IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < hiclim ) THEN 
    514457            a_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim  
    515458            ht_i(ii,ij,1) = hiclim 
    516             v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem-useless 
    517459         ENDIF 
    518460      END DO !ji 
     
    542484      CALL wrk_dealloc( jpi,jpj, zremap_flag )    ! integer 
    543485      CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor )   ! integer 
    544       CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_o, dummy_es ) 
     486      CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 
    545487      CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    546488      CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
     
    799741            !-------------- 
    800742 
    801             zdvsnow          = v_s(ii,ij,jl1) * zworka(ii,ij) 
     743            zdvsnow        = v_s(ii,ij,jl1) * zworka(ii,ij) 
    802744            v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 
    803745            v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow  
     
    807749            !-------------------- 
    808750 
    809             zdesnow              = e_s(ii,ij,1,jl1) * zworka(ii,ij) 
     751            zdesnow            = e_s(ii,ij,1,jl1) * zworka(ii,ij) 
    810752            e_s(ii,ij,1,jl1)   = e_s(ii,ij,1,jl1) - zdesnow 
    811753            e_s(ii,ij,1,jl2)   = e_s(ii,ij,1,jl2) + zdesnow 
     
    815757            !-------------- 
    816758 
    817             zdo_aice             = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 
     759            zdo_aice           = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 
    818760            oa_i(ii,ij,jl1)    = oa_i(ii,ij,jl1) - zdo_aice 
    819761            oa_i(ii,ij,jl2)    = oa_i(ii,ij,jl2) + zdo_aice 
     
    823765            !-------------- 
    824766 
    825             zdsm_vice            = smv_i(ii,ij,jl1) * zworka(ii,ij) 
     767            zdsm_vice          = smv_i(ii,ij,jl1) * zworka(ii,ij) 
    826768            smv_i(ii,ij,jl1)   = smv_i(ii,ij,jl1) - zdsm_vice 
    827769            smv_i(ii,ij,jl2)   = smv_i(ii,ij,jl2) + zdsm_vice 
     
    831773            !--------------------- 
    832774 
    833             zdaTsf               = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 
     775            zdaTsf             = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 
    834776            zaTsfn(ii,ij,jl1)  = zaTsfn(ii,ij,jl1) - zdaTsf 
    835777            zaTsfn(ii,ij,jl2)  = zaTsfn(ii,ij,jl2) + zdaTsf  
     
    888830    
    889831 
    890    SUBROUTINE lim_itd_th_reb( klbnd, kubnd, ntyp ) 
     832   SUBROUTINE lim_itd_th_reb( klbnd, kubnd ) 
    891833      !!------------------------------------------------------------------ 
    892834      !!                ***  ROUTINE lim_itd_th_reb *** 
     
    898840      INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point 
    899841      INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied 
    900       INTEGER , INTENT (in) ::   ntyp    ! number of the ice type involved in the rebinning process 
    901842      ! 
    902843      INTEGER ::   ji,jj, jl   ! dummy loop indices 
     
    910851      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   ! snow volume summed over categories 
    911852      !!------------------------------------------------------------------ 
     853      !! clem 2014/04: be carefull, rebining does not conserve salt(maybe?) => the difference is taken into account in limupdate 
    912854       
    913855      CALL wrk_alloc( jpi,jpj,jpl, zdonor )   ! interger 
     
    937879 
    938880      !------------------------------------------------------------------------------ 
    939       ! 2) Make sure thickness of cat klbnd is at least hi_max_typ(klbnd) 
     881      ! 2) Make sure thickness of cat klbnd is at least hi_max(klbnd) 
    940882      !------------------------------------------------------------------------------ 
    941883      DO jj = 1, jpj  
    942884         DO ji = 1, jpi  
    943885            IF( a_i(ji,jj,klbnd) > epsi10 ) THEN 
    944                IF( ht_i(ji,jj,klbnd) <= hi_max_typ(0,ntyp) .AND. hi_max_typ(0,ntyp) > 0._wp ) THEN 
    945                   a_i(ji,jj,klbnd)  = v_i(ji,jj,klbnd) / hi_max_typ(0,ntyp)  
    946                   ht_i(ji,jj,klbnd) = hi_max_typ(0,ntyp) 
     886               IF( ht_i(ji,jj,klbnd) <= hi_max(0) .AND. hi_max(0) > 0._wp ) THEN 
     887                  a_i(ji,jj,klbnd)  = v_i(ji,jj,klbnd) / hi_max(0)  
     888                  ht_i(ji,jj,klbnd) = hi_max(0) 
    947889               ENDIF 
    948890            ENDIF 
     
    1015957 
    1016958!clem-change 
     959         DO jj = 1, jpj 
     960            DO ji = 1, jpi 
     961               IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 
     962                  ! 
     963                  zshiftflag = 1 
     964                  zdonor(ji,jj,jl) = jl + 1 
     965                  zdaice(ji,jj,jl) = a_i(ji,jj,jl+1)  
     966                  zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 
     967               ENDIF 
     968            END DO                 ! ji 
     969         END DO                 ! jj 
     970 
     971         IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
     972          
     973         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
     974            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 
     975            ! Reset shift parameters 
     976            zdonor(:,:,jl) = 0 
     977            zdaice(:,:,jl) = 0._wp 
     978            zdvice(:,:,jl) = 0._wp 
     979         ENDIF 
     980!clem-change 
     981 
     982!         ! clem-change begin: why not doing that? 
    1017983!         DO jj = 1, jpj 
    1018984!            DO ji = 1, jpi 
    1019 !               IF( a_i(ji,jj,jl+1) >  epsi10 .AND.   & 
    1020 !                  ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 
    1021 !                  ! 
    1022 !                  zshiftflag = 1 
    1023 !                  zdonor(ji,jj,jl) = jl + 1 
    1024 !                  zdaice(ji,jj,jl) = a_i(ji,jj,jl+1)  
    1025 !                  zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 
     985!               IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 
     986!                  ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10 
     987!                  a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1)  
    1026988!               ENDIF 
    1027989!            END DO                 ! ji 
    1028990!         END DO                 ! jj 
    1029 ! 
    1030 !         IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
    1031 !          
    1032 !         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
    1033 !            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 
    1034 !            ! Reset shift parameters 
    1035 !            zdonor(:,:,jl) = 0 
    1036 !            zdaice(:,:,jl) = 0._wp 
    1037 !            zdvice(:,:,jl) = 0._wp 
    1038 !         ENDIF 
    1039 !clem-change 
    1040  
    1041          ! clem-change begin: why not doing that? 
    1042          DO jj = 1, jpj 
    1043             DO ji = 1, jpi 
    1044                IF( a_i(ji,jj,jl+1) >  epsi10 .AND.   & 
    1045                   ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 
    1046                   ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10 
    1047                   a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1)  
    1048                ENDIF 
    1049             END DO                 ! ji 
    1050          END DO                 ! jj 
    1051991         ! clem-change end 
    1052992 
Note: See TracChangeset for help on using the changeset viewer.