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 5313 for branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 – NEMO

Ignore:
Timestamp:
2015-05-29T11:46:03+02:00 (9 years ago)
Author:
timgraham
Message:

Merged head of trunk (r5302) into branch

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r5312 r5313  
    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 
     
    7069 
    7170      REAL(wp) ::   ztmelts             ! local scalar 
    72       REAL(wp) ::   zdh, zfdum  ! 
     71      REAL(wp) ::   zfdum        
    7372      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment 
    7473      REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads 
     
    8786      REAL(wp) ::   zsstK        ! SST in Kelvin 
    8887 
    89       REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness 
    9088      REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow                       (J.m-3) 
    9189      REAL(wp), POINTER, DIMENSION(:) ::   zq_su       ! heat for surface ablation                   (J.m-2) 
    9290      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) 
    9491      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  
     92      REAL(wp), POINTER, DIMENSION(:) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2) 
    9793 
    9894      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt  
     
    10298      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah 
    10399      REAL(wp), POINTER, DIMENSION(:,:) ::   zh_i      ! ice layer thickness 
     100      INTEGER , POINTER, DIMENSION(:,:) ::   icount    ! number of layers vanished by melting  
    104101 
    105102      REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2) 
     
    107104      REAL(wp), POINTER, DIMENSION(:) ::   zq_s        ! total snow enthalpy     (J.m-3) 
    108105 
    109       ! mass and salt flux (clem) 
    110       REAL(wp) :: zdvres, zswitch_sal 
     106      REAL(wp) :: zswitch_sal 
    111107 
    112108      ! Heat conservation  
     
    115111      !!------------------------------------------------------------------ 
    116112 
    117       ! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values) 
    118       SELECT CASE( num_sal )                       ! varying salinity or not 
     113      ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values) 
     114      SELECT CASE( nn_icesal )                       ! varying salinity or not 
    119115         CASE( 1, 3, 4 ) ;   zswitch_sal = 0       ! prescribed salinity profile 
    120116         CASE( 2 )       ;   zswitch_sal = 1       ! varying salinity profile 
    121117      END SELECT 
    122118 
    123       CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) 
     119      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 
    124120      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 ) 
     121      CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 
     122      CALL wrk_alloc( jpij, nlay_i, icount ) 
    127123       
    128124      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp 
     
    130126  
    131127      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        
     128      zq_rema(:) = 0._wp 
     129 
    135130      zdh_s_pre(:) = 0._wp 
    136131      zdh_s_mel(:) = 0._wp 
     
    141136      zh_i      (:,:) = 0._wp        
    142137      zdeltah   (:,:) = 0._wp        
    143       icount    (:)   = 0 
     138      icount    (:,:) = 0 
     139 
     140      ! Initialize enthalpy at nlay_i+1 
     141      DO ji = kideb, kiut 
     142         q_i_1d(ji,nlay_i+1) = 0._wp 
     143      END DO 
    144144 
    145145      ! initialize layer thicknesses and enthalpies 
     
    148148      DO jk = 1, nlay_i 
    149149         DO ji = kideb, kiut 
    150             h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 
     150            h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i 
    151151            qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) 
    152152         ENDDO 
     
    158158      ! 
    159159      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  
    163160         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)  
    164161         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)  
    165162 
    166          zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts ) ) 
     163         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
    167164         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 
    168165      END DO 
     
    174171      !------------------------------------------------------------------------------! 
    175172      DO ji = kideb, kiut 
    176          IF( t_s_1d(ji,1) > rtt ) THEN !!! Internal melting 
     173         IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting 
    177174            ! Contribution to heat flux to the ocean [W.m-2], < 0   
    178175            hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 
     
    182179            ht_s_1d(ji)   = 0._wp 
    183180            q_s_1d (ji,1) = 0._wp 
    184             t_s_1d (ji,1) = rtt 
     181            t_s_1d (ji,1) = rt0 
    185182         END IF 
    186183      END DO 
     
    190187      !------------------------------------------------------------! 
    191188      ! 
    192       DO ji = kideb, kiut      
    193          zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 
    194       END DO 
    195       ! 
    196189      DO jk = 1, nlay_s 
    197190         DO ji = kideb, kiut 
    198             zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji) 
     191            zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s 
    199192         END DO 
    200193      END DO 
     
    202195      DO jk = 1, nlay_i 
    203196         DO ji = kideb, kiut 
    204             zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 
     197            zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i 
    205198            zqh_i(ji)   = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) 
    206199         END DO 
     
    225218      ! Martin Vancoppenolle, December 2006 
    226219 
     220      zdeltah(:,:) = 0._wp 
    227221      DO ji = kideb, kiut 
    228222         !----------- 
     
    230224         !----------- 
    231225         ! 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 
     226         zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**rn_betas ) / at_i_1d(ji)  
     227         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice * r1_rhosn 
    234228         ! 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 )    
     229         zqprec   (ji) = rhosn * ( cpic * ( rt0 - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )    
    236230         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 
    237231         ! heat flux from snow precip (>0, W.m-2) 
     
    239233         ! mass flux, <0 
    240234         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) ) 
    243235 
    244236         !--------------------- 
     
    246238         !--------------------- 
    247239         ! 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  
     240         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 
     241         zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 
     242         zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting  
    252243         ! 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 
     244         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 
    254245         ! 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 
     246         wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice     
     247         ! updates available heat + precipitations after melting 
     248         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )       
     249         zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
     250 
     251         ! update thickness 
     252         ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 
     253      END DO 
     254 
     255      ! If heat still available (zq_su > 0), then melt more snow 
     256      zdeltah(:,:) = 0._wp 
    267257      DO jk = 1, nlay_s 
    268258         DO ji = kideb, kiut 
    269259            ! thickness change 
    270260            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 ) ) )  
     261            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) )  
    272262            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 
     263            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting 
    274264            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)     
    275265            ! heat used to melt snow(W.m-2, >0) 
     
    277267            ! snow melting only = water into the ocean (then without snow precip) 
    278268            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
    279  
    280269            ! updates available heat + thickness 
    281             zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 
     270            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 
    282271            ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 
    283  
    284272         END DO 
    285273      END DO 
     
    289277      !---------------------- 
    290278      ! 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) 
     279      ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean) 
    292280      ! clem comment: ice should also sublimate 
     281      zdeltah(:,:) = 0._wp 
    293282      IF( lk_cpl ) THEN 
    294283         ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) 
     
    297286         ! forced  mode: snow thickness change due to sublimation 
    298287         DO ji = kideb, kiut 
    299             zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 
     288            zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 
    300289            ! Heat flux by sublimation [W.m-2], < 0 
    301290            !      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 
     291            zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
     292            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)  & 
     293               &                              ) * a_i_1d(ji) * r1_rdtice 
    306294            ! Mass flux by sublimation 
    307295            wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 
    308296            ! new snow thickness 
    309             ht_s_1d(ji)     =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 
     297            ht_s_1d(ji)    =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 
     298            ! update precipitations after sublimation and correct sublimation 
     299            zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
     300            zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 
    310301         END DO 
    311302      ENDIF 
     
    313304      ! --- Update snow diags --- ! 
    314305      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 
     306         dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 
     307      END DO 
    318308 
    319309      !------------------------------------------- 
     
    324314      DO jk = 1, nlay_s 
    325315         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 ) ) 
     316            rswitch       = MAX(  0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 )  ) 
     317            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *                          & 
     318              &            ( (   zdh_s_pre(ji)             ) * zqprec(ji) +  & 
     319              &              (   ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
    330320            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk) 
    331321         END DO 
     
    337327      zdeltah(:,:) = 0._wp ! important 
    338328      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             
     329         DO ji = kideb, kiut 
     330            ztmelts           = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K] 
     331             
     332            IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 
     333 
     334               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]        
     335               zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0) 
     336                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted) 
     337               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing      
     338                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this 
     339               zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0 
     340                          
     341               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
     342                
     343               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     344 
     345               ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)  
     346               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 
     347                
     348               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
     349               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
     350                
     351               ! Contribution to mass flux 
     352               wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     353 
     354            ELSE                                !!! Surface melting 
     355                
     356               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0] 
     357               zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0] 
     358               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
     359                
     360               zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0] 
     361                
     362               zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0] 
     363                
     364               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] 
     365                
     366               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 
     367                
     368               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
     369                
     370               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     371                
     372               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
     373                
     374               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
     375               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
     376                
     377               ! Contribution to heat flux [W.m-2], < 0 
     378               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
     379                
     380               ! Total heat flux used in this process [W.m-2], > 0   
     381               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
     382                
     383               ! Contribution to mass flux 
     384               wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     385                
     386            END IF 
    374387            ! record which layers have disappeared (for bottom melting)  
    375388            !    => icount=0 : no layer has vanished 
    376389            !    => 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) ) 
     390            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )  
     391            icount(ji,jk) = NINT( rswitch ) 
     392            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 
    380393 
    381394            ! update heat content (J.m-2) and layer thickness 
     
    408421      ! -> need for an iterative procedure, which converges quickly 
    409422 
    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 
     423      num_iter_max = 1 
     424      IF( nn_icesal == 2 ) num_iter_max = 5 
    420425 
    421426      ! Iterative procedure 
     
    440445                                  + ( 1. - zswitch_sal ) * sm_i_1d(ji)  
    441446               ! New ice growth 
    442                ztmelts            = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) 
     447               ztmelts            = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K) 
    443448 
    444449               zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    445450                
    446451               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 )           
     452                  &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   & 
     453                  &               + rcp  * ( ztmelts-rt0 )           
    449454 
    450455               zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     
    456461               q_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0) 
    457462                
    458             ENDIF ! fc_bo_i 
    459          END DO ! ji 
    460       END DO ! iter 
     463            ENDIF 
     464         END DO 
     465      END DO 
    461466 
    462467      ! Contribution to Energy and Salt Fluxes 
     
    467472            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0) 
    468473 
    469             ztmelts        = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) 
     474            ztmelts        = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K) 
    470475             
    471476            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    472477             
    473478            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 )           
     479               &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   & 
     480               &               + rcp  * ( ztmelts-rt0 )           
    476481             
    477482            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     
    486491             
    487492            ! 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 
     493            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice 
    489494 
    490495            ! Contribution to mass flux, <0 
     
    503508      DO jk = nlay_i, 1, -1 
    504509         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) 
     510            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting  
     511 
     512               ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K) 
    508513 
    509514               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 
    510515 
    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  
     516                  zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0) 
    515517                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
    516518                                                                    ! 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 
     519                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing      
     520                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this 
    520521 
    521522                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk) 
    522523 
    523                   zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
     524                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0 
    524525 
    525526                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)  
     
    527528 
    528529                  ! 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 
     530                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
    530531                                     
    531532                  ! Contribution to mass flux 
     
    538539               ELSE                               !!! Basal melting 
    539540 
    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 
     541                  zEi             = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
     542                  zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0) 
     543                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0) 
     544 
     545                  zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0) 
     546 
     547                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change 
     548 
     549                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 
    551550                   
    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 
     551                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors 
     552 
     553                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt 
     554 
     555                  zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
     556 
     557                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean 
    559558 
    560559                  ! 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 
     560                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
    562561 
    563562                  ! 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 
     563                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
    565564                   
    566565                  ! 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 
     566                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
    568567                   
    569568                  ! Contribution to mass flux 
    570                   wfx_bom_1d(ji) =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     569                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
    571570 
    572571                  ! update heat content (J.m-2) and layer thickness 
     
    576575            
    577576            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 
     577         END DO 
     578      END DO 
    605579 
    606580      !------------------------------------------- 
     
    619593      DO ji = kideb, kiut 
    620594         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 !     
     595         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow 
     596         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) ) 
     597         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 ) 
     598         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 
     599         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1) 
     600         ht_s_1d   (ji)  = ht_s_1d(ji)   + zdeltah(ji,1) 
     601         
     602         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1)                ! update available heat (J.m-2) 
     603         ! heat used to melt snow 
     604         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) 
     605         ! Contribution to mass flux 
     606         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
     607         !     
    635608         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
    636609         ! 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) 
     610         hfx_out(ii,ij)  = hfx_out(ii,ij) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 
     611 
     612         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 
    640613      END DO 
    641614       
     
    650623         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic )  ) 
    651624 
    652          ht_i_1d(ji)     = ht_i_1d(ji) + dh_snowice(ji) 
    653          ht_s_1d(ji)     = ht_s_1d(ji) - dh_snowice(ji) 
     625         ht_i_1d(ji)    = ht_i_1d(ji) + dh_snowice(ji) 
     626         ht_s_1d(ji)    = ht_s_1d(ji) - dh_snowice(ji) 
    654627 
    655628         ! Salinity of snow ice 
    656629         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) 
     630         zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
    658631 
    659632         ! 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 ) ) 
     633         ! new salinity difference stored (to be used in limthd_sal.F90) 
     634         IF (  nn_icesal == 2  ) THEN 
     635            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 
    663636            ! 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      
     637            dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch      
    665638            ! salinity dif due to bottom growth  
    666639            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 
     640               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch 
    668641            ENDIF 
    669642         ENDIF 
     
    691664         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 
    692665          
    693          ! Total ablation (to debug) 
    694          IF( ht_i_1d(ji) <= 0._wp )   a_i_1d(ji) = 0._wp 
    695  
    696       END DO !ji 
     666      END DO 
    697667 
    698668      ! 
     
    700670      ! Update temperature, energy 
    701671      !------------------------------------------- 
    702       !clem bug: we should take snow into account here 
    703672      DO ji = kideb, kiut 
    704673         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 
     674         t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0 
     675      END DO 
    707676 
    708677      DO jk = 1, nlay_s 
    709678         DO ji = kideb,kiut 
    710679            ! 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) 
     680            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  ) 
     681            q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk) 
    713682            ! 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 ) 
    715          END DO 
    716       END DO 
    717  
    718       CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) 
     683            t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 
     684         END DO 
     685      END DO 
     686 
     687      ! --- ensure that a_i = 0 where ht_i = 0 --- 
     688      WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 
     689       
     690      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 
    719691      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 ) 
     692      CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 
     693      CALL wrk_dealloc( jpij, nlay_i, icount ) 
    722694      ! 
    723695      ! 
Note: See TracChangeset for help on using the changeset viewer.