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 2715 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 – NEMO

Ignore:
Timestamp:
2011-03-30T17:58:35+02:00 (13 years ago)
Author:
rblod
Message:

First attempt to put dynamic allocation on the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r2528 r2715  
    55   !!                   computation of changes in g(h)       
    66   !!====================================================================== 
     7   !! History :   -   !          (W. H. Lipscomb and E.C. Hunke) CICE (c) original code 
     8   !!            3.0  ! 2005-12  (M. Vancoppenolle) adaptation to LIM-3 
     9   !!             -   ! 2006-06  (M. Vancoppenolle) adaptation to include salt, age and types 
     10   !!             -   ! 2007-04  (M. Vancoppenolle) Mass conservation checked 
     11   !!---------------------------------------------------------------------- 
    712#if defined key_lim3 
    813   !!---------------------------------------------------------------------- 
    914   !!   'key_lim3' :                                   LIM3 sea-ice model 
    1015   !!---------------------------------------------------------------------- 
     16   !!   lim_itd_th       : thermodynamics of ice thickness distribution 
     17   !!   lim_itd_th_rem   : 
     18   !!   lim_itd_th_reb   : 
     19   !!   lim_itd_fitline  : 
     20   !!   lim_itd_shiftice : 
    1121   !!---------------------------------------------------------------------- 
    12    USE dom_ice 
     22   USE dom_ice          ! LIM-3 domain 
    1323   USE par_oce          ! ocean parameters 
    14    USE dom_oce 
     24   USE dom_oce          ! ocean domain 
    1525   USE phycst           ! physical constants (ocean directory)  
    16    USE thd_ice 
    17    USE ice 
    18    USE par_ice 
    19    USE limthd_lac 
    20    USE limvar 
    21    USE limcons 
     26   USE thd_ice          ! LIM-3 thermodynamic variables 
     27   USE ice              ! LIM-3 variables 
     28   USE par_ice          ! LIM-3 parameters 
     29   USE limthd_lac       ! LIM-3 lateral accretion 
     30   USE limvar           ! LIM-3 variables 
     31   USE limcons          ! LIM-3 conservation 
    2232   USE prtctl           ! Print control 
    23    USE in_out_manager 
    24    USE lib_mpp  
     33   USE in_out_manager   ! I/O manager 
     34   USE lib_mpp          ! MPP library 
    2535 
    2636   IMPLICIT NONE 
    2737   PRIVATE 
    2838 
    29    PUBLIC lim_itd_th        ! called by ice_stp 
    30    PUBLIC lim_itd_th_rem 
    31    PUBLIC lim_itd_th_reb 
    32    PUBLIC lim_itd_fitline 
    33    PUBLIC lim_itd_shiftice 
    34  
    35    REAL(wp)  ::           &  ! constant values 
    36       epsi20 = 1e-20   ,  & 
    37       epsi13 = 1e-13   ,  & 
    38       zzero  = 0.e0    ,  & 
    39       zone   = 1.e0 
     39   PUBLIC   lim_itd_th        ! called by ice_stp 
     40   PUBLIC   lim_itd_th_rem 
     41   PUBLIC   lim_itd_th_reb 
     42   PUBLIC   lim_itd_fitline 
     43   PUBLIC   lim_itd_shiftice 
     44 
     45   REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values 
     46   REAL(wp) ::   epsi13 = 1e-13_wp   ! 
     47   REAL(wp) ::   epsi10 = 1e-10_wp   ! 
    4048 
    4149   !!---------------------------------------------------------------------- 
    42    !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
     50   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
    4351   !! $Id$ 
    44    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4553   !!---------------------------------------------------------------------- 
    46  
    47  
    4854CONTAINS 
    4955 
     
    5157      !!------------------------------------------------------------------ 
    5258      !!                ***  ROUTINE lim_itd_th *** 
    53       !! ** Purpose : 
    54       !!        This routine computes the thermodynamics of ice thickness 
    55       !!         distribution 
     59      !! 
     60      !! ** Purpose :   computes the thermodynamics of ice thickness distribution 
     61      !! 
    5662      !! ** Method  : 
    57       !! 
    58       !! ** Arguments : 
    59       !!           kideb , kiut : Starting and ending points on which the  
    60       !!                         the computation is applied 
    61       !! 
    62       !! ** Inputs / Ouputs : (global commons) 
    63       !! 
    64       !! ** External :  
    65       !! 
    66       !! ** References : 
    67       !! 
    68       !! ** History : 
    69       !!           (12-2005) Martin Vancoppenolle  
    70       !! 
    71       !!------------------------------------------------------------------ 
    72       !! * Arguments 
    73       INTEGER, INTENT(in) :: kt 
    74       !! * Local variables 
    75       INTEGER ::   jl, ja,   &   ! ice category, layers 
    76          jm,       &   ! ice types    dummy loop index 
    77          jbnd1,    & 
    78          jbnd2 
    79  
    80       REAL(wp)  ::           &  ! constant values 
    81          zeps      =  1.0e-10, & 
    82          epsi10    =  1.0e-10 
     63      !!------------------------------------------------------------------ 
     64      INTEGER, INTENT(in) ::   kt   ! time step index 
     65      ! 
     66      INTEGER ::   jl, ja, jm, jbnd1, jbnd2   ! ice types    dummy loop index          
     67 
     68      !!------------------------------------------------------------------ 
    8369 
    8470      IF( kt == nit000 .AND. lwp ) THEN 
     
    9682         jbnd1 = ice_cat_bounds(jm,1) 
    9783         jbnd2 = ice_cat_bounds(jm,2) 
    98          IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_rem( jbnd1, jbnd2, jm, kt ) 
     84         IF( ice_ncat_types(jm) > 1 )  CALL lim_itd_th_rem( jbnd1, jbnd2, jm, kt ) 
    9985      END DO 
    100  
    101       CALL lim_var_glo2eqv ! only for info 
     86      ! 
     87      CALL lim_var_glo2eqv    ! only for info 
    10288      CALL lim_var_agg(1) 
    10389 
     
    10793 
    10894      CALL lim_thd_lac 
    109       CALL lim_var_glo2eqv ! only for info 
     95      CALL lim_var_glo2eqv    ! only for info 
    11096 
    11197      !---------------------------------------------------------------------------------------- 
     
    120106      d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 
    121107 
    122       d_smv_i_thd(:,:,:) = 0.0 
    123       IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 
    124          d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
     108      d_smv_i_thd(:,:,:) = 0._wp 
     109      IF( num_sal == 2 .OR. num_sal == 4 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
    125110 
    126111      IF(ln_ctl) THEN   ! Control print 
     
    157142 
    158143      !- Recover Old values 
    159       a_i(:,:,:)         = old_a_i (:,:,:) 
    160       v_s(:,:,:)         = old_v_s (:,:,:) 
    161       v_i(:,:,:)         = old_v_i (:,:,:) 
    162       e_s(:,:,:,:)       = old_e_s (:,:,:,:) 
    163       e_i(:,:,:,:)       = old_e_i (:,:,:,:) 
    164  
    165       IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 
    166          smv_i(:,:,:)       = old_smv_i (:,:,:) 
    167  
     144      a_i(:,:,:)   = old_a_i (:,:,:) 
     145      v_s(:,:,:)   = old_v_s (:,:,:) 
     146      v_i(:,:,:)   = old_v_i (:,:,:) 
     147      e_s(:,:,:,:) = old_e_s (:,:,:,:) 
     148      e_i(:,:,:,:) = old_e_i (:,:,:,:) 
     149      ! 
     150      IF( num_sal == 2 .OR. num_sal == 4 )   smv_i(:,:,:)       = old_smv_i (:,:,:) 
     151      ! 
    168152   END SUBROUTINE lim_itd_th 
    169153   ! 
     
    172156      !!------------------------------------------------------------------ 
    173157      !!                ***  ROUTINE lim_itd_th_rem *** 
    174       !! ** Purpose : 
    175       !!        This routine computes the redistribution of ice thickness 
    176       !!        after thermodynamic growth of ice thickness 
     158      !! 
     159      !! ** Purpose :  computes the redistribution of ice thickness 
     160      !!              after thermodynamic growth of ice thickness 
    177161      !! 
    178162      !! ** Method  : Linear remapping  
    179163      !! 
    180       !! ** Arguments : 
    181       !!           klbnd, kubnd : Starting and ending category index on which the  
    182       !!                         the computation is applied 
    183       !! 
    184       !! ** Inputs / Ouputs : (global commons) 
    185       !! 
    186       !! ** External :  
    187       !! 
    188       !! ** References : W.H. Lipscomb, JGR 2001 
    189       !! 
    190       !! ** History : 
    191       !!           largely inspired from CICE (c) W. H. Lipscomb and E.C. Hunke 
    192       !!  
    193       !!           (01-2006) Martin Vancoppenolle, UCL-ASTR, translation from 
    194       !!                     CICE 
    195       !!           (06-2006) Adaptation to include salt, age and types 
    196       !!           (04-2007) Mass conservation checked 
    197       !!------------------------------------------------------------------ 
    198       !! * Arguments 
    199  
    200       INTEGER , INTENT (IN) ::  & 
    201          klbnd ,  &  ! Start thickness category index point 
    202          kubnd ,  &  ! End point on which the  the computation is applied 
    203          ntyp  ,  &  ! Number of the type used 
    204          kt          ! Ocean time step  
    205  
    206       !! * Local variables 
    207       INTEGER ::   ji,       &   ! spatial dummy loop index 
    208          jj,       &   ! spatial dummy loop index 
    209          jl,       &   ! ice category dummy loop index 
    210          zji, zjj, &   ! dummy indices used when changing coordinates 
    211          nd            ! used for thickness categories 
    212  
    213       INTEGER , DIMENSION(jpi,jpj,jpl-1) :: &  
    214          zdonor        ! donor category index 
    215  
    216       REAL(wp)  ::           &   ! constant values 
    217          zeps      =  1.0e-10 
    218  
    219       REAL(wp)  ::           &  ! constant values for ice enthalpy 
    220          zindb     ,         & 
    221          zareamin  ,         &  ! minimum tolerated area in a thickness category 
    222          zwk1, zwk2,         &  ! all the following are dummy arguments 
    223          zx1, zx2, zx3,      &  ! 
    224          zetamin   ,         &  ! minimum value of eta 
    225          zetamax   ,         &  ! maximum value of eta 
    226          zdh0      ,         &  !  
    227          zda0      ,         &  ! 
    228          zdamax    ,         &  ! 
    229          zhimin 
     164      !! References : W.H. Lipscomb, JGR 2001 
     165      !!------------------------------------------------------------------ 
     166      INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point 
     167      INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied 
     168      INTEGER , INTENT (in) ::   ntyp    ! Number of the type used 
     169      INTEGER , INTENT (in) ::   kt      ! Ocean time step  
     170      ! 
     171      INTEGER  ::   ji, jj, jl     ! dummy loop index 
     172      INTEGER  ::   zji, zjj, nd   ! local integer 
     173      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars 
     174      REAL(wp) ::   zx2, zwk2, zda0, zetamax, zhimin   !   -      - 
     175      REAL(wp) ::   zx3,             zareamin, zindb   !   -      - 
     176      CHARACTER (len = 15) :: fieldid 
     177 
     178      INTEGER , DIMENSION(jpi,jpj,jpl-1) ::   zdonor   ! donor category index 
    230179 
    231180      REAL(wp), DIMENSION(jpi,jpj,jpl) :: & 
     
    238187         dummy_es 
    239188 
    240       REAL(wp), DIMENSION(jpi,jpj,jpl-1) :: & 
    241          zdaice           ,  &  ! local increment of ice area  
    242          zdvice                 ! local increment of ice volume 
    243  
    244       REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: & 
    245          zhbnew                 ! new boundaries of ice categories 
    246  
    247       REAL(wp), DIMENSION(jpi,jpj) :: & 
    248          zhb0, zhb1             ! category boundaries for thinnes categories 
    249  
    250       REAL, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 
    251          zvetamin, zvetamax     ! maximum values for etas 
    252  
    253       INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 
    254          nind_i      ,  &  ! compressed indices for i/j directions 
    255          nind_j 
    256  
    257       INTEGER :: & 
    258          nbrem             ! number of cells with ice to transfer 
    259  
    260       LOGICAL, DIMENSION(jpi,jpj) ::   &  !: 
    261          zremap_flag             ! compute remapping or not ???? 
    262  
    263       REAL(wp)  ::           &  ! constant values for ice enthalpy 
    264          zslope                 ! used to compute local thermodynamic "speeds" 
    265  
    266       REAL (wp), DIMENSION(jpi,jpj) :: &  !  
    267          vt_i_init, vt_i_final,   &  !  ice volume summed over categories 
    268          vt_s_init, vt_s_final,   &  !  snow volume summed over categories 
    269          et_i_init, et_i_final,   &  !  ice energy summed over categories 
    270          et_s_init, et_s_final       !  snow energy summed over categories 
    271  
    272       CHARACTER (len = 15) :: fieldid 
    273  
    274       !!-- End of declarations 
    275       !!---------------------------------------------------------------------------------------------- 
    276       zhimin = 0.1      !minimum ice thickness tolerated by the model 
    277       zareamin = zeps   !minimum area in thickness categories tolerated by the conceptors of the model 
     189      REAL(wp), DIMENSION(jpi,jpj,jpl-1) ::   zdaice, zdvice   ! local increment of ice area and volume 
     190 
     191      REAL(wp), DIMENSION(jpi,jpj,0:jpl) ::   zhbnew           ! new boundaries of ice categories 
     192 
     193 
     194      REAL, DIMENSION(1:(jpi+1)*(jpj+1)) ::   zvetamin, zvetamax     ! maximum values for etas 
     195 
     196      INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) ::   nind_i, nind_j  ! compressed indices for i/j directions 
     197 
     198      INTEGER ::   nbrem             ! number of cells with ice to transfer 
     199 
     200      LOGICAL, DIMENSION(jpi,jpj) ::   zremap_flag             ! compute remapping or not ???? 
     201 
     202      REAL(wp)  ::   zslope                 ! used to compute local thermodynamic "speeds" 
     203 
     204      REAL(wp), DIMENSION(jpi,jpj) ::   zhb0, zhb1             ! category boundaries for thinnes categories 
     205      REAL(wp), DIMENSION(jpi,jpj) ::   vt_i_init, vt_i_final   !  ice volume summed over categories 
     206      REAL(wp), DIMENSION(jpi,jpj) ::   vt_s_init, vt_s_final   !  snow volume summed over categories 
     207      REAL(wp), DIMENSION(jpi,jpj) ::   et_i_init, et_i_final   !  ice energy summed over categories 
     208      REAL(wp), DIMENSION(jpi,jpj) ::   et_s_init, et_s_final   !  snow energy summed over categories 
     209      !!------------------------------------------------------------------ 
     210 
     211      zhimin   = 0.1      !minimum ice thickness tolerated by the model 
     212      zareamin = epsi10   !minimum area in thickness categories tolerated by the conceptors of the model 
    278213 
    279214      !!---------------------------------------------------------------------------------------------- 
    280215      !! 0) Conservation checkand changes in each ice category 
    281216      !!---------------------------------------------------------------------------------------------- 
    282       IF ( con_i ) THEN 
     217      IF( con_i ) THEN 
    283218         CALL lim_column_sum (jpl,   v_i, vt_i_init) 
    284219         CALL lim_column_sum (jpl,   v_s, vt_s_init) 
     
    291226      !! 1) Compute thickness and changes in each ice category 
    292227      !!---------------------------------------------------------------------------------------------- 
    293       IF (kt == nit000 .AND. lwp) THEN 
     228      IF( kt == nit000 .AND. lwp) THEN 
    294229         WRITE(numout,*) 
    295230         WRITE(numout,*) 'lim_itd_th_rem  : Remapping the ice thickness distribution' 
     
    300235      ENDIF 
    301236 
    302       zdhice(:,:,:) = 0.0 
     237      zdhice(:,:,:) = 0._wp 
    303238      DO jl = klbnd, kubnd 
    304239         DO jj = 1, jpj 
    305240            DO ji = 1, jpi 
    306241               zindb             = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)))     !0 if no ice and 1 if yes 
    307                ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),zeps) * zindb 
     242               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),epsi10) * zindb 
    308243               zindb             = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl))) !0 if no ice and 1 if yes 
    309                zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),zeps) * zindb 
    310                IF (a_i(ji,jj,jl).gt.1e-6) THEN 
    311                   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl)  
    312                ENDIF 
     244               zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),epsi10) * zindb 
     245               IF( a_i(ji,jj,jl) > 1e-6 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl)  
    313246            END DO 
    314247         END DO 
     
    318251      !  2) Compute fractional ice area in each grid cell 
    319252      !----------------------------------------------------------------------------------------------- 
    320       at_i(:,:) = 0.0 
     253      at_i(:,:) = 0._wp 
    321254      DO jl = klbnd, kubnd 
    322          DO jj = 1, jpj 
    323             DO ji = 1, jpi 
    324                at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) 
    325             END DO 
    326          END DO 
     255         at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    327256      END DO 
    328257 
     
    351280      ! will be soon removed, CT 
    352281      ! hi_max(kubnd) = 999.99 
    353       zhbnew(:,:,:) = 0.0 
     282      zhbnew(:,:,:) = 0._wp 
    354283 
    355284      DO jl = klbnd, kubnd - 1 
    356          ! jl 
    357285         DO ji = 1, nbrem 
    358             ! jl, ji 
    359286            zji = nind_i(ji) 
    360287            zjj = nind_j(ji) 
    361288            ! 
    362             IF ( ( zht_i_o(zji,zjj,jl)  .GT.zeps ) .AND. &  
    363                ( zht_i_o(zji,zjj,jl+1).GT.zeps ) ) THEN 
     289            IF ( ( zht_i_o(zji,zjj,jl)  .GT.epsi10 ) .AND. &  
     290               ( zht_i_o(zji,zjj,jl+1).GT.epsi10 ) ) THEN 
    364291               !interpolate between adjacent category growth rates 
    365292               zslope = ( zdhice(zji,zjj,jl+1)     - zdhice(zji,zjj,jl) ) / & 
     
    367294               zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl) + & 
    368295                  zslope * ( hi_max(jl) - zht_i_o(zji,zjj,jl) ) 
    369             ELSEIF (zht_i_o(zji,zjj,jl).gt.zeps) THEN 
     296            ELSEIF (zht_i_o(zji,zjj,jl).gt.epsi10) THEN 
    370297               zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl) 
    371             ELSEIF (zht_i_o(zji,zjj,jl+1).gt.zeps) THEN 
     298            ELSEIF (zht_i_o(zji,zjj,jl+1).gt.epsi10) THEN 
    372299               zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl+1) 
    373300            ELSE 
    374301               zhbnew(zji,zjj,jl) = hi_max(jl) 
    375302            ENDIF 
    376             ! jl, ji 
    377          END DO !ji 
    378          ! jl 
     303         END DO 
    379304 
    380305         !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness 
     
    384309            zjj = nind_j(ji) 
    385310            ! jl, ji 
    386             IF ( ( a_i(zji,zjj,jl) .GT.zeps) .AND. &  
     311            IF ( ( a_i(zji,zjj,jl) .GT.epsi10) .AND. &  
    387312               ( ht_i(zji,zjj,jl).GE. zhbnew(zji,zjj,jl) ) & 
    388313               ) THEN 
    389314               zremap_flag(zji,zjj) = .false. 
    390             ELSEIF ( ( a_i(zji,zjj,jl+1) .GT. zeps ) .AND. & 
     315            ELSEIF ( ( a_i(zji,zjj,jl+1) .GT. epsi10 ) .AND. & 
    391316               ( ht_i(zji,zjj,jl+1).LE. zhbnew(zji,zjj,jl) ) & 
    392317               ) THEN 
     
    430355            zhb1(ji,jj) = hi_max_typ(1,ntyp) ! 1er 
    431356 
    432             zhbnew(ji,jj,klbnd-1) = 0.0 
    433  
    434             IF ( a_i(ji,jj,kubnd) .GT. zeps ) THEN 
    435                zhbnew(ji,jj,kubnd) = 3.0*ht_i(ji,jj,kubnd) - 2.0*zhbnew(ji,jj,kubnd-1) 
     357            zhbnew(ji,jj,klbnd-1) = 0._wp 
     358 
     359            IF( a_i(ji,jj,kubnd) > epsi10 ) THEN 
     360               zhbnew(ji,jj,kubnd) = 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) 
    436361            ELSE 
    437362               zhbnew(ji,jj,kubnd) = hi_max(kubnd) 
    438363            ENDIF 
    439364 
    440             IF ( zhbnew(ji,jj,kubnd) .LT. hi_max(kubnd-1) ) & 
    441                zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 
     365            IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) )   zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 
    442366 
    443367         END DO !jj 
     
    448372      !----------------------------------------------------------------------------------------------- 
    449373      !- 7.1 g(h) for category 1 at start of time step 
    450       CALL lim_itd_fitline(klbnd, zhb0, zhb1, zht_i_o(:,:,klbnd), & 
    451          g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 
    452          hR(:,:,klbnd), zremap_flag) 
     374      CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_o(:,:,klbnd),        & 
     375         &                  g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),  & 
     376         &                  hR(:,:,klbnd), zremap_flag ) 
    453377 
    454378      !- 7.2 Area lost due to melting of thin ice (first category,  klbnd) 
     
    458382 
    459383         !ji 
    460          IF (a_i(zji,zjj,klbnd) .gt. zeps) THEN 
     384         IF (a_i(zji,zjj,klbnd) .gt. epsi10) THEN 
    461385            zdh0 = zdhice(zji,zjj,klbnd) !decrease of ice thickness in the lower category 
    462             ! ji, a_i > zeps 
     386            ! ji, a_i > epsi10 
    463387            IF (zdh0 .lt. 0.0) THEN !remove area from category 1 
    464                ! ji, a_i > zeps; zdh0 < 0 
     388               ! ji, a_i > epsi10; zdh0 < 0 
    465389               zdh0 = MIN(-zdh0,hi_max(klbnd)) 
    466390 
     
    483407                  v_i(zji,zjj,klbnd)  = a_i(zji,zjj,klbnd)*ht_i(zji,zjj,klbnd) 
    484408               ENDIF     ! zetamax > 0 
    485                ! ji, a_i > zeps 
     409               ! ji, a_i > epsi10 
    486410 
    487411            ELSE ! if ice accretion 
    488                ! ji, a_i > zeps; zdh0 > 0 
     412               ! ji, a_i > epsi10; zdh0 > 0 
    489413               IF ( ntyp .EQ. 1 ) zhbnew(zji,zjj,klbnd-1) = MIN(zdh0,hi_max(klbnd))  
    490414               ! zhbnew was 0, and is shifted to the right to account for thin ice 
     
    495419            ENDIF ! zdh0  
    496420 
    497             ! a_i > zeps 
    498          ENDIF ! a_i > zeps 
     421            ! a_i > epsi10 
     422         ENDIF ! a_i > epsi10 
    499423 
    500424      END DO ! ji 
     
    571495         zjj = nind_j(ji) 
    572496         IF ( ( zhimin .GT. 0.0 ) .AND. &  
    573             ( ( a_i(zji,zjj,1) .GT. zeps ) .AND. ( ht_i(zji,zjj,1) .LT. zhimin ) ) & 
     497            ( ( a_i(zji,zjj,1) .GT. epsi10 ) .AND. ( ht_i(zji,zjj,1) .LT. zhimin ) ) & 
    574498            ) THEN 
    575499            a_i(zji,zjj,1)  = a_i(zji,zjj,1) * ht_i(zji,zjj,1) / zhimin  
     
    602526 
    603527   END SUBROUTINE lim_itd_th_rem 
    604    ! 
    605  
    606    SUBROUTINE lim_itd_fitline(num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 
    607  
     528 
     529 
     530   SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice,   & 
     531      &                        g0, g1, hL, hR, zremap_flag ) 
    608532      !!------------------------------------------------------------------ 
    609533      !!                ***  ROUTINE lim_itd_fitline *** 
    610       !! ** Purpose : 
    611       !! fit g(h) with a line using area, volume constraints 
    612534      !! 
    613       !! ** Method  : 
    614       !! Fit g(h) with a line, satisfying area and volume constraints. 
    615       !! To reduce roundoff errors caused by large values of g0 and g1, 
    616       !! we actually compute g(eta), where eta = h - hL, and hL is the 
    617       !! left boundary. 
     535      !! ** Purpose :   fit g(h) with a line using area, volume constraints 
    618536      !! 
    619       !! ** Arguments : 
    620       !! 
    621       !! ** Inputs / Ouputs : (global commons) 
    622       !! 
    623       !! ** External :  
    624       !! 
    625       !! ** References : 
    626       !! 
    627       !! ** History : 
    628       !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL 
    629       !!          (01-2006) Martin Vancoppenolle  
    630       !! 
    631       !!------------------------------------------------------------------ 
    632       !! * Arguments 
    633  
    634       INTEGER, INTENT(in) :: num_cat      ! category index 
    635  
    636       REAL(wp), DIMENSION(jpi,jpj), INTENT(IN)   ::   &  !: 
    637          HbL, HbR        ! left and right category boundaries 
    638  
    639       REAL(wp), DIMENSION(jpi,jpj), INTENT(IN)   ::   &  !: 
    640          hice            ! ice thickness 
    641  
    642       REAL(wp), DIMENSION(jpi,jpj), INTENT(OUT)  ::   &  !: 
    643          g0, g1      , & ! coefficients in linear equation for g(eta) 
    644          hL          , & ! min value of range over which g(h) > 0 
    645          hR              ! max value of range over which g(h) > 0 
    646  
    647       LOGICAL, DIMENSION(jpi,jpj), INTENT(IN)    ::   &  !: 
    648          zremap_flag 
    649  
    650       INTEGER :: &               
    651          ji,jj           ! horizontal indices 
    652  
    653       REAL(wp) :: &            
    654          zh13        , & ! HbL + 1/3 * (HbR - HbL) 
    655          zh23        , & ! HbL + 2/3 * (HbR - HbL) 
    656          zdhr        , & ! 1 / (hR - hL) 
    657          zwk1, zwk2  , & ! temporary variables 
    658          zacrith         ! critical minimum concentration in an ice category 
    659  
    660       REAL(wp)  ::           &  ! constant values 
    661          zeps      =  1.0e-10 
    662  
     537      !! ** Method  :   Fit g(h) with a line, satisfying area and volume constraints. 
     538      !!              To reduce roundoff errors caused by large values of g0 and g1, 
     539      !!              we actually compute g(eta), where eta = h - hL, and hL is the 
     540      !!              left boundary. 
     541      !!------------------------------------------------------------------ 
     542      INTEGER                     , INTENT(in   ) ::   num_cat      ! category index 
     543      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   HbL, HbR     ! left and right category boundaries 
     544      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   hice         ! ice thickness 
     545      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   g0, g1       ! coefficients in linear equation for g(eta) 
     546      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   hL           ! min value of range over which g(h) > 0 
     547      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   hR           ! max value of range over which g(h) > 0 
     548      LOGICAL , DIMENSION(jpi,jpj), INTENT(in   ) ::   zremap_flag  ! 
     549      ! 
     550      INTEGER ::   ji,jj           ! horizontal indices 
     551      REAL(wp) ::   zh13         ! HbL + 1/3 * (HbR - HbL) 
     552      REAL(wp) ::   zh23         ! HbL + 2/3 * (HbR - HbL) 
     553      REAL(wp) ::   zdhr         ! 1 / (hR - hL) 
     554      REAL(wp) ::   zwk1, zwk2   ! temporary variables 
     555      REAL(wp) ::   zacrith      ! critical minimum concentration in an ice category 
     556      !!------------------------------------------------------------------ 
     557      ! 
    663558      zacrith       = 1.0e-6 
    664       !!-- End of declarations 
    665       !!---------------------------------------------------------------------------------------------- 
    666  
     559      ! 
    667560      DO jj = 1, jpj 
    668561         DO ji = 1, jpi 
    669  
    670             IF ( zremap_flag(ji,jj) .AND. a_i(ji,jj,num_cat) .gt. zacrith & 
    671                .AND. hice(ji,jj) .GT. 0.0 ) THEN 
     562            ! 
     563            IF( zremap_flag(ji,jj) .AND. a_i(ji,jj,num_cat) > zacrith  & 
     564               &                   .AND. hice(ji,jj)        > 0._wp    ) THEN 
    672565 
    673566               ! Initialize hL and hR 
     
    681574               zh23 = 1.0/3.0 * (hL(ji,jj) + 2.0*hR(ji,jj)) 
    682575 
    683                IF (hice(ji,jj) < zh13) THEN 
    684                   hR(ji,jj) = 3.0*hice(ji,jj) - 2.0*hL(ji,jj) 
    685                ELSEIF (hice(ji,jj) > zh23) THEN 
    686                   hL(ji,jj) = 3.0*hice(ji,jj) - 2.0*hR(ji,jj) 
     576               IF    ( hice(ji,jj) < zh13 ) THEN   ;   hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) 
     577               ELSEIF( hice(ji,jj) > zh23 ) THEN   ;   hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj) 
    687578               ENDIF 
    688579 
    689580               ! Compute coefficients of g(eta) = g0 + g1*eta 
    690581 
    691                zdhr = 1.0 / (hR(ji,jj) - hL(ji,jj)) 
    692                zwk1 = 6.0 * a_i(ji,jj,num_cat) * zdhr 
    693                zwk2 = (hice(ji,jj) - hL(ji,jj)) * zdhr 
    694                g0(ji,jj) = zwk1 * (2.0/3.0 - zwk2) 
    695                g1(ji,jj) = 2.0*zdhr * zwk1 * (zwk2 - 0.5) 
    696  
    697             ELSE                   ! remap_flag = .false. or a_i < zeps  
    698  
    699                hL(ji,jj) = 0.0 
    700                hR(ji,jj) = 0.0 
    701                g0(ji,jj) = 0.0 
    702                g1(ji,jj) = 0.0 
    703  
    704             ENDIF                  ! a_i > zeps 
    705  
    706          END DO !ji 
    707       END DO ! jj 
    708  
     582               zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj)) 
     583               zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 
     584               zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr 
     585               g0(ji,jj) = zwk1 * ( 2._wp/3._wp - zwk2 ) 
     586               g1(ji,jj) = 2._wp * zdhr * zwk1 * (zwk2 - 0.5) 
     587               ! 
     588            ELSE                   ! remap_flag = .false. or a_i < epsi10  
     589               hL(ji,jj) = 0._wp 
     590               hR(ji,jj) = 0._wp 
     591               g0(ji,jj) = 0._wp 
     592               g1(ji,jj) = 0._wp 
     593            ENDIF                  ! a_i > epsi10 
     594            ! 
     595         END DO 
     596      END DO 
     597      ! 
    709598   END SUBROUTINE lim_itd_fitline 
    710    ! 
    711  
    712    SUBROUTINE lim_itd_shiftice (klbnd, kubnd, zdonor, zdaice, zdvice) 
     599 
     600 
     601   SUBROUTINE lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 
    713602      !!------------------------------------------------------------------ 
    714603      !!                ***  ROUTINE lim_itd_shiftice *** 
    715       !! ** Purpose : shift ice across category boundaries, conserving everything 
     604      !! 
     605      !! ** Purpose :   shift ice across category boundaries, conserving everything 
    716606      !!              ( area, volume, energy, age*vol, and mass of salt ) 
    717607      !! 
    718608      !! ** Method  : 
    719       !! 
    720       !! ** Arguments : 
    721       !! 
    722       !! ** Inputs / Ouputs : (global commons) 
    723       !! 
    724       !! ** External :  
    725       !! 
    726       !! ** References : 
    727       !! 
    728       !! ** History : 
    729       !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL 
    730       !!          (01-2006) Martin Vancoppenolle  
    731       !! 
    732       !!------------------------------------------------------------------ 
    733       !! * Arguments 
    734  
    735       INTEGER , INTENT (IN) ::  & 
    736          klbnd ,  &  ! Start thickness category index point 
    737          kubnd       ! End point on which the  the computation is applied 
    738  
    739       INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(IN) :: &  
    740          zdonor             ! donor category index 
    741  
    742       REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(INOUT) :: &  
    743          zdaice     ,  &   ! ice area transferred across boundary 
    744          zdvice            ! ice volume transferred across boundary 
    745  
    746       INTEGER :: & 
    747          ji,jj,jl,      &  ! horizontal indices, thickness category index 
    748          jl2,           &  ! receiver category 
    749          jl1,           &  ! donor category 
    750          jk,            &  ! ice layer index 
    751          zji, zjj          ! indices when changing from 2D-1D is done 
    752  
    753       REAL(wp), DIMENSION(jpi,jpj,jpl) :: & 
    754          zaTsfn 
    755  
    756       REAL(wp), DIMENSION(jpi,jpj) :: & 
    757          zworka            ! temporary array used here 
    758  
    759       REAL(wp) :: &           
    760          zdvsnow     ,  &  ! snow volume transferred 
    761          zdesnow     ,  &  ! snow energy transferred 
    762          zdeice      ,  &  ! ice energy transferred 
    763          zdsm_vice      ,  &  ! ice salinity times volume transferred 
    764          zdo_aice      ,  &  ! ice age times volume transferred 
    765          zdaTsf      ,  &  ! aicen*Tsfcn transferred 
    766          zindsn      ,  &  ! snow or not 
    767          zindb             ! ice or not 
    768  
    769       INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 
    770          nind_i      ,  &  ! compressed indices for i/j directions 
    771          nind_j 
    772  
    773       INTEGER :: & 
    774          nbrem             ! number of cells with ice to transfer 
    775  
    776       LOGICAL :: & 
    777          zdaice_negative       , & ! true if daice < -puny 
    778          zdvice_negative       , & ! true if dvice < -puny 
    779          zdaice_greater_aicen  , & ! true if daice > aicen 
    780          zdvice_greater_vicen      ! true if dvice > vicen 
    781  
    782       REAL(wp)  ::           &  ! constant values 
    783          zeps      =  1.0e-10 
    784  
    785       !!-- End of declarations 
     609      !!------------------------------------------------------------------ 
     610      INTEGER , INTENT(in   ) ::   klbnd   ! Start thickness category index point 
     611      INTEGER , INTENT(in   ) ::   kubnd   ! End point on which the  the computation is applied 
     612 
     613      INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(in   ) ::   zdonor   ! donor category index 
     614 
     615      REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) ::   zdaice   ! ice area transferred across boundary 
     616      REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) ::   zdvice   ! ice volume transferred across boundary 
     617 
     618      INTEGER ::   ji, jj, jl, jl2, jl1, jk   ! dummy loop indices 
     619      INTEGER ::   zji, zjj          ! indices when changing from 2D-1D is done 
     620 
     621      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zaTsfn 
     622 
     623      REAL(wp), DIMENSION(jpi,jpj) ::   zworka            ! temporary array used here 
     624 
     625      REAL(wp) ::   zdvsnow, zdesnow   ! snow volume and energy transferred 
     626      REAL(wp) ::   zdeice             ! ice energy transferred 
     627      REAL(wp) ::   zdsm_vice          ! ice salinity times volume transferred 
     628      REAL(wp) ::   zdo_aice           ! ice age times volume transferred 
     629      REAL(wp) ::   zdaTsf             ! aicen*Tsfcn transferred 
     630      REAL(wp) ::   zindsn             ! snow or not 
     631      REAL(wp) ::   zindb              ! ice or not 
     632 
     633      INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) ::   nind_i, nind_j   ! compressed indices for i/j directions 
     634 
     635      INTEGER ::   nbrem             ! number of cells with ice to transfer 
     636 
     637      LOGICAL ::   zdaice_negative         ! true if daice < -puny 
     638      LOGICAL ::   zdvice_negative         ! true if dvice < -puny 
     639      LOGICAL ::   zdaice_greater_aicen    ! true if daice > aicen 
     640      LOGICAL ::   zdvice_greater_vicen    ! true if dvice > vicen 
     641      !!------------------------------------------------------------------ 
    786642 
    787643      !---------------------------------------------------------------------------------------------- 
     
    790646 
    791647      DO jl = klbnd, kubnd 
    792          DO jj = 1, jpj 
    793             DO ji = 1, jpi 
    794                zaTsfn(ji,jj,jl) = a_i(ji,jj,jl)*t_su(ji,jj,jl) 
    795             END DO ! ji 
    796          END DO ! jj 
    797       END DO ! jl 
     648         zaTsfn(:,:,jl) = a_i(:,:,jl)*t_su(:,:,jl) 
     649      END DO 
    798650 
    799651      !---------------------------------------------------------------------------------------------- 
     
    821673 
    822674                  IF (zdaice(ji,jj,jl) .LT. 0.0) THEN 
    823                      IF (zdaice(ji,jj,jl) .GT. -zeps) THEN 
     675                     IF (zdaice(ji,jj,jl) .GT. -epsi10) THEN 
    824676                        IF ( ( jl1.EQ.jl   .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) )           & 
    825677                           .OR.                                      & 
     
    838690 
    839691                  IF (zdvice(ji,jj,jl) .LT. 0.0) THEN 
    840                      IF (zdvice(ji,jj,jl) .GT. -zeps ) THEN 
     692                     IF (zdvice(ji,jj,jl) .GT. -epsi10 ) THEN 
    841693                        IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) )     & 
    842694                           .OR.                                     & 
     
    855707 
    856708                  ! If daice is close to aicen, set daice = aicen. 
    857                   IF (zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - zeps ) THEN 
    858                      IF (zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1)+zeps) THEN 
     709                  IF (zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - epsi10 ) THEN 
     710                     IF (zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1)+epsi10) THEN 
    859711                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 
    860712                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
     
    864716                  ENDIF 
    865717 
    866                   IF (zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1)-zeps) THEN 
    867                      IF (zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1)+zeps) THEN 
     718                  IF (zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1)-epsi10) THEN 
     719                     IF (zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1)+epsi10) THEN 
    868720                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 
    869721                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
     
    900752 
    901753            jl1 = zdonor(zji,zjj,jl) 
    902             zindb             = MAX( 0.0 , SIGN( 1.0 , v_i(zji,zjj,jl1) - zeps ) ) 
    903             zworka(zji,zjj)   = zdvice(zji,zjj,jl) / MAX(v_i(zji,zjj,jl1),zeps) * zindb 
    904             IF (jl1 .eq. jl) THEN 
    905                jl2 = jl1+1 
    906             ELSE                ! n1 = n+1 
    907                jl2 = jl  
     754            zindb             = MAX( 0.0 , SIGN( 1.0 , v_i(zji,zjj,jl1) - epsi10 ) ) 
     755            zworka(zji,zjj)   = zdvice(zji,zjj,jl) / MAX(v_i(zji,zjj,jl1),epsi10) * zindb 
     756            IF( jl1 == jl) THEN   ;   jl2 = jl1+1 
     757            ELSE                    ;   jl2 = jl  
    908758            ENDIF 
    909759 
     
    996846         DO jj = 1, jpj 
    997847            DO ji = 1, jpi  
    998                IF ( a_i(ji,jj,jl) .GT. zeps ) THEN  
    999                   ht_i(ji,jj,jl)  =  v_i(ji,jj,jl) / a_i(ji,jj,jl)  
     848               IF ( a_i(ji,jj,jl) > epsi10 ) THEN  
     849                  ht_i(ji,jj,jl)  =  v_i   (ji,jj,jl) / a_i(ji,jj,jl)  
    1000850                  t_su(ji,jj,jl)  =  zaTsfn(ji,jj,jl) / a_i(ji,jj,jl)  
    1001851                  zindsn          =  1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl))) !0 if no ice and 1 if yes 
    1002852               ELSE 
    1003                   ht_i(ji,jj,jl)  = 0.0 
     853                  ht_i(ji,jj,jl)  = 0._wp 
    1004854                  t_su(ji,jj,jl)  = rtt 
    1005855               ENDIF 
     
    1007857         END DO                 ! jj 
    1008858      END DO                    ! jl 
    1009  
     859      ! 
    1010860   END SUBROUTINE lim_itd_shiftice 
    1011    ! 
    1012  
    1013    SUBROUTINE lim_itd_th_reb(klbnd, kubnd, ntyp) 
     861    
     862 
     863   SUBROUTINE lim_itd_th_reb( klbnd, kubnd, ntyp ) 
    1014864      !!------------------------------------------------------------------ 
    1015865      !!                ***  ROUTINE lim_itd_th_reb *** 
     866      !! 
    1016867      !! ** Purpose : rebin - rebins thicknesses into defined categories 
    1017868      !! 
    1018869      !! ** Method  : 
    1019       !! 
    1020       !! ** Arguments : 
    1021       !! 
    1022       !! ** Inputs / Ouputs : (global commons) 
    1023       !! 
    1024       !! ** External :  
    1025       !! 
    1026       !! ** References : 
    1027       !! 
    1028       !! ** History : (2005) Translation from CICE 
    1029       !!              (2006) Adaptation to include salt, age and types 
    1030       !!              (2007) Mass conservation checked 
    1031       !! 
    1032       !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL 
    1033       !!          (01-2006) Martin Vancoppenolle (adaptation) 
    1034       !! 
    1035       !!------------------------------------------------------------------ 
    1036       !! * Arguments 
    1037       INTEGER , INTENT (in) ::  & 
    1038          klbnd ,  &  ! Start thickness category index point 
    1039          kubnd ,  &  ! End point on which the  the computation is applied 
    1040          ntyp        ! number of the ice type involved in the rebinning process 
    1041  
    1042       INTEGER :: & 
    1043          ji,jj,          &  ! horizontal indices 
    1044          jl                 ! category index 
    1045  
    1046       INTEGER ::   &  !: 
    1047          zshiftflag          ! = .true. if ice must be shifted 
    1048  
    1049       INTEGER, DIMENSION(jpi,jpj,jpl) :: & 
    1050          zdonor             ! donor category index 
    1051  
    1052       REAL(wp), DIMENSION(jpi, jpj, jpl) :: & 
    1053          zdaice         , & ! ice area transferred 
    1054          zdvice             ! ice volume transferred 
    1055  
    1056       REAL(wp)  ::           &  ! constant values 
    1057          zeps      =  1.0e-10, & 
    1058          epsi10    =  1.0e-10 
    1059  
    1060       REAL (wp), DIMENSION(jpi,jpj) :: &  !  
    1061          vt_i_init, vt_i_final,   &  !  ice volume summed over categories 
    1062          vt_s_init, vt_s_final       !  snow volume summed over categories 
    1063  
     870      !!------------------------------------------------------------------ 
     871      INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point 
     872      INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied 
     873      INTEGER , INTENT (in) ::   ntyp    ! number of the ice type involved in the rebinning process 
     874      ! 
     875      INTEGER ::   ji,jj, jl   ! dummy loop indices 
     876      INTEGER ::   zshiftflag          ! = .true. if ice must be shifted 
    1064877      CHARACTER (len = 15) :: fieldid 
    1065878 
    1066       !!-- End of declarations 
    1067       !------------------------------------------------------------------------------ 
    1068  
    1069       !     ! conservation check 
    1070       IF ( con_i ) THEN 
     879      INTEGER , DIMENSION(jpi,jpj,jpl) ::   zdonor           ! donor category index 
     880      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zdaice, zdvice   ! ice area and volume transferred 
     881 
     882      REAL (wp), DIMENSION(jpi,jpj) ::   vt_i_init, vt_i_final   ! ice volume summed over categories 
     883      REAL (wp), DIMENSION(jpi,jpj) ::   vt_s_init, vt_s_final   ! snow volume summed over categories 
     884      !!------------------------------------------------------------------ 
     885      !      
     886      IF( con_i ) THEN                 ! conservation check 
    1071887         CALL lim_column_sum (jpl,   v_i, vt_i_init) 
    1072888         CALL lim_column_sum (jpl,   v_s, vt_s_init) 
     
    1080896         DO jj = 1, jpj 
    1081897            DO ji = 1, jpi  
    1082                IF (a_i(ji,jj,jl) .GT. zeps) THEN  
     898               IF( a_i(ji,jj,jl) > epsi10 ) THEN  
    1083899                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    1084900               ELSE 
    1085                   ht_i(ji,jj,jl) = 0.0 
     901                  ht_i(ji,jj,jl) = 0._wp 
    1086902               ENDIF 
    1087             END DO                 ! i 
    1088          END DO                 ! j 
    1089       END DO                    ! n 
     903            END DO 
     904         END DO 
     905      END DO 
    1090906 
    1091907      !------------------------------------------------------------------------------ 
     
    1094910      DO jj = 1, jpj  
    1095911         DO ji = 1, jpi  
    1096  
    1097             IF (a_i(ji,jj,klbnd) > zeps) THEN 
    1098                IF (ht_i(ji,jj,klbnd) .LE. hi_max_typ(0,ntyp) .AND. hi_max_typ(0,ntyp) .GT. 0.0 ) THEN 
     912            IF( a_i(ji,jj,klbnd) > epsi10 ) THEN 
     913               IF( ht_i(ji,jj,klbnd) <= hi_max_typ(0,ntyp) .AND. hi_max_typ(0,ntyp) > 0._wp ) THEN 
    1099914                  a_i(ji,jj,klbnd)  = v_i(ji,jj,klbnd) / hi_max_typ(0,ntyp)  
    1100915                  ht_i(ji,jj,klbnd) = hi_max_typ(0,ntyp) 
    1101916               ENDIF 
    1102917            ENDIF 
    1103          END DO                    ! i 
    1104       END DO                    ! j 
     918         END DO 
     919      END DO 
    1105920 
    1106921      !------------------------------------------------------------------------------ 
     
    1111926      ! Initialize shift arrays 
    1112927      !------------------------- 
    1113  
    1114928      DO jl = klbnd, kubnd 
    1115          DO jj = 1, jpj  
    1116             DO ji = 1, jpi 
    1117                zdonor(ji,jj,jl) = 0 
    1118                zdaice(ji,jj,jl) = 0.0 
    1119                zdvice(ji,jj,jl) = 0.0 
    1120             END DO 
    1121          END DO 
     929         zdonor(:,:,jl) = 0 
     930         zdaice(:,:,jl) = 0._wp 
     931         zdvice(:,:,jl) = 0._wp 
    1122932      END DO 
    1123933 
     
    1135945         DO jj = 1, jpj  
    1136946            DO ji = 1, jpi  
    1137                IF (a_i(ji,jj,jl) .GT. zeps .AND. ht_i(ji,jj,jl) .GT. hi_max(jl) ) THEN  
     947               IF( a_i(ji,jj,jl) > epsi10 .AND. ht_i(ji,jj,jl) > hi_max(jl) ) THEN  
    1138948                  zshiftflag        = 1 
    1139949                  zdonor(ji,jj,jl)  = jl  
     
    1143953            END DO                 ! ji 
    1144954         END DO                 ! jj 
    1145          IF( lk_mpp ) CALL mpp_max(zshiftflag) 
    1146  
    1147          IF ( zshiftflag == 1 ) THEN 
    1148  
    1149             !------------------------------ 
    1150             ! Shift ice between categories 
    1151             !------------------------------ 
    1152             CALL lim_itd_shiftice (klbnd, kubnd, zdonor, zdaice, zdvice) 
    1153  
    1154             !------------------------ 
     955         IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
     956 
     957         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
     958            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 
    1155959            ! Reset shift parameters 
    1156             !------------------------ 
    1157             DO jj = 1, jpj 
    1158                DO ji = 1, jpi 
    1159                   zdonor(ji,jj,jl) = 0 
    1160                   zdaice(ji,jj,jl) = 0.0 
    1161                   zdvice(ji,jj,jl) = 0.0 
    1162                END DO 
    1163             END DO 
    1164  
    1165          ENDIF                  ! zshiftflag 
    1166  
     960            zdonor(:,:,jl) = 0 
     961            zdaice(:,:,jl) = 0._wp 
     962            zdvice(:,:,jl) = 0._wp 
     963         ENDIF 
     964         ! 
    1167965      END DO                    ! jl 
    1168966 
     
    1180978         DO jj = 1, jpj 
    1181979            DO ji = 1, jpi 
    1182                IF (a_i(ji,jj,jl+1) .GT. zeps .AND. & 
    1183                   ht_i(ji,jj,jl+1) .LE. hi_max(jl)) THEN 
    1184  
     980               IF( a_i(ji,jj,jl+1) >  epsi10 .AND.  & 
     981                  ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 
     982                  ! 
    1185983                  zshiftflag = 1 
    1186984                  zdonor(ji,jj,jl) = jl + 1 
     
    1191989         END DO                 ! jj 
    1192990 
    1193          IF(lk_mpp) CALL mpp_max(zshiftflag) 
    1194          IF (zshiftflag==1) THEN 
    1195  
    1196             !------------------------------ 
    1197             ! Shift ice between categories 
    1198             !------------------------------ 
    1199             CALL lim_itd_shiftice (klbnd, kubnd, zdonor, zdaice, zdvice) 
    1200  
    1201             !------------------------ 
     991         IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
     992          
     993         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
     994            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 
    1202995            ! Reset shift parameters 
    1203             !------------------------ 
    1204             DO jj = 1, jpj  
    1205                DO ji = 1, jpi  
    1206                   zdonor(ji,jj,jl)  = 0 
    1207                   zdaice(ji,jj,jl)  = 0.0 
    1208                   zdvice(ji,jj,jl)  = 0.0 
    1209                END DO 
    1210             END DO 
    1211  
    1212          ENDIF                  ! zshiftflag 
     996            zdonor(:,:,jl) = 0 
     997            zdaice(:,:,jl) = 0._wp 
     998            zdvice(:,:,jl) = 0._wp 
     999         ENDIF 
    12131000 
    12141001      END DO                    ! jl 
     
    12181005      !------------------------------------------------------------------------------ 
    12191006 
    1220       IF ( con_i ) THEN 
     1007      IF( con_i ) THEN 
    12211008         CALL lim_column_sum (jpl,   v_i, vt_i_final) 
    12221009         fieldid = ' v_i : limitd_reb ' 
     
    12271014         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)  
    12281015      ENDIF 
    1229  
     1016      ! 
    12301017   END SUBROUTINE lim_itd_th_reb 
    12311018 
    12321019#else 
    1233    !!====================================================================== 
    1234    !!                       ***  MODULE limitd_th    *** 
    1235    !!                              no sea ice model 
    1236    !!====================================================================== 
     1020   !!---------------------------------------------------------------------- 
     1021   !!   Default option            Dummy module         NO LIM sea-ice model 
     1022   !!---------------------------------------------------------------------- 
    12371023CONTAINS 
    12381024   SUBROUTINE lim_itd_th           ! Empty routines 
     
    12491035   END SUBROUTINE lim_itd_th_reb 
    12501036#endif 
     1037   !!====================================================================== 
    12511038END MODULE limitd_th 
Note: See TracChangeset for help on using the changeset viewer.