Ignore:
Timestamp:
2015-08-12T17:46:45+02:00 (6 years ago)
Author:
mattmartin
Message:

OBS simplification changes committed to branch after running SETTE tests to make sure we get the same results as the trunk for ORCA2_LIM_OBS.

File:
1 edited

Legend:

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

    r4990 r5682  
    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/sbcclio/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, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw ) 
    124125      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 ) 
    127        
     126      CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 
     127      CALL wrk_alloc( jpij, nlay_i, icount ) 
     128        
    128129      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp 
    129130      dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp    
    130   
    131       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        
    135       zdh_s_pre(:) = 0._wp 
    136       zdh_s_mel(:) = 0._wp 
    137       zdh_s_sub(:) = 0._wp 
    138       zqh_s    (:) = 0._wp       
    139       zqh_i    (:) = 0._wp    
    140  
    141       zh_i      (:,:) = 0._wp        
    142       zdeltah   (:,:) = 0._wp        
    143       icount    (:)   = 0 
     131 
     132      zqprec   (:) = 0._wp ; zq_su    (:) = 0._wp ; zq_bo    (:) = 0._wp ; zf_tt(:) = 0._wp 
     133      zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp 
     134      zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zqh_i(:) = 0._wp 
     135      zqh_s    (:) = 0._wp ; zq_s     (:) = 0._wp      
     136 
     137      zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp        
     138      icount (:,:) = 0 
     139 
     140 
     141      ! Initialize enthalpy at nlay_i+1 
     142      DO ji = kideb, kiut 
     143         q_i_1d(ji,nlay_i+1) = 0._wp 
     144      END DO 
    144145 
    145146      ! initialize layer thicknesses and enthalpies 
     
    148149      DO jk = 1, nlay_i 
    149150         DO ji = kideb, kiut 
    150             h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 
     151            h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i 
    151152            qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) 
    152153         ENDDO 
     
    158159      ! 
    159160      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  
    163161         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)  
    164162         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)  
    165163 
    166          zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts ) ) 
     164         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
    167165         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 
    168166      END DO 
     
    174172      !------------------------------------------------------------------------------! 
    175173      DO ji = kideb, kiut 
    176          IF( t_s_1d(ji,1) > rtt ) THEN !!! Internal melting 
     174         IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting 
    177175            ! Contribution to heat flux to the ocean [W.m-2], < 0   
    178176            hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 
     
    182180            ht_s_1d(ji)   = 0._wp 
    183181            q_s_1d (ji,1) = 0._wp 
    184             t_s_1d (ji,1) = rtt 
     182            t_s_1d (ji,1) = rt0 
    185183         END IF 
    186184      END DO 
     
    190188      !------------------------------------------------------------! 
    191189      ! 
    192       DO ji = kideb, kiut      
    193          zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 
    194       END DO 
    195       ! 
    196190      DO jk = 1, nlay_s 
    197191         DO ji = kideb, kiut 
    198             zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji) 
     192            zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s 
    199193         END DO 
    200194      END DO 
     
    202196      DO jk = 1, nlay_i 
    203197         DO ji = kideb, kiut 
    204             zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 
     198            zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i 
    205199            zqh_i(ji)   = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) 
    206200         END DO 
     
    225219      ! Martin Vancoppenolle, December 2006 
    226220 
     221      CALL lim_thd_snwblow( 1. - at_i_1d(kideb:kiut), zsnw(kideb:kiut) ) ! snow distribution over ice after wind blowing 
     222 
     223      zdeltah(:,:) = 0._wp 
    227224      DO ji = kideb, kiut 
    228225         !----------- 
     
    230227         !----------- 
    231228         ! 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 )    
     229         zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji) 
     230         ! enthalpy of the precip (>0, J.m-3) 
     231         zqprec   (ji) = - qprec_ice_1d(ji)    
    236232         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 
    237233         ! heat flux from snow precip (>0, W.m-2) 
     
    239235         ! mass flux, <0 
    240236         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) ) 
    243237 
    244238         !--------------------- 
     
    246240         !--------------------- 
    247241         ! 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  
     242         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 
     243         zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 
     244         zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting  
    252245         ! 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 
     246         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 
    254247         ! 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 
     248         wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice     
     249         ! updates available heat + precipitations after melting 
     250         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )       
     251         zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
     252 
     253         ! update thickness 
     254         ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 
     255      END DO 
     256 
     257      ! If heat still available (zq_su > 0), then melt more snow 
     258      zdeltah(:,:) = 0._wp 
    267259      DO jk = 1, nlay_s 
    268260         DO ji = kideb, kiut 
    269261            ! thickness change 
    270262            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 ) ) )  
     263            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) )  
    272264            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 
     265            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting 
    274266            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)     
    275267            ! heat used to melt snow(W.m-2, >0) 
     
    277269            ! snow melting only = water into the ocean (then without snow precip) 
    278270            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
    279  
    280271            ! updates available heat + thickness 
    281             zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 
     272            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 
    282273            ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 
    283  
    284274         END DO 
    285275      END DO 
     
    289279      !---------------------- 
    290280      ! 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) 
     281      ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean) 
    292282      ! 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  
     283      zdeltah(:,:) = 0._wp 
     284      ! coupled mode: sublimation is set to 0 (evap_ice = 0) until further notice 
     285      ! forced  mode: snow thickness change due to sublimation 
     286      DO ji = kideb, kiut 
     287         zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 
     288         ! Heat flux by sublimation [W.m-2], < 0 
     289         !      sublimate first snow that had fallen, then pre-existing snow 
     290         zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
     291         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)  & 
     292            &                              ) * a_i_1d(ji) * r1_rdtice 
     293         ! Mass flux by sublimation 
     294         wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 
     295         ! new snow thickness 
     296         ht_s_1d(ji)    =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 
     297         ! update precipitations after sublimation and correct sublimation 
     298         zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
     299         zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 
     300      END DO 
     301       
    313302      ! --- Update snow diags --- ! 
    314303      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 
     304         dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 
     305      END DO 
    318306 
    319307      !------------------------------------------- 
     
    324312      DO jk = 1, nlay_s 
    325313         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 ) ) 
     314            rswitch       = MAX(  0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 )  ) 
     315            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *                          & 
     316              &            ( (   zdh_s_pre(ji)             ) * zqprec(ji) +  & 
     317              &              (   ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
    330318            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk) 
    331319         END DO 
     
    337325      zdeltah(:,:) = 0._wp ! important 
    338326      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             
     327         DO ji = kideb, kiut 
     328            ztmelts           = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K] 
     329             
     330            IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 
     331 
     332               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]        
     333               zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0) 
     334                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted) 
     335               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing      
     336                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this 
     337               zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0 
     338                          
     339               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
     340                
     341               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     342 
     343               ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)  
     344               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 
     345                
     346               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
     347               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
     348                
     349               ! Contribution to mass flux 
     350               wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     351 
     352            ELSE                                !!! Surface melting 
     353                
     354               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0] 
     355               zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0] 
     356               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
     357                
     358               zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0] 
     359                
     360               zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0] 
     361                
     362               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] 
     363                
     364               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 
     365                
     366               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
     367                
     368               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     369                
     370               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
     371                
     372               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
     373               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
     374                
     375               ! Contribution to heat flux [W.m-2], < 0 
     376               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
     377                
     378               ! Total heat flux used in this process [W.m-2], > 0   
     379               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
     380                
     381               ! Contribution to mass flux 
     382               wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     383                
     384            END IF 
    374385            ! record which layers have disappeared (for bottom melting)  
    375386            !    => icount=0 : no layer has vanished 
    376387            !    => 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) ) 
     388            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )  
     389            icount(ji,jk) = NINT( rswitch ) 
     390            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 
    380391 
    381392            ! update heat content (J.m-2) and layer thickness 
     
    408419      ! -> need for an iterative procedure, which converges quickly 
    409420 
    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 
     421      num_iter_max = 1 
     422      IF( nn_icesal == 2 ) num_iter_max = 5 
    420423 
    421424      ! Iterative procedure 
     
    440443                                  + ( 1. - zswitch_sal ) * sm_i_1d(ji)  
    441444               ! New ice growth 
    442                ztmelts            = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) 
     445               ztmelts            = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K) 
    443446 
    444447               zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    445448                
    446449               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 )           
     450                  &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   & 
     451                  &               + rcp  * ( ztmelts-rt0 )           
    449452 
    450453               zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     
    456459               q_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0) 
    457460                
    458             ENDIF ! fc_bo_i 
    459          END DO ! ji 
    460       END DO ! iter 
     461            ENDIF 
     462         END DO 
     463      END DO 
    461464 
    462465      ! Contribution to Energy and Salt Fluxes 
     
    467470            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0) 
    468471 
    469             ztmelts        = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) 
     472            ztmelts        = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K) 
    470473             
    471474            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    472475             
    473476            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 )           
     477               &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   & 
     478               &               + rcp  * ( ztmelts-rt0 )           
    476479             
    477480            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     
    486489             
    487490            ! 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 
     491            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice 
    489492 
    490493            ! Contribution to mass flux, <0 
     
    503506      DO jk = nlay_i, 1, -1 
    504507         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) 
     508            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting  
     509 
     510               ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K) 
    508511 
    509512               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 
    510513 
    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  
     514                  zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0) 
    515515                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
    516516                                                                    ! 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 
     517                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing      
     518                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this 
    520519 
    521520                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk) 
    522521 
    523                   zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
     522                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0 
    524523 
    525524                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)  
     
    527526 
    528527                  ! 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 
     528                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
    530529                                     
    531530                  ! Contribution to mass flux 
     
    538537               ELSE                               !!! Basal melting 
    539538 
    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 
     539                  zEi             = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
     540                  zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0) 
     541                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0) 
     542 
     543                  zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0) 
     544 
     545                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change 
     546 
     547                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 
    551548                   
    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 
     549                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors 
     550 
     551                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt 
     552 
     553                  zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
     554 
     555                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean 
    559556 
    560557                  ! 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 
     558                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
    562559 
    563560                  ! 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 
     561                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
    565562                   
    566563                  ! 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 
     564                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
    568565                   
    569566                  ! Contribution to mass flux 
    570                   wfx_bom_1d(ji) =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     567                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
    571568 
    572569                  ! update heat content (J.m-2) and layer thickness 
     
    576573            
    577574            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 
     575         END DO 
     576      END DO 
    605577 
    606578      !------------------------------------------- 
     
    619591      DO ji = kideb, kiut 
    620592         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 !     
     593         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow 
     594         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) ) 
     595         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 ) 
     596         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 
     597         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1) 
     598         ht_s_1d   (ji)  = ht_s_1d(ji)   + zdeltah(ji,1) 
     599         
     600         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1)                ! update available heat (J.m-2) 
     601         ! heat used to melt snow 
     602         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) 
     603         ! Contribution to mass flux 
     604         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
     605         !     
    635606         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
    636607         ! 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) 
     608         hfx_out(ii,ij)  = hfx_out(ii,ij) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 
     609 
     610         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 
    640611      END DO 
    641612       
     
    650621         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic )  ) 
    651622 
    652          ht_i_1d(ji)     = ht_i_1d(ji) + dh_snowice(ji) 
    653          ht_s_1d(ji)     = ht_s_1d(ji) - dh_snowice(ji) 
     623         ht_i_1d(ji)    = ht_i_1d(ji) + dh_snowice(ji) 
     624         ht_s_1d(ji)    = ht_s_1d(ji) - dh_snowice(ji) 
    654625 
    655626         ! Salinity of snow ice 
    656627         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) 
     628         zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
    658629 
    659630         ! 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 ) ) 
     631         ! new salinity difference stored (to be used in limthd_sal.F90) 
     632         IF (  nn_icesal == 2  ) THEN 
     633            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 
    663634            ! 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      
     635            dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch      
    665636            ! salinity dif due to bottom growth  
    666637            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 
     638               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch 
    668639            ENDIF 
    669640         ENDIF 
     
    691662         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 
    692663          
    693          ! Total ablation (to debug) 
    694          IF( ht_i_1d(ji) <= 0._wp )   a_i_1d(ji) = 0._wp 
    695  
    696       END DO !ji 
     664      END DO 
    697665 
    698666      ! 
     
    700668      ! Update temperature, energy 
    701669      !------------------------------------------- 
    702       !clem bug: we should take snow into account here 
    703670      DO ji = kideb, kiut 
    704671         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 
     672         t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0 
     673      END DO 
    707674 
    708675      DO jk = 1, nlay_s 
    709676         DO ji = kideb,kiut 
    710677            ! 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) 
     678            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  ) 
     679            q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk) 
    713680            ! 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 ) 
     681            t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 
    715682         END DO 
    716683      END DO 
    717684 
    718       CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) 
     685      ! --- ensure that a_i = 0 where ht_i = 0 --- 
     686      WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 
     687       
     688      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw ) 
    719689      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 ) 
     690      CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 
     691      CALL wrk_dealloc( jpij, nlay_i, icount ) 
    722692      ! 
    723693      ! 
    724694   END SUBROUTINE lim_thd_dh 
     695 
     696 
     697   !!-------------------------------------------------------------------------- 
     698   !! INTERFACE lim_thd_snwblow 
     699   !! ** Purpose :   Compute distribution of precip over the ice 
     700   !!-------------------------------------------------------------------------- 
     701   SUBROUTINE lim_thd_snwblow_2d( pin, pout ) 
     702      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( pfrld or (1. - a_i_b) ) 
     703      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout 
     704      pout = ( 1._wp - ( pin )**rn_betas ) 
     705   END SUBROUTINE lim_thd_snwblow_2d 
     706 
     707   SUBROUTINE lim_thd_snwblow_1d( pin, pout ) 
     708      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin 
     709      REAL(wp), DIMENSION(:), INTENT(inout) :: pout 
     710      pout = ( 1._wp - ( pin )**rn_betas ) 
     711   END SUBROUTINE lim_thd_snwblow_1d 
     712 
    725713    
    726714#else 
Note: See TracChangeset for help on using the changeset viewer.