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 5034 for branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 – NEMO

Ignore:
Timestamp:
2015-01-15T14:48:42+01:00 (9 years ago)
Author:
andrewryan
Message:

merge with trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r4333 r5034  
    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 
     
    4445   PUBLIC   lim_itd_fitline 
    4546   PUBLIC   lim_itd_shiftice 
    46  
    47    REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
    48    REAL(wp) ::   epsi06 = 1.e-6_wp   ! 
    4947 
    5048   !!---------------------------------------------------------------------- 
     
    6563      INTEGER, INTENT(in) ::   kt   ! time step index 
    6664      ! 
    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) 
     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  
    7068      !!------------------------------------------------------------------ 
    7169      IF( nn_timing == 1 )  CALL timing_start('limitd_th') 
    7270 
    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       ! ------------------------------- 
     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) 
    8373 
    8474      IF( kt == nit000 .AND. lwp ) THEN 
     
    9383      ! Given thermodynamic growth rates, transport ice between 
    9484      ! 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 
     85      IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt ) 
    10086      ! 
    10187      CALL lim_var_glo2eqv    ! only for info 
     
    10591      !  3) Add frazil ice growing in leads. 
    10692      !------------------------------------------------------------------------------| 
    107  
    10893      CALL lim_thd_lac 
    10994      CALL lim_var_glo2eqv    ! only for info 
    110  
    111      IF(ln_ctl) THEN   ! Control print 
     95      
     96      IF(ln_ctl) THEN   ! Control print 
    11297         CALL prt_ctl_info(' ') 
    11398         CALL prt_ctl_info(' - Cell values : ') 
     
    131116            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : sm_i     : ') 
    132117            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_th  : smv_i    : ') 
    133             DO ja = 1, nlay_i 
     118            DO jk = 1, nlay_i 
    134119               CALL prt_ctl_info(' ') 
    135                CALL prt_ctl_info(' - Layer : ', ivar1=ja) 
     120               CALL prt_ctl_info(' - Layer : ', ivar1=jk) 
    136121               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      : ') 
     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      : ') 
    139124            END DO 
    140125         END DO 
    141126      ENDIF 
    142127      ! 
    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       ! ------------------------------- 
     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) 
    166130      ! 
    167131     IF( nn_timing == 1 )  CALL timing_stop('limitd_th') 
     
    169133   ! 
    170134 
    171    SUBROUTINE lim_itd_th_rem( klbnd, kubnd, ntyp, kt ) 
     135   SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) 
    172136      !!------------------------------------------------------------------ 
    173137      !!                ***  ROUTINE lim_itd_th_rem *** 
     
    182146      INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point 
    183147      INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied 
    184       INTEGER , INTENT (in) ::   ntyp    ! Number of the type used 
    185148      INTEGER , INTENT (in) ::   kt      ! Ocean time step  
    186149      ! 
     
    190153      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars 
    191154      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      - 
    192       REAL(wp) ::   zx3,             zareamin, zindb   !   -      - 
     155      REAL(wp) ::   zx3,             zareamin          !   -      - 
    193156      CHARACTER (len = 15) :: fieldid 
    194157 
     
    200163      REAL(wp), POINTER, DIMENSION(:,:,:) ::   hL          ! left boundary for the ITD for each thickness 
    201164      REAL(wp), POINTER, DIMENSION(:,:,:) ::   hR          ! left boundary for the ITD for each thickness 
    202       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zht_i_o     ! old ice thickness 
     165      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zht_i_b     ! old ice thickness 
    203166      REAL(wp), POINTER, DIMENSION(:,:,:) ::   dummy_es 
    204167      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice          ! local increment of ice area and volume 
     
    218181      CALL wrk_alloc( jpi,jpj, zremap_flag )    ! integer 
    219182      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 ) 
     183      CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 
    221184      CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    222185      CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
     
    247210         WRITE(numout,*) ' klbnd :       ', klbnd 
    248211         WRITE(numout,*) ' kubnd :       ', kubnd 
    249          WRITE(numout,*) ' ntyp  :       ', ntyp  
    250212      ENDIF 
    251213 
     
    254216         DO jj = 1, jpj 
    255217            DO ji = 1, jpi 
    256                zindb             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) )     !0 if no ice and 1 if yes 
    257                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)  
     218               rswitch             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) )     !0 if no ice and 1 if yes 
     219               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 
     221               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)  
    261223            END DO 
    262224         END DO 
     
    302264            ij = nind_j(ji) 
    303265            ! 
    304             IF ( ( zht_i_o(ii,ij,jl) .GT. epsi10 ) .AND. &  
    305                ( zht_i_o(ii,ij,jl+1) .GT. epsi10 ) ) THEN 
     266            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 
    306268               !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 
     269               zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 
     270               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 
    312272               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 
    313             ELSEIF (zht_i_o(ii,ij,jl+1).gt.epsi10) THEN 
     273            ELSEIF ( a_i_b(ii,ij,jl+1) > epsi10) THEN 
    314274               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 
    315             ELSE 
    316                zhbnew(ii,ij,jl) = hi_max(jl) 
    317275            ENDIF 
    318276         END DO 
     
    320278         !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness 
    321279         DO ji = 1, nbrem 
    322             ! jl, ji 
    323280            ii = nind_i(ji) 
    324281            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 
     282            IF( a_i(ii,ij,jl) > epsi10 .AND. ht_i(ii,ij,jl) >= zhbnew(ii,ij,jl) ) THEN 
    329283               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 
     284            ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) <= zhbnew(ii,ij,jl) ) THEN 
    333285               zremap_flag(ii,ij) = 0 
    334286            ENDIF 
    335287 
    336288            !- 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 
     289            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 
    348293      END DO !jl 
    349294 
     
    354299      DO jj = 1, jpj 
    355300         DO ji = 1, jpi 
    356             IF ( zremap_flag(ji,jj) == 1 ) THEN 
     301            IF( zremap_flag(ji,jj) == 1 ) THEN 
    357302               nbrem         = nbrem + 1 
    358303               nind_i(nbrem) = ji 
    359304               nind_j(nbrem) = jj 
    360305            ENDIF 
    361          END DO !ji 
    362       END DO !jj 
     306         END DO  
     307      END DO  
    363308 
    364309      !----------------------------------------------------------------------------------------------- 
     
    367312      DO jj = 1, jpj 
    368313         DO ji = 1, jpi 
    369             zhb0(ji,jj) = hi_max_typ(0,ntyp) ! 0eme 
    370             zhb1(ji,jj) = hi_max_typ(1,ntyp) ! 1er 
     314            zhb0(ji,jj) = hi_max(0) ! 0eme 
     315            zhb1(ji,jj) = hi_max(1) ! 1er 
    371316 
    372317            zhbnew(ji,jj,klbnd-1) = 0._wp 
     
    380325            ENDIF 
    381326 
    382             IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) )   zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 
     327            IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 
    383328 
    384329         END DO !jj 
     
    389334      !----------------------------------------------------------------------------------------------- 
    390335      !- 7.1 g(h) for category 1 at start of time step 
    391       CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_o(:,:,klbnd),         & 
     336      CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd),         & 
    392337         &                  g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),   & 
    393338         &                  hR(:,:,klbnd), zremap_flag ) 
     
    414359                  ! Constrain new thickness <= ht_i 
    415360                  zdamax = a_i(ii,ij,klbnd) * &  
    416                      (1.0 - ht_i(ii,ij,klbnd)/zht_i_o(ii,ij,klbnd)) ! zdamax > 0 
     361                     (1.0 - ht_i(ii,ij,klbnd)/zht_i_b(ii,ij,klbnd)) ! zdamax > 0 
    417362                  !ice area lost due to melting of thin ice 
    418363                  zda0   = MIN(zda0, zdamax) 
     
    428373            ELSE ! if ice accretion 
    429374               ! ji, a_i > epsi10; zdh0 > 0 
    430                IF ( ntyp .EQ. 1 ) zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd))  
     375               zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd))  
    431376               ! zhbnew was 0, and is shifted to the right to account for thin ice 
    432377               ! 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) 
    436378            ENDIF ! zdh0  
    437379 
     
    444386      DO jl = klbnd, kubnd 
    445387         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) 
     388            g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag) 
    448389      END DO 
    449390 
     
    493434            nd   = zdonor(ii,ij,jl) 
    494435            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) 
     436            zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 
    497437 
    498438         END DO ! ji 
     
    511451         ii = nind_i(ji) 
    512452         ij = nind_j(ji) 
    513          IF ( ( a_i(ii,ij,1) > epsi10 ) .AND. ( ht_i(ii,ij,1) < hiclim ) ) THEN 
     453         IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < hiclim ) THEN 
    514454            a_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim  
    515455            ht_i(ii,ij,1) = hiclim 
    516             v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem-useless 
    517456         ENDIF 
    518457      END DO !ji 
     
    542481      CALL wrk_dealloc( jpi,jpj, zremap_flag )    ! integer 
    543482      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 ) 
     483      CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 
    545484      CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    546485      CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
     
    647586      REAL(wp) ::   zdo_aice           ! ice age times volume transferred 
    648587      REAL(wp) ::   zdaTsf             ! aicen*Tsfcn transferred 
    649       REAL(wp) ::   zindsn             ! snow or not 
    650       REAL(wp) ::   zindb              ! ice or not 
    651588 
    652589      INTEGER, POINTER, DIMENSION(:) ::   nind_i, nind_j   ! compressed indices for i/j directions 
     
    775712 
    776713            jl1 = zdonor(ii,ij,jl) 
    777             zindb             = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) ) 
    778             zworka(ii,ij)   = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * zindb 
     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 
    779716            IF( jl1 == jl) THEN   ;   jl2 = jl1+1 
    780717            ELSE                    ;   jl2 = jl  
     
    799736            !-------------- 
    800737 
    801             zdvsnow          = v_s(ii,ij,jl1) * zworka(ii,ij) 
     738            zdvsnow        = v_s(ii,ij,jl1) * zworka(ii,ij) 
    802739            v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 
    803740            v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow  
     
    807744            !-------------------- 
    808745 
    809             zdesnow              = e_s(ii,ij,1,jl1) * zworka(ii,ij) 
     746            zdesnow            = e_s(ii,ij,1,jl1) * zworka(ii,ij) 
    810747            e_s(ii,ij,1,jl1)   = e_s(ii,ij,1,jl1) - zdesnow 
    811748            e_s(ii,ij,1,jl2)   = e_s(ii,ij,1,jl2) + zdesnow 
     
    815752            !-------------- 
    816753 
    817             zdo_aice             = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 
     754            zdo_aice           = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 
    818755            oa_i(ii,ij,jl1)    = oa_i(ii,ij,jl1) - zdo_aice 
    819756            oa_i(ii,ij,jl2)    = oa_i(ii,ij,jl2) + zdo_aice 
     
    823760            !-------------- 
    824761 
    825             zdsm_vice            = smv_i(ii,ij,jl1) * zworka(ii,ij) 
     762            zdsm_vice          = smv_i(ii,ij,jl1) * zworka(ii,ij) 
    826763            smv_i(ii,ij,jl1)   = smv_i(ii,ij,jl1) - zdsm_vice 
    827764            smv_i(ii,ij,jl2)   = smv_i(ii,ij,jl2) + zdsm_vice 
     
    831768            !--------------------- 
    832769 
    833             zdaTsf               = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 
     770            zdaTsf             = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 
    834771            zaTsfn(ii,ij,jl1)  = zaTsfn(ii,ij,jl1) - zdaTsf 
    835772            zaTsfn(ii,ij,jl2)  = zaTsfn(ii,ij,jl2) + zdaTsf  
     
    872809                  ht_i(ji,jj,jl)  =  v_i   (ji,jj,jl) / a_i(ji,jj,jl)  
    873810                  t_su(ji,jj,jl)  =  zaTsfn(ji,jj,jl) / a_i(ji,jj,jl)  
    874                   zindsn          =  1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 
     811                  rswitch         =  1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 
    875812               ELSE 
    876813                  ht_i(ji,jj,jl)  = 0._wp 
     
    888825    
    889826 
    890    SUBROUTINE lim_itd_th_reb( klbnd, kubnd, ntyp ) 
     827   SUBROUTINE lim_itd_th_reb( klbnd, kubnd ) 
    891828      !!------------------------------------------------------------------ 
    892829      !!                ***  ROUTINE lim_itd_th_reb *** 
     
    898835      INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point 
    899836      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 
    901837      ! 
    902838      INTEGER ::   ji,jj, jl   ! dummy loop indices 
     
    910846      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   ! snow volume summed over categories 
    911847      !!------------------------------------------------------------------ 
     848      !! clem 2014/04: be carefull, rebining does not conserve salt(maybe?) => the difference is taken into account in limupdate 
    912849       
    913850      CALL wrk_alloc( jpi,jpj,jpl, zdonor )   ! interger 
     
    937874 
    938875      !------------------------------------------------------------------------------ 
    939       ! 2) Make sure thickness of cat klbnd is at least hi_max_typ(klbnd) 
     876      ! 2) Make sure thickness of cat klbnd is at least hi_max(klbnd) 
    940877      !------------------------------------------------------------------------------ 
    941878      DO jj = 1, jpj  
    942879         DO ji = 1, jpi  
    943880            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) 
     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) 
    947884               ENDIF 
    948885            ENDIF 
     
    1015952 
    1016953!clem-change 
     954         DO jj = 1, jpj 
     955            DO ji = 1, jpi 
     956               IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 
     957                  ! 
     958                  zshiftflag = 1 
     959                  zdonor(ji,jj,jl) = jl + 1 
     960                  zdaice(ji,jj,jl) = a_i(ji,jj,jl+1)  
     961                  zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 
     962               ENDIF 
     963            END DO                 ! ji 
     964         END DO                 ! jj 
     965 
     966         IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
     967          
     968         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
     969            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 
     970            ! Reset shift parameters 
     971            zdonor(:,:,jl) = 0 
     972            zdaice(:,:,jl) = 0._wp 
     973            zdvice(:,:,jl) = 0._wp 
     974         ENDIF 
     975!clem-change 
     976 
     977!         ! clem-change begin: why not doing that? 
    1017978!         DO jj = 1, jpj 
    1018979!            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) 
     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)  
    1026983!               ENDIF 
    1027984!            END DO                 ! ji 
    1028985!         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 
    1051986         ! clem-change end 
    1052987 
Note: See TracChangeset for help on using the changeset viewer.