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

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

File:
1 edited

Legend:

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

    r5234 r5443  
    2020   USE sbc_oce        ! Surface boundary condition: ocean fields 
    2121   USE ice            ! LIM variables 
    22    USE par_ice        ! LIM parameters 
    2322   USE thd_ice        ! LIM thermodynamics 
    2423   USE in_out_manager ! I/O manager 
     
    3029   PRIVATE 
    3130 
    32    PUBLIC   lim_thd_dh   ! called by lim_thd 
     31   PUBLIC   lim_thd_dh      ! called by lim_thd 
     32   PUBLIC   lim_thd_snwblow ! called in sbcblk/sbccpl and here 
     33 
     34   INTERFACE lim_thd_snwblow 
     35      MODULE PROCEDURE lim_thd_snwblow_1d, lim_thd_snwblow_2d 
     36   END INTERFACE 
    3337 
    3438   !!---------------------------------------------------------------------- 
     
    7074 
    7175      REAL(wp) ::   ztmelts             ! local scalar 
    72       REAL(wp) ::   zdh, zfdum  ! 
     76      REAL(wp) ::   zfdum        
    7377      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment 
    74       REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads 
    75       REAL(wp) ::   zs_snic  ! snow-ice salinity 
     78      REAL(wp) ::   zs_snic      ! snow-ice salinity 
    7679      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity 
    7780      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity 
     
    8790      REAL(wp) ::   zsstK        ! SST in Kelvin 
    8891 
    89       REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness 
    9092      REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow                       (J.m-3) 
    9193      REAL(wp), POINTER, DIMENSION(:) ::   zq_su       ! heat for surface ablation                   (J.m-2) 
    9294      REAL(wp), POINTER, DIMENSION(:) ::   zq_bo       ! heat for bottom ablation                    (J.m-2) 
    93       REAL(wp), POINTER, DIMENSION(:) ::   zq_1cat     ! corrected heat in case 1-cat and hmelt>15cm (J.m-2) 
    9495      REAL(wp), POINTER, DIMENSION(:) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2) 
    95       REAL(wp), POINTER, DIMENSION(:) ::   zf_tt     ! Heat budget to determine melting or freezing(W.m-2) 
    96       INTEGER , POINTER, DIMENSION(:) ::   icount      ! number of layers vanished by melting  
     96      REAL(wp), POINTER, DIMENSION(:) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2) 
    9797 
    9898      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt  
     
    102102      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah 
    103103      REAL(wp), POINTER, DIMENSION(:,:) ::   zh_i      ! ice layer thickness 
     104      INTEGER , POINTER, DIMENSION(:,:) ::   icount    ! number of layers vanished by melting  
    104105 
    105106      REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2) 
    106107      REAL(wp), POINTER, DIMENSION(:) ::   zqh_s       ! total snow heat content (J.m-2) 
    107108      REAL(wp), POINTER, DIMENSION(:) ::   zq_s        ! total snow enthalpy     (J.m-3) 
    108  
    109       ! mass and salt flux (clem) 
    110       REAL(wp) :: zdvres, zswitch_sal 
     109      REAL(wp), POINTER, DIMENSION(:) ::   zsnw        ! distribution of snow after wind blowing 
     110 
     111      REAL(wp) :: zswitch_sal 
    111112 
    112113      ! Heat conservation  
     
    115116      !!------------------------------------------------------------------ 
    116117 
    117       ! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values) 
    118       SELECT CASE( num_sal )                       ! varying salinity or not 
     118      ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values) 
     119      SELECT CASE( nn_icesal )                       ! varying salinity or not 
    119120         CASE( 1, 3, 4 ) ;   zswitch_sal = 0       ! prescribed salinity profile 
    120121         CASE( 2 )       ;   zswitch_sal = 1       ! varying salinity profile 
    121122      END SELECT 
    122123 
    123       CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) 
    124       CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 
    125       CALL wrk_alloc( jpij, nlay_i+1, zdeltah, zh_i ) 
    126       CALL wrk_alloc( jpij, icount ) 
     124      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 
     125      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s, zsnw ) 
     126      CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 
     127      CALL wrk_alloc( jpij, nlay_i, icount ) 
    127128       
    128129      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp 
     
    130131  
    131132      zqprec (:) = 0._wp ; zq_su  (:) = 0._wp ; zq_bo  (:) = 0._wp ; zf_tt  (:) = 0._wp 
    132       zq_1cat(:) = 0._wp ; zq_rema(:) = 0._wp 
    133  
    134       zh_s     (:) = 0._wp        
     133      zq_rema(:) = 0._wp 
     134 
    135135      zdh_s_pre(:) = 0._wp 
    136136      zdh_s_mel(:) = 0._wp 
     
    141141      zh_i      (:,:) = 0._wp        
    142142      zdeltah   (:,:) = 0._wp        
    143       icount    (:)   = 0 
     143      icount    (:,:) = 0 
     144 
     145      ! Initialize enthalpy at nlay_i+1 
     146      DO ji = kideb, kiut 
     147         q_i_1d(ji,nlay_i+1) = 0._wp 
     148      END DO 
    144149 
    145150      ! initialize layer thicknesses and enthalpies 
     
    148153      DO jk = 1, nlay_i 
    149154         DO ji = kideb, kiut 
    150             h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 
     155            h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i 
    151156            qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) 
    152157         ENDDO 
     
    158163      ! 
    159164      DO ji = kideb, kiut 
    160          rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 
    161          ztmelts       = rswitch * rtt + ( 1._wp - rswitch ) * rtt 
    162  
    163165         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)  
    164166         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)  
    165167 
    166          zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts ) ) 
     168         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
    167169         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 
    168170      END DO 
     
    174176      !------------------------------------------------------------------------------! 
    175177      DO ji = kideb, kiut 
    176          IF( t_s_1d(ji,1) > rtt ) THEN !!! Internal melting 
     178         IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting 
    177179            ! Contribution to heat flux to the ocean [W.m-2], < 0   
    178180            hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 
     
    182184            ht_s_1d(ji)   = 0._wp 
    183185            q_s_1d (ji,1) = 0._wp 
    184             t_s_1d (ji,1) = rtt 
     186            t_s_1d (ji,1) = rt0 
    185187         END IF 
    186188      END DO 
     
    190192      !------------------------------------------------------------! 
    191193      ! 
    192       DO ji = kideb, kiut      
    193          zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 
    194       END DO 
    195       ! 
    196194      DO jk = 1, nlay_s 
    197195         DO ji = kideb, kiut 
    198             zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji) 
     196            zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s 
    199197         END DO 
    200198      END DO 
     
    202200      DO jk = 1, nlay_i 
    203201         DO ji = kideb, kiut 
    204             zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 
     202            zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i 
    205203            zqh_i(ji)   = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) 
    206204         END DO 
     
    225223      ! Martin Vancoppenolle, December 2006 
    226224 
     225      zdeltah(:,:) = 0._wp 
     226      CALL lim_thd_snwblow( 1. - at_i_1d, zsnw ) ! snow distribution over ice after wind blowing 
    227227      DO ji = kideb, kiut 
    228228         !----------- 
     
    230230         !----------- 
    231231         ! thickness change 
    232          zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**betas ) / at_i_1d(ji)  
    233          zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn 
    234          ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) 
    235          zqprec   (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )    
     232         zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji) 
     233         ! enthalpy of the precip (>0, J.m-3) 
     234         zqprec   (ji) = - qprec_ice_1d(ji)    
    236235         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 
    237236         ! heat flux from snow precip (>0, W.m-2) 
     
    239238         ! mass flux, <0 
    240239         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice 
    241          ! update thickness 
    242          ht_s_1d    (ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 
    243240 
    244241         !--------------------- 
     
    246243         !--------------------- 
    247244         ! thickness change 
    248          IF( zdh_s_pre(ji) > 0._wp ) THEN 
    249          rswitch        = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) ) 
    250          zdh_s_mel (ji) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 
    251          zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting  
     245         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 
     246         zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 
     247         zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting  
    252248         ! heat used to melt snow (W.m-2, >0) 
    253          hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 
     249         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 
    254250         ! snow melting only = water into the ocean (then without snow precip), >0 
    255          wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_mel(ji) * r1_rdtice 
    256           
    257          ! updates available heat + thickness 
    258          zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) )       
    259          ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) ) 
    260          zh_s  (ji) = ht_s_1d(ji) / REAL( nlay_s ) 
    261  
    262          ENDIF 
    263       END DO 
    264  
    265       ! If heat still available, then melt more snow 
    266       zdeltah(:,:) = 0._wp ! important 
     251         wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice     
     252         ! updates available heat + precipitations after melting 
     253         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )       
     254         zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
     255 
     256         ! update thickness 
     257         ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 
     258      END DO 
     259 
     260      ! If heat still available (zq_su > 0), then melt more snow 
     261      zdeltah(:,:) = 0._wp 
    267262      DO jk = 1, nlay_s 
    268263         DO ji = kideb, kiut 
    269264            ! thickness change 
    270265            rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )  
    271             rswitch          = rswitch * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) + epsi20 ) ) )  
     266            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) )  
    272267            zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) 
    273             zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting 
     268            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting 
    274269            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)     
    275270            ! heat used to melt snow(W.m-2, >0) 
     
    277272            ! snow melting only = water into the ocean (then without snow precip) 
    278273            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
    279  
    280274            ! updates available heat + thickness 
    281             zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 
     275            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 
    282276            ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 
    283  
    284277         END DO 
    285278      END DO 
     
    289282      !---------------------- 
    290283      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 
    291       ! clem comment: not counted in mass exchange in limsbc since this is an exchange with atm. (not ocean) 
     284      ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean) 
    292285      ! clem comment: ice should also sublimate 
    293       IF( lk_cpl ) THEN 
    294          ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) 
    295          zdh_s_sub(:)      =  0._wp  
    296       ELSE 
    297          ! forced  mode: snow thickness change due to sublimation 
    298          DO ji = kideb, kiut 
    299             zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 
    300             ! Heat flux by sublimation [W.m-2], < 0 
    301             !      sublimate first snow that had fallen, then pre-existing snow 
    302             zcoeff         =      ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) )   * zqprec(ji) +   & 
    303                &  ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_1d(ji,1) )  & 
    304                &  * a_i_1d(ji) * r1_rdtice 
    305             hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff 
    306             ! Mass flux by sublimation 
    307             wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 
    308             ! new snow thickness 
    309             ht_s_1d(ji)     =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 
    310          END DO 
    311       ENDIF 
    312  
     286      zdeltah(:,:) = 0._wp 
     287      ! coupled mode: sublimation is set to 0 (evap_ice = 0) until further notice 
     288      ! forced  mode: snow thickness change due to sublimation 
     289      DO ji = kideb, kiut 
     290         zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 
     291         ! Heat flux by sublimation [W.m-2], < 0 
     292         !      sublimate first snow that had fallen, then pre-existing snow 
     293         zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
     294         hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1)  & 
     295            &                              ) * a_i_1d(ji) * r1_rdtice 
     296         ! Mass flux by sublimation 
     297         wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 
     298         ! new snow thickness 
     299         ht_s_1d(ji)    =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 
     300         ! update precipitations after sublimation and correct sublimation 
     301         zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
     302         zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 
     303      END DO 
     304       
    313305      ! --- Update snow diags --- ! 
    314306      DO ji = kideb, kiut 
    315          dh_s_tot(ji)   = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 
    316          zh_s(ji)       = ht_s_1d(ji) / REAL( nlay_s ) 
    317       END DO ! ji 
     307         dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 
     308      END DO 
    318309 
    319310      !------------------------------------------- 
     
    324315      DO jk = 1, nlay_s 
    325316         DO ji = kideb,kiut 
    326             rswitch       =  MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) + epsi20 )  ) 
    327             q_s_1d(ji,jk) = ( 1._wp - rswitch ) / MAX( ht_s_1d(ji), epsi20 ) *             & 
    328               &            ( (   MAX( 0._wp, dh_s_tot(ji) )              ) * zqprec(ji) +  & 
    329               &              ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) ) 
     317            rswitch       = MAX(  0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 )  ) 
     318            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *                          & 
     319              &            ( (   zdh_s_pre(ji)             ) * zqprec(ji) +  & 
     320              &              (   ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
    330321            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk) 
    331322         END DO 
     
    337328      zdeltah(:,:) = 0._wp ! important 
    338329      DO jk = 1, nlay_i 
    339          DO ji = kideb, kiut  
    340             zEi            = - q_i_1d(ji,jk) / rhoic                ! Specific enthalpy of layer k [J/kg, <0] 
    341  
    342             ztmelts        = - tmut * s_i_1d(ji,jk) + rtt           ! Melting point of layer k [K] 
    343  
    344             zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0] 
    345  
    346             zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
    347  
    348             zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0] 
    349  
    350             zdeltah(ji,jk) = - zfmdt / rhoic                       ! Melt of layer jk [m, <0] 
    351  
    352             zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
    353  
    354             zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 
    355  
    356             dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
    357  
    358             zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
    359  
    360             zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
    361  
    362             ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
    363             sfx_sum_1d(ji)   = sfx_sum_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 
    364  
    365             ! Contribution to heat flux [W.m-2], < 0 
    366             hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
    367  
    368             ! Total heat flux used in this process [W.m-2], > 0   
    369             hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
    370  
    371             ! Contribution to mass flux 
    372             wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
    373             
     330         DO ji = kideb, kiut 
     331            ztmelts           = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K] 
     332             
     333            IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 
     334 
     335               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]        
     336               zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0) 
     337                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted) 
     338               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing      
     339                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this 
     340               zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0 
     341                          
     342               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
     343                
     344               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     345 
     346               ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)  
     347               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 
     348                
     349               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
     350               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
     351                
     352               ! Contribution to mass flux 
     353               wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     354 
     355            ELSE                                !!! Surface melting 
     356                
     357               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0] 
     358               zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0] 
     359               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
     360                
     361               zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0] 
     362                
     363               zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0] 
     364                
     365               zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
     366                
     367               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 
     368                
     369               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
     370                
     371               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     372                
     373               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
     374                
     375               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
     376               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
     377                
     378               ! Contribution to heat flux [W.m-2], < 0 
     379               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
     380                
     381               ! Total heat flux used in this process [W.m-2], > 0   
     382               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
     383                
     384               ! Contribution to mass flux 
     385               wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     386                
     387            END IF 
    374388            ! record which layers have disappeared (for bottom melting)  
    375389            !    => icount=0 : no layer has vanished 
    376390            !    => icount=5 : 5 layers have vanished 
    377             rswitch     = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )  
    378             icount(ji)  = icount(ji) + NINT( rswitch ) 
    379             zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 
     391            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )  
     392            icount(ji,jk) = NINT( rswitch ) 
     393            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 
    380394 
    381395            ! update heat content (J.m-2) and layer thickness 
     
    408422      ! -> need for an iterative procedure, which converges quickly 
    409423 
    410       IF ( num_sal == 2 ) THEN 
    411          num_iter_max = 5 
    412       ELSE 
    413          num_iter_max = 1 
    414       ENDIF 
    415  
    416       !clem debug. Just to be sure that enthalpy at nlay_i+1 is null 
    417       DO ji = kideb, kiut 
    418          q_i_1d(ji,nlay_i+1) = 0._wp 
    419       END DO 
     424      num_iter_max = 1 
     425      IF( nn_icesal == 2 ) num_iter_max = 5 
    420426 
    421427      ! Iterative procedure 
     
    440446                                  + ( 1. - zswitch_sal ) * sm_i_1d(ji)  
    441447               ! New ice growth 
    442                ztmelts            = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) 
     448               ztmelts            = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K) 
    443449 
    444450               zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    445451                
    446452               zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)       
    447                   &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   & 
    448                   &               + rcp  * ( ztmelts-rtt )           
     453                  &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   & 
     454                  &               + rcp  * ( ztmelts-rt0 )           
    449455 
    450456               zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     
    456462               q_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0) 
    457463                
    458             ENDIF ! fc_bo_i 
    459          END DO ! ji 
    460       END DO ! iter 
     464            ENDIF 
     465         END DO 
     466      END DO 
    461467 
    462468      ! Contribution to Energy and Salt Fluxes 
     
    467473            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0) 
    468474 
    469             ztmelts        = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) 
     475            ztmelts        = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K) 
    470476             
    471477            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    472478             
    473479            zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)       
    474                &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   & 
    475                &               + rcp  * ( ztmelts-rtt )           
     480               &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   & 
     481               &               + rcp  * ( ztmelts-rt0 )           
    476482             
    477483            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     
    486492             
    487493            ! Contribution to salt flux, <0 
    488             sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_1d(ji) * zfmdt * r1_rdtice 
     494            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice 
    489495 
    490496            ! Contribution to mass flux, <0 
     
    503509      DO jk = nlay_i, 1, -1 
    504510         DO ji = kideb, kiut 
    505             IF(  zf_tt(ji)  >=  0._wp  .AND. jk > icount(ji) ) THEN   ! do not calculate where layer has already disappeared from surface melting  
    506  
    507                ztmelts = - tmut * s_i_1d(ji,jk) + rtt  ! Melting point of layer jk (K) 
     511            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting  
     512 
     513               ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K) 
    508514 
    509515               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 
    510516 
    511                   zEi               = - q_i_1d(ji,jk) / rhoic        ! Specific enthalpy of melting ice (J/kg, <0) 
    512  
    513                   !!zEw               = rcp * ( t_i_1d(ji,jk) - rtt )  ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) 
    514  
     517                  zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0) 
    515518                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
    516519                                                                    ! set up at 0 since no energy is needed to melt water...(it is already melted) 
    517  
    518                   zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing      
    519                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this 
     520                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing      
     521                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this 
    520522 
    521523                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk) 
    522524 
    523                   zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
     525                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0 
    524526 
    525527                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)  
     
    527529 
    528530                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
    529                   sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 
     531                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
    530532                                     
    531533                  ! Contribution to mass flux 
     
    538540               ELSE                               !!! Basal melting 
    539541 
    540                   zEi               = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
    541  
    542                   zEw               = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0) 
    543  
    544                   zdE               = zEi - zEw              ! Specific enthalpy difference   (J/kg, <0) 
    545  
    546                   zfmdt             = - zq_bo(ji) / zdE  ! Mass flux x time step (kg/m2, >0) 
    547  
    548                   zdeltah(ji,jk)    = - zfmdt / rhoic        ! Gross thickness change 
    549  
    550                   zdeltah(ji,jk)    = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 
     542                  zEi             = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
     543                  zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0) 
     544                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0) 
     545 
     546                  zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0) 
     547 
     548                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change 
     549 
     550                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 
    551551                   
    552                   zq_bo(ji)         = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors 
    553  
    554                   dh_i_bott(ji)     = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt 
    555  
    556                   zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
    557  
    558                   zQm               = zfmdt * zEw         ! Heat exchanged with ocean 
     552                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors 
     553 
     554                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt 
     555 
     556                  zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
     557 
     558                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean 
    559559 
    560560                  ! Contribution to heat flux to the ocean [W.m-2], <0   
    561                   hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
     561                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
    562562 
    563563                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
    564                   sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 
     564                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
    565565                   
    566566                  ! Total heat flux used in this process [W.m-2], >0   
    567                   hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
     567                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
    568568                   
    569569                  ! Contribution to mass flux 
    570                   wfx_bom_1d(ji) =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     570                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
    571571 
    572572                  ! update heat content (J.m-2) and layer thickness 
     
    576576            
    577577            ENDIF 
    578          END DO ! ji 
    579       END DO ! jk 
    580  
    581       !------------------------------------------------------------------------------! 
    582       ! Excessive ablation in a 1-category model 
    583       !     in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15) 
    584       !------------------------------------------------------------------------------! 
    585       ! ??? keep ??? 
    586       ! clem bug: I think this should be included above, so we would not have to  
    587       !           track heat/salt/mass fluxes backwards 
    588 !      IF( jpl == 1 ) THEN 
    589 !         DO ji = kideb, kiut 
    590 !            IF(  zf_tt(ji)  >=  0._wp  ) THEN 
    591 !               zdh            = MAX( hmelt , dh_i_bott(ji) ) 
    592 !               zdvres         = zdh - dh_i_bott(ji) ! >=0 
    593 !               dh_i_bott(ji)  = zdh 
    594 ! 
    595 !               ! excessive energy is sent to lateral ablation 
    596 !               rswitch = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_1d(ji) - epsi20 ) ) 
    597 !               zq_1cat(ji) =  rswitch * rhoic * lfus * at_i_1d(ji) / MAX( 1._wp - at_i_1d(ji) , epsi20 ) * zdvres ! J.m-2 >=0 
    598 ! 
    599 !               ! correct salt and mass fluxes 
    600 !               sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation 
    601 !               wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdvres * r1_rdtice 
    602 !            ENDIF 
    603 !         END DO 
    604 !      ENDIF 
     578         END DO 
     579      END DO 
    605580 
    606581      !------------------------------------------- 
     
    619594      DO ji = kideb, kiut 
    620595         zq_rema(ji)     = zq_su(ji) + zq_bo(ji)  
    621 !         zindh           = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow 
    622 !         zindq           = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zq_s(ji) + epsi20 ) ) 
    623 !         zdeltah  (ji,1) = - zindh * zindq * zq_rema(ji) / MAX( zq_s(ji), epsi20 ) 
    624 !         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 
    625 !         zdh_s_mel(ji)   = zdh_s_mel(ji) + zdeltah(ji,1)     
    626 !         dh_s_tot (ji)   = dh_s_tot(ji) + zdeltah(ji,1) 
    627 !         ht_s_1d   (ji)   = ht_s_1d(ji)   + zdeltah(ji,1) 
    628 !         
    629 !         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji)                ! update available heat (J.m-2) 
    630 !         ! heat used to melt snow 
    631 !         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0) 
    632 !         ! Contribution to mass flux 
    633 !         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
    634 !     
     596         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow 
     597         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) ) 
     598         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 ) 
     599         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 
     600         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1) 
     601         ht_s_1d   (ji)  = ht_s_1d(ji)   + zdeltah(ji,1) 
     602         
     603         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1)                ! update available heat (J.m-2) 
     604         ! heat used to melt snow 
     605         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * q_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0) 
     606         ! Contribution to mass flux 
     607         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
     608         !     
    635609         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
    636610         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 
    637          hfx_out(ii,ij)  = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 
    638  
    639          IF( ln_nicep .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 
     611         hfx_out(ii,ij)  = hfx_out(ii,ij) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 
     612 
     613         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 
    640614      END DO 
    641615       
     
    650624         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic )  ) 
    651625 
    652          ht_i_1d(ji)     = ht_i_1d(ji) + dh_snowice(ji) 
    653          ht_s_1d(ji)     = ht_s_1d(ji) - dh_snowice(ji) 
     626         ht_i_1d(ji)    = ht_i_1d(ji) + dh_snowice(ji) 
     627         ht_s_1d(ji)    = ht_s_1d(ji) - dh_snowice(ji) 
    654628 
    655629         ! Salinity of snow ice 
    656630         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
    657          zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
     631         zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
    658632 
    659633         ! entrapment during snow ice formation 
    660          ! new salinity difference stored (to be used in limthd_ent.F90) 
    661          IF (  num_sal == 2  ) THEN 
    662             rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) ) 
     634         ! new salinity difference stored (to be used in limthd_sal.F90) 
     635         IF (  nn_icesal == 2  ) THEN 
     636            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 
    663637            ! salinity dif due to snow-ice formation 
    664             dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi10 ) * rswitch      
     638            dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch      
    665639            ! salinity dif due to bottom growth  
    666640            IF (  zf_tt(ji)  < 0._wp ) THEN 
    667                dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi10 ) * rswitch 
     641               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch 
    668642            ENDIF 
    669643         ENDIF 
     
    691665         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 
    692666          
    693          ! Total ablation (to debug) 
    694          IF( ht_i_1d(ji) <= 0._wp )   a_i_1d(ji) = 0._wp 
    695  
    696       END DO !ji 
     667      END DO 
    697668 
    698669      ! 
     
    700671      ! Update temperature, energy 
    701672      !------------------------------------------- 
    702       !clem bug: we should take snow into account here 
    703673      DO ji = kideb, kiut 
    704674         rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )  
    705          t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rtt 
    706       END DO  ! ji 
     675         t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0 
     676      END DO 
    707677 
    708678      DO jk = 1, nlay_s 
    709679         DO ji = kideb,kiut 
    710680            ! mask enthalpy 
    711             rswitch       = MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  ) 
    712             q_s_1d(ji,jk) = ( 1.0 - rswitch ) * q_s_1d(ji,jk) 
     681            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  ) 
     682            q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk) 
    713683            ! recalculate t_s_1d from q_s_1d 
    714             t_s_1d(ji,jk) = rtt + ( 1._wp - rswitch ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 
     684            t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 
    715685         END DO 
    716686      END DO 
    717687 
    718       CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) 
    719       CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 
    720       CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i ) 
    721       CALL wrk_dealloc( jpij, icount ) 
     688      ! --- ensure that a_i = 0 where ht_i = 0 --- 
     689      WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 
     690       
     691      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 
     692      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s, zsnw ) 
     693      CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 
     694      CALL wrk_dealloc( jpij, nlay_i, icount ) 
    722695      ! 
    723696      ! 
    724697   END SUBROUTINE lim_thd_dh 
     698 
     699 
     700   !!-------------------------------------------------------------------------- 
     701   !! INTERFACE lim_thd_snwblow 
     702   !! ** Purpose :   Compute distribution of precip over the ice 
     703   !!-------------------------------------------------------------------------- 
     704   SUBROUTINE lim_thd_snwblow_2d( pin, pout ) 
     705      REAL(wp), DIMENSION(:,:), INTENT(in)  :: pin   ! previous fraction lead ( pfrld or (1. - a_i_b) ) 
     706      REAL(wp), DIMENSION(:,:), INTENT(out) :: pout 
     707      pout = ( 1._wp - ( pin )**rn_betas ) 
     708   END SUBROUTINE lim_thd_snwblow_2d 
     709 
     710   SUBROUTINE lim_thd_snwblow_1d( pin, pout ) 
     711      REAL(wp), DIMENSION(:), INTENT(in)  :: pin 
     712      REAL(wp), DIMENSION(:), INTENT(out) :: pout 
     713      pout = ( 1._wp - ( pin )**rn_betas ) 
     714   END SUBROUTINE lim_thd_snwblow_1d 
     715 
    725716    
    726717#else 
Note: See TracChangeset for help on using the changeset viewer.