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

Ignore:
Timestamp:
2015-02-17T10:06:39+01:00 (9 years ago)
Author:
timgraham
Message:

Merged head of trunk into branch in preparation for putting code back onto the trunk
In working copy ran the command:
svn merge svn+sshtimgraham@…/ipsl/forge/projets/nemo/svn/trunk

Also recompiled NEMO_book.pdf with merged input files

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r4333 r5086  
    66   !! History :  LIM  ! 2003-05 (M. Vancoppenolle) Original code in 1D 
    77   !!                 ! 2005-06 (M. Vancoppenolle) 3D version  
    8    !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm_snw & rdm_ice 
     8   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw & wfx_ice 
    99   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation 
    1010   !!            3.5  ! 2012-10 (G. Madec & co) salt flux + bug fixes  
    … …  
    2626   USE wrk_nemo       ! work arrays 
    2727   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    28  
     28    
    2929   IMPLICIT NONE 
    3030   PRIVATE 
    3131 
    3232   PUBLIC   lim_thd_dh   ! called by lim_thd 
    33  
    34    REAL(wp) ::   epsi20 = 1.e-20   ! constant values 
    35    REAL(wp) ::   epsi10 = 1.e-10   ! 
    36    REAL(wp) ::   epsi13 = 1.e-13   ! 
    37    REAL(wp) ::   zzero  = 0._wp    ! 
    38    REAL(wp) ::   zone   = 1._wp    ! 
    3933 
    4034   !!---------------------------------------------------------------------- 
    … …  
    4539CONTAINS 
    4640 
    47    SUBROUTINE lim_thd_dh( kideb, kiut, jl ) 
     41   SUBROUTINE lim_thd_dh( kideb, kiut ) 
    4842      !!------------------------------------------------------------------ 
    4943      !!                ***  ROUTINE lim_thd_dh  *** 
    … …  
    7064      !!------------------------------------------------------------------ 
    7165      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied 
    72       INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number 
    7366      !!  
    7467      INTEGER  ::   ji , jk        ! dummy loop indices 
    7568      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji 
    76       INTEGER  ::   isnow          ! switch for presence (1) or absence (0) of snow 
    77       INTEGER  ::   isnowic        ! snow ice formation not 
    78       INTEGER  ::   i_ice_switch   ! ice thickness above a certain treshold or not 
    7969      INTEGER  ::   iter 
    8070 
    81       REAL(wp) ::   zzfmass_i, zihgnew                     ! local scalar 
    82       REAL(wp) ::   zzfmass_s, zhsnew, ztmelts             ! local scalar 
    83       REAL(wp) ::   zhn, zdhcf, zdhbf, zhni, zhnfi, zihg   ! 
    84       REAL(wp) ::   zdhnm, zhnnew, zhisn, zihic, zzc       ! 
     71      REAL(wp) ::   ztmelts             ! local scalar 
     72      REAL(wp) ::   zdh, zfdum  ! 
    8573      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment 
    8674      REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads 
    87       REAL(wp) ::   zsm_snowice  ! snow-ice salinity 
     75      REAL(wp) ::   zs_snic  ! snow-ice salinity 
    8876      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity 
    8977      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity 
    9078      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity 
    9179      REAL(wp) ::   zgrr         ! bottom growth rate 
    92       REAL(wp) ::   ztform       ! bottom formation temperature 
    93       ! 
    94       REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness 
     80      REAL(wp) ::   zt_i_new     ! bottom formation temperature 
     81 
     82      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean 
     83      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg) 
     84      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg) 
     85      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg) 
     86      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean 
     87      REAL(wp) ::   zsstK        ! SST in Kelvin 
     88 
    9589      REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness 
    96       REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! melting point 
    97       REAL(wp), POINTER, DIMENSION(:) ::   zhsold      ! old snow thickness 
    98       REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow 
    99       REAL(wp), POINTER, DIMENSION(:) ::   zqfont_su   ! incoming, remaining surface energy 
    100       REAL(wp), POINTER, DIMENSION(:) ::   zqfont_bo   ! incoming, bottom energy 
    101       REAL(wp), POINTER, DIMENSION(:) ::   z_f_surf    ! surface heat for ablation 
    102       REAL(wp), POINTER, DIMENSION(:) ::   zhgnew      ! new ice thickness 
    103       REAL(wp), POINTER, DIMENSION(:) ::   zfmass_i    !  
     90      REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow                       (J.m-3) 
     91      REAL(wp), POINTER, DIMENSION(:) ::   zq_su       ! heat for surface ablation                   (J.m-2) 
     92      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) 
     94      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  
    10497 
    10598      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt  
    … …  
    108101 
    109102      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah 
    110  
    111       ! Pathological cases 
    112       REAL(wp), POINTER, DIMENSION(:) ::   zfdt_init   ! total incoming heat for ice melt 
    113       REAL(wp), POINTER, DIMENSION(:) ::   zfdt_final  ! total remaing heat for ice melt 
    114       REAL(wp), POINTER, DIMENSION(:) ::   zqt_i       ! total ice heat content 
    115       REAL(wp), POINTER, DIMENSION(:) ::   zqt_s       ! total snow heat content 
    116       REAL(wp), POINTER, DIMENSION(:) ::   zqt_dummy   ! dummy heat content 
    117  
    118       REAL(wp), POINTER, DIMENSION(:,:) ::   zqt_i_lay   ! total ice heat content 
     103      REAL(wp), POINTER, DIMENSION(:,:) ::   zh_i      ! ice layer thickness 
     104 
     105      REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2) 
     106      REAL(wp), POINTER, DIMENSION(:) ::   zqh_s       ! total snow heat content (J.m-2) 
     107      REAL(wp), POINTER, DIMENSION(:) ::   zq_s        ! total snow enthalpy     (J.m-3) 
    119108 
    120109      ! mass and salt flux (clem) 
    121       REAL(wp) :: zdvres, zdvsur, zdvbot 
    122       REAL(wp), POINTER, DIMENSION(:) ::   zviold, zvsold   ! old ice volume... 
     110      REAL(wp) :: zdvres, zswitch_sal 
    123111 
    124112      ! Heat conservation  
    125       INTEGER  ::   num_iter_max, numce_dh 
    126       REAL(wp) ::   meance_dh 
    127       REAL(wp) ::   zinda  
    128       REAL(wp), POINTER, DIMENSION(:) ::   zinnermelt 
    129       REAL(wp), POINTER, DIMENSION(:) ::   zfbase, zdq_i 
     113      INTEGER  ::   num_iter_max 
     114 
    130115      !!------------------------------------------------------------------ 
    131116 
    132       CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 
    133       CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 
    134       CALL wrk_alloc( jpij, zinnermelt, zfbase, zdq_i ) 
    135       CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay ) 
    136  
    137       CALL wrk_alloc( jpij, zviold, zvsold ) ! clem 
     117      ! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values) 
     118      SELECT CASE( num_sal )                       ! varying salinity or not 
     119         CASE( 1, 3, 4 ) ;   zswitch_sal = 0       ! prescribed salinity profile 
     120         CASE( 2 )       ;   zswitch_sal = 1       ! varying salinity profile 
     121      END SELECT 
     122 
     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 ) 
    138127       
    139       ftotal_fin(:) = 0._wp 
    140       zfdt_init (:) = 0._wp 
    141       zfdt_final(:) = 0._wp 
    142  
    143       dh_i_surf (:) = 0._wp 
    144       dh_i_bott (:) = 0._wp 
    145       dh_snowice(:) = 0._wp 
    146  
    147       DO ji = kideb, kiut 
    148          old_ht_i_b(ji) = ht_i_b(ji) 
    149          old_ht_s_b(ji) = ht_s_b(ji) 
    150          zviold(ji) = a_i_b(ji) * ht_i_b(ji) ! clem 
    151          zvsold(ji) = a_i_b(ji) * ht_s_b(ji) ! clem 
    152       END DO 
     128      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp 
     129      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 
     144 
     145      ! initialize layer thicknesses and enthalpies 
     146      h_i_old (:,0:nlay_i+1) = 0._wp 
     147      qh_i_old(:,0:nlay_i+1) = 0._wp 
     148      DO jk = 1, nlay_i 
     149         DO ji = kideb, kiut 
     150            h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 
     151            qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) 
     152         ENDDO 
     153      ENDDO 
    153154      ! 
    154155      !------------------------------------------------------------------------------! 
    155       !  1) Calculate available heat for surface ablation                            ! 
     156      !  1) Calculate available heat for surface and bottom ablation                 ! 
    156157      !------------------------------------------------------------------------------! 
    157158      ! 
    158159      DO ji = kideb, kiut 
    159          isnow         = INT(  1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_s_b(ji) )  )  ) 
    160          ztfs     (ji) = isnow * rtt + ( 1.0 - isnow ) * rtt 
    161          z_f_surf (ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
    162          z_f_surf (ji) = MAX(  zzero , z_f_surf(ji)  ) * MAX(  zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) )  ) 
    163          zfdt_init(ji) = ( z_f_surf(ji) + MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) * rdt_ice 
    164       END DO ! ji 
    165  
    166       zqfont_su  (:) = 0._wp 
    167       zqfont_bo  (:) = 0._wp 
    168       dsm_i_se_1d(:) = 0._wp      
    169       dsm_i_si_1d(:) = 0._wp    
     160         rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 
     161         ztmelts       = rswitch * rtt + ( 1._wp - rswitch ) * rtt 
     162 
     163         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)  
     164         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)  
     165 
     166         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts ) ) 
     167         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 
     168      END DO 
     169 
    170170      ! 
    171171      !------------------------------------------------------------------------------! 
    172       !  2) Computing layer thicknesses and  snow and sea-ice enthalpies.            ! 
     172      ! If snow temperature is above freezing point, then snow melts  
     173      ! (should not happen but sometimes it does) 
    173174      !------------------------------------------------------------------------------! 
    174       ! 
    175       DO ji = kideb, kiut     ! Layer thickness 
    176          zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 
    177          zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 
    178       END DO 
    179       ! 
    180       zqt_s(:) = 0._wp        ! Total enthalpy of the snow 
     175      DO ji = kideb, kiut 
     176         IF( t_s_1d(ji,1) > rtt ) THEN !!! Internal melting 
     177            ! Contribution to heat flux to the ocean [W.m-2], < 0   
     178            hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 
     179            ! Contribution to mass flux 
     180            wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 
     181            ! updates 
     182            ht_s_1d(ji)   = 0._wp 
     183            q_s_1d (ji,1) = 0._wp 
     184            t_s_1d (ji,1) = rtt 
     185         END IF 
     186      END DO 
     187 
     188      !------------------------------------------------------------! 
     189      !  2) Computing layer thicknesses and enthalpies.            ! 
     190      !------------------------------------------------------------! 
     191      ! 
     192      DO ji = kideb, kiut      
     193         zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 
     194      END DO 
     195      ! 
    181196      DO jk = 1, nlay_s 
    182197         DO ji = kideb, kiut 
    183             zqt_s(ji) =  zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) 
     198            zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji) 
    184199         END DO 
    185200      END DO 
    186201      ! 
    187       zqt_i(:) = 0._wp        ! Total enthalpy of the ice 
    188202      DO jk = 1, nlay_i 
    189203         DO ji = kideb, kiut 
    190             zzc = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 
    191             zqt_i(ji)        =  zqt_i(ji) + zzc 
    192             zqt_i_lay(ji,jk) =              zzc 
     204            zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 
     205            zqh_i(ji)   = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) 
    193206         END DO 
    194207      END DO 
    … …  
    212225      ! Martin Vancoppenolle, December 2006 
    213226 
    214       ! Snow fall 
    215       DO ji = kideb, kiut 
    216          zcoeff = ( 1.0 - ( 1.0 - at_i_b(ji) )**betas ) / at_i_b(ji)  
     227      DO ji = kideb, kiut 
     228         !----------- 
     229         ! Snow fall 
     230         !----------- 
     231         ! thickness change 
     232         zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**betas ) / at_i_1d(ji)  
    217233         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn 
    218       END DO 
    219       zdh_s_mel(:) =  0._wp 
    220  
    221       ! Melt of fallen snow 
    222       DO ji = kideb, kiut 
    223          ! tatm_ice is now in K 
    224          zqprec   (ji)   =  rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus )   
    225          zqfont_su(ji)   =  z_f_surf(ji) * rdt_ice 
    226          zdeltah  (ji,1) =  MIN( 0.e0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) ) 
    227          zqfont_su(ji)   =  MAX( 0.e0 , - zdh_s_pre(ji) - zdeltah(ji,1)              ) * zqprec(ji) 
    228          zdeltah  (ji,1) =  MAX( - zdh_s_pre(ji) , zdeltah(ji,1) ) 
    229          zdh_s_mel(ji)   =  zdh_s_mel(ji) + zdeltah(ji,1) 
    230          ! heat conservation 
    231          qt_s_in(ji,jl)  =  qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji) 
    232          zqt_s  (ji)     =  zqt_s  (ji)    + zqprec(ji) * zdh_s_pre(ji) 
    233          zqt_s  (ji)     =  MAX( zqt_s(ji) - zqfont_su(ji) , 0.e0 )  
    234       END DO 
    235  
    236  
    237       ! Snow melt due to surface heat imbalance 
     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 )    
     236         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 
     237         ! heat flux from snow precip (>0, W.m-2) 
     238         hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 
     239         ! mass flux, <0 
     240         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) ) 
     243 
     244         !--------------------- 
     245         ! Melt of falling snow 
     246         !--------------------- 
     247         ! 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  
     252         ! 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 
     254         ! 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 
    238267      DO jk = 1, nlay_s 
    239268         DO ji = kideb, kiut 
    240             zdeltah  (ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk) 
    241             zqfont_su(ji)    =  MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * q_s_b(ji,jk)  
    242             zdeltah  (ji,jk) =  MAX( zdeltah(ji,jk) , - zh_s(ji) ) 
    243             zdh_s_mel(ji)    =  zdh_s_mel(ji) + zdeltah(ji,jk)        ! resulting melt of snow     
     269            ! thickness change 
     270            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 ) ) )  
     272            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 
     274            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)     
     275            ! heat used to melt snow(W.m-2, >0) 
     276            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice  
     277            ! snow melting only = water into the ocean (then without snow precip) 
     278            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     279 
     280            ! updates available heat + thickness 
     281            zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 
     282            ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 
     283 
    244284         END DO 
    245285      END DO 
    246286 
    247       ! Apply snow melt to snow depth 
    248       DO ji = kideb, kiut 
    249          dh_s_tot(ji)   =  zdh_s_mel(ji) + zdh_s_pre(ji) 
    250          ! Old and new snow depths 
    251          zhsold(ji)     =  ht_s_b(ji) 
    252          zhsnew         =  ht_s_b(ji) + dh_s_tot(ji) 
    253          ! If snow is still present zhn = 1, else zhn = 0 
    254          zhn            =  1.0 - MAX(  zzero , SIGN( zone , - zhsnew )  ) 
    255          ht_s_b(ji)     =  MAX( zzero , zhsnew ) 
    256          ! we recompute dh_s_tot (clem)  
    257          dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji) 
    258          ! Volume and mass variations of snow 
    259          dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) ) 
    260          dvsbq_1d  (ji) =  MIN( zzero, dvsbq_1d(ji) ) 
    261          !clem rdm_snw_1d(ji) =  rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji) 
     287      !---------------------- 
     288      ! 3.2 Snow sublimation  
     289      !---------------------- 
     290      ! 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) 
     292      ! 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 
     313      ! --- Update snow diags --- ! 
     314      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 ) 
    262317      END DO ! ji 
    263318 
     319      !------------------------------------------- 
     320      ! 3.3 Update temperature, energy 
     321      !------------------------------------------- 
     322      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow) 
     323      zq_s(:) = 0._wp  
     324      DO jk = 1, nlay_s 
     325         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 ) ) 
     330            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk) 
     331         END DO 
     332      END DO 
     333 
    264334      !-------------------------- 
    265       ! 3.2 Surface ice ablation  
     335      ! 3.4 Surface ice ablation  
    266336      !-------------------------- 
    267       DO ji = kideb, kiut  
    268          z_f_surf (ji) =  zqfont_su(ji) * r1_rdtice   ! heat conservation test 
    269          zdq_i    (ji) =  0._wp 
    270       END DO ! ji 
    271  
     337      zdeltah(:,:) = 0._wp ! important 
    272338      DO jk = 1, nlay_i 
    273339         DO ji = kideb, kiut  
    274             !                                                    ! melt of layer jk 
    275             zdeltah  (ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk) 
    276             !                                                    ! recompute heat available 
    277             zqfont_su(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk)  
    278             !                                                    ! melt of layer jk cannot be higher than its thickness 
    279             zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) ) 
    280             !                                                    ! update surface melt 
    281             dh_i_surf(ji   ) = dh_i_surf(ji) + zdeltah(ji,jk)  
    282             !                                                    ! for energy conservation 
    283             zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 
    284             ! 
    285             ! clem 
    286             sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    & 
    287                &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 
     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            
     374            ! record which layers have disappeared (for bottom melting)  
     375            !    => icount=0 : no layer has vanished 
     376            !    => 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) ) 
     380 
     381            ! update heat content (J.m-2) and layer thickness 
     382            qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 
     383            h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 
    288384         END DO 
    289385      END DO 
    290  
    291       !                                          !------------------- 
    292       IF( con_i .AND. jiindex_1d > 0 ) THEN      ! Conservation test 
    293          !                                       !------------------- 
    294          numce_dh  = 0 
    295          meance_dh = 0._wp 
    296          DO ji = kideb, kiut 
    297             IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN 
    298                numce_dh  = numce_dh + 1 
    299                meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji) 
    300             ENDIF 
    301             IF( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN! 
    302                WRITE(numout,*) ' ALERTE heat loss for surface melt ' 
    303                WRITE(numout,*) ' ii, ij, jl :', ii, ij, jl 
    304                WRITE(numout,*) ' ht_i_b       : ', ht_i_b(ji) 
    305                WRITE(numout,*) ' z_f_surf     : ', z_f_surf(ji) 
    306                WRITE(numout,*) ' zdq_i        : ', zdq_i(ji) 
    307                WRITE(numout,*) ' ht_i_b       : ', ht_i_b(ji) 
    308                WRITE(numout,*) ' fc_bo_i      : ', fc_bo_i(ji) 
    309                WRITE(numout,*) ' fbif_1d      : ', fbif_1d(ji) 
    310                WRITE(numout,*) ' qlbbq_1d     : ', qlbbq_1d(ji) 
    311                WRITE(numout,*) ' s_i_new      : ', s_i_new(ji) 
    312                WRITE(numout,*) ' sss_m        : ', sss_m(ii,ij) 
    313             ENDIF 
    314          END DO 
    315          ! 
    316          IF( numce_dh > 0 )   meance_dh = meance_dh / numce_dh 
    317          WRITE(numout,*) ' Error report - Category : ', jl 
    318          WRITE(numout,*) ' ~~~~~~~~~~~~ ' 
    319          WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh 
    320          WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh 
    321          ! 
    322       ENDIF 
    323  
    324       !---------------------- 
    325       ! 3.3 Snow sublimation 
    326       !---------------------- 
    327  
    328       DO ji = kideb, kiut 
    329          ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 
    330 #if defined key_coupled 
    331          zdh_s_sub(ji)    =  0._wp      ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) 
    332 #else 
    333          !                              ! forced  mode: snow thickness change due to sublimation 
    334          zdh_s_sub(ji)    =  - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice 
    335 #endif 
    336          dh_s_tot (ji)    =  dh_s_tot(ji) + zdh_s_sub(ji) 
    337          zdhcf            =  ht_s_b(ji) + zdh_s_sub(ji)  
    338          ht_s_b   (ji)    =  MAX( zzero , zdhcf ) 
    339          ! we recompute dh_s_tot  
    340          dh_s_tot (ji)    =  ht_s_b(ji) - zhsold(ji) 
    341          qt_s_in  (ji,jl) =  qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1) 
    342       END DO 
    343  
    344       zqt_dummy(:) = 0.e0 
    345       DO jk = 1, nlay_s 
    346          DO ji = kideb,kiut 
    347             q_s_b    (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 
    348             zqt_dummy(ji)    =  zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s )            ! heat conservation 
    349          END DO 
    350       END DO 
    351  
    352       DO jk = 1, nlay_s 
    353          DO ji = kideb, kiut 
    354             ! In case of disparition of the snow, we have to update the snow temperatures 
    355             zhisn  =  MAX(  zzero , SIGN( zone, - ht_s_b(ji) )  ) 
    356             t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt 
    357             q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk) 
    358          END DO 
     386      ! update ice thickness 
     387      DO ji = kideb, kiut 
     388         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) ) 
    359389      END DO 
    360390 
    … …  
    364394      !------------------------------------------------------------------------------! 
    365395      ! 
    366       ! Ice basal growth / melt is given by the ratio of heat budget over basal 
    367       ! ice heat content.  Basal heat budget is given by the difference between 
    368       ! the inner conductive flux  (fc_bo_i), from the open water heat flux  
    369       ! (qlbbqb) and the turbulent ocean flux (fbif).  
    370       ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice  
    371  
    372       !----------------------------------------------------- 
    373       ! 4.1 Basal growth - (a) salinity not varying in time  
    374       !----------------------------------------------------- 
    375       IF(  num_sal /= 2  ) THEN   ! ice salinity constant in time 
     396      !------------------ 
     397      ! 4.1 Basal growth  
     398      !------------------ 
     399      ! Basal growth is driven by heat imbalance at the ice-ocean interface, 
     400      ! between the inner conductive flux  (fc_bo_i), from the open water heat flux  
     401      ! (fhld) and the turbulent ocean flux (fhtur).  
     402      ! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice  
     403 
     404      ! If salinity varies in time, an iterative procedure is required, because 
     405      ! the involved quantities are inter-dependent. 
     406      ! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi), 
     407      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott) 
     408      ! -> need for an iterative procedure, which converges quickly 
     409 
     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 
     420 
     421      ! Iterative procedure 
     422      DO iter = 1, num_iter_max 
    376423         DO ji = kideb, kiut 
    377             IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0._wp  ) THEN 
    378                s_i_new(ji)         =  sm_i_b(ji) 
    379                ! Melting point in K 
    380                ztmelts             =  - tmut * s_i_new(ji) + rtt  
    381                ! New ice heat content (Bitz and Lipscomb, 1999) 
    382                ztform              =  t_i_b(ji,nlay_i)  ! t_bo_b crashes in the 
    383                ! Baltic 
    384                q_i_b(ji,nlay_i+1)  = rhoic * (  cpic * ( ztmelts - ztform )                                & 
    385                   &                           + lfus * (  1.0 - ( ztmelts - rtt ) / ( ztform - rtt )  )    & 
    386                   &                           - rcp  * ( ztmelts - rtt )                                 ) 
    387                ! Basal growth rate = - F*dt / q 
    388                dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)  
    389                sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 
    390             ENDIF 
    391          END DO 
    392       ENDIF 
    393  
    394       !------------------------------------------------- 
    395       ! 4.1 Basal growth - (b) salinity varying in time  
    396       !------------------------------------------------- 
    397       IF(  num_sal == 2  ) THEN 
    398          ! the growth rate (dh_i_bott) is function of the new ice heat content (q_i_b(nlay_i+1)).  
    399          ! q_i_b depends on the new ice salinity (snewice).  
    400          ! snewice depends on dh_i_bott ; it converges quickly, so, no problem 
    401          ! See Vancoppenolle et al., OM08 for more info on this 
    402  
    403          ! Initial value (tested 1D, can be anything between 1 and 20) 
    404          num_iter_max = 4 
    405          s_i_new(:)   = 4.0 
    406  
    407          ! Iterative procedure 
    408          DO iter = 1, num_iter_max 
    409             DO ji = kideb, kiut 
    410                IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN 
    411                   ii = MOD( npb(ji) - 1, jpi ) + 1 
    412                   ij = ( npb(ji) - 1 ) / jpi + 1 
    413                   ! Melting point in K 
    414                   ztmelts             =   - tmut * s_i_new(ji) + rtt  
    415                   ! New ice heat content (Bitz and Lipscomb, 1999) 
    416                   q_i_b(ji,nlay_i+1)  =  rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             & 
    417                      &                            + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   & 
    418                      &                            - rcp * ( ztmelts-rtt )                                     ) 
    419                   ! Bottom growth rate = - F*dt / q 
    420                   dh_i_bott(ji) =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
    421                   ! New ice salinity ( Cox and Weeks, JGR, 1988 ) 
    422                   ! zswi2  (1) if dh_i_bott/rdt .GT. 3.6e-7 
    423                   ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8 
    424                   ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8 
    425                   zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi13 ) ) 
    426                   zswi2  = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) )  
    427                   zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 
    428                   zswi1  = 1. - zswi2 * zswi12  
    429                   zfracs = zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   & 
    430                      &                   + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  
    431                   zfracs = MIN( 0.5 , zfracs ) 
    432                   s_i_new(ji) = zfracs * sss_m(ii,ij) 
    433                ENDIF ! fc_bo_i 
    434             END DO ! ji 
    435          END DO ! iter 
    436  
    437          ! Final values 
    438          DO ji = kideb, kiut 
    439             IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN 
    440                ! New ice salinity must not exceed 20 psu 
    441                s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) 
    442                ! Metling point in K 
    443                ztmelts     =   - tmut * s_i_new(ji) + rtt  
    444                ! New ice heat content (Bitz and Lipscomb, 1999) 
    445                q_i_b(ji,nlay_i+1)  =  rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                              & 
    446                   &                            + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )    & 
    447                   &                            - rcp * ( ztmelts - rtt )                                    ) 
    448                ! Basal growth rate = - F*dt / q 
    449                dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
    450                ! Salinity update 
    451                ! entrapment during bottom growth 
    452                sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 
    453             ENDIF ! heat budget 
    454          END DO 
    455       ENDIF 
     424            IF(  zf_tt(ji) < 0._wp  ) THEN 
     425 
     426               ! New bottom ice salinity (Cox & Weeks, JGR88 ) 
     427               !--- zswi1  if dh/dt < 2.0e-8 
     428               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7  
     429               !--- zswi2  if dh/dt > 3.6e-7 
     430               zgrr               = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) ) 
     431               zswi2              = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) ) 
     432               zswi12             = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 
     433               zswi1              = 1. - zswi2 * zswi12 
     434               zfracs             = MIN ( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   & 
     435                  &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 ) 
     436 
     437               ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
     438 
     439               s_i_new(ji)        = zswitch_sal * zfracs * sss_m(ii,ij)  &  ! New ice salinity 
     440                                  + ( 1. - zswitch_sal ) * sm_i_1d(ji)  
     441               ! New ice growth 
     442               ztmelts            = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) 
     443 
     444               zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
     445                
     446               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 )           
     449 
     450               zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     451 
     452               zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0) 
     453 
     454               dh_i_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 
     455 
     456               q_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0) 
     457                
     458            ENDIF ! fc_bo_i 
     459         END DO ! ji 
     460      END DO ! iter 
     461 
     462      ! Contribution to Energy and Salt Fluxes 
     463      DO ji = kideb, kiut 
     464         IF(  zf_tt(ji) < 0._wp  ) THEN 
     465            ! New ice growth 
     466                                     
     467            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0) 
     468 
     469            ztmelts        = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) 
     470             
     471            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
     472             
     473            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 )           
     476             
     477            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     478             
     479            zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0) 
     480             
     481            ! Contribution to heat flux to the ocean [W.m-2], >0   
     482            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
     483 
     484            ! Total heat flux used in this process [W.m-2], <0   
     485            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
     486             
     487            ! 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 
     489 
     490            ! Contribution to mass flux, <0 
     491            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice 
     492 
     493            ! update heat content (J.m-2) and layer thickness 
     494            qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_1d(ji,nlay_i+1) 
     495            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji) 
     496         ENDIF 
     497      END DO 
    456498 
    457499      !---------------- 
    458500      ! 4.2 Basal melt 
    459501      !---------------- 
    460       meance_dh = 0._wp 
    461       numce_dh  = 0 
    462       zinnermelt(:) = 0._wp 
    463  
    464       DO ji = kideb, kiut 
    465          ! heat convergence at the surface > 0 
    466          IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp  ) THEN 
    467             s_i_new(ji)   =  s_i_b(ji,nlay_i) 
    468             zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) 
    469             zfbase(ji)    =  zqfont_bo(ji) * r1_rdtice     ! heat conservation test 
    470             zdq_i(ji)     =  0._wp 
    471             dh_i_bott(ji) =  0._wp 
    472          ENDIF 
    473       END DO 
    474  
     502      zdeltah(:,:) = 0._wp ! important 
    475503      DO jk = nlay_i, 1, -1 
    476504         DO ji = kideb, kiut 
    477             IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  >=  0._wp  ) THEN 
    478                ztmelts = - tmut * s_i_b(ji,jk) + rtt  
    479                IF( t_i_b(ji,jk) >= ztmelts ) THEN   !!gm : a comment is needed 
    480                   zdeltah   (ji,jk) = - zh_i(ji) 
    481                   dh_i_bott (ji   ) = dh_i_bott(ji) + zdeltah(ji,jk) 
    482                   zinnermelt(ji   ) = 1._wp 
    483                ELSE                                  ! normal ablation 
    484                   zdeltah  (ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk) 
    485                   zqfont_bo(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 
    486                   zdeltah  (ji,jk) = MAX(zdeltah(ji,jk), - zh_i(ji) ) 
    487                   dh_i_bott(ji   ) = dh_i_bott(ji) + zdeltah(ji,jk) 
    488                   zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 
     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 
     509               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 
     510 
     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 
     515                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
     516                                                                    ! 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 
     521                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk) 
     522 
     523                  zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
     524 
     525                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)  
     526                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 
     527 
     528                  ! 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                                     
     531                  ! Contribution to mass flux 
     532                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     533 
     534                  ! update heat content (J.m-2) and layer thickness 
     535                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 
     536                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 
     537 
     538               ELSE                               !!! Basal melting 
     539 
     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 
     551                   
     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 
     559 
     560                  ! 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 
     562 
     563                  ! 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 
     565                   
     566                  ! 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 
     568                   
     569                  ! Contribution to mass flux 
     570                  wfx_bom_1d(ji) =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     571 
     572                  ! update heat content (J.m-2) and layer thickness 
     573                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 
     574                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 
    489575               ENDIF 
    490                ! clem: contribution to salt flux 
    491                sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    & 
    492                     &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 
     576            
    493577            ENDIF 
    494578         END DO ! ji 
    495579      END DO ! jk 
    496580 
    497       !                                          !------------------- 
    498       IF( con_i .AND. jiindex_1d > 0 ) THEN      ! Conservation test 
    499       !                                          !------------------- 
    500          DO ji = kideb, kiut 
    501             IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0  ) THEN 
    502                IF( ( zfbase(ji) + zdq_i(ji) ) >= 1.e-3 ) THEN 
    503                   numce_dh  = numce_dh + 1 
    504                   meance_dh = meance_dh + zfbase(ji) + zdq_i(ji) 
    505                ENDIF 
    506                IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN 
    507                   WRITE(numout,*) ' ALERTE heat loss for basal melt : ii, ij, jl :', ii, ij, jl 
    508                   WRITE(numout,*) ' ht_i_b    : ', ht_i_b(ji) 
    509                   WRITE(numout,*) ' zfbase    : ', zfbase(ji) 
    510                   WRITE(numout,*) ' zdq_i     : ', zdq_i(ji) 
    511                   WRITE(numout,*) ' ht_i_b    : ', ht_i_b(ji) 
    512                   WRITE(numout,*) ' fc_bo_i   : ', fc_bo_i(ji) 
    513                   WRITE(numout,*) ' fbif_1d   : ', fbif_1d(ji) 
    514                   WRITE(numout,*) ' qlbbq_1d  : ', qlbbq_1d(ji) 
    515                   WRITE(numout,*) ' s_i_new   : ', s_i_new(ji) 
    516                   WRITE(numout,*) ' sss_m     : ', sss_m(ii,ij) 
    517                   WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 
    518                   WRITE(numout,*) ' innermelt : ', INT( zinnermelt(ji) ) 
    519                ENDIF 
    520             ENDIF 
    521          END DO 
    522          IF( numce_dh > 0 )   meance_dh = meance_dh / numce_dh 
    523          WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh 
    524          WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh 
    525          WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d) 
    526          ! 
    527       ENDIF 
    528  
    529       ! 
    530581      !------------------------------------------------------------------------------! 
    531       !  5) Pathological cases                                                       ! 
     582      ! Excessive ablation in a 1-category model 
     583      !     in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15) 
    532584      !------------------------------------------------------------------------------! 
    533       ! 
    534       !---------------------------------------------- 
    535       ! 5.1 Excessive ablation in a 1-category model 
    536       !---------------------------------------------- 
    537  
    538       DO ji = kideb, kiut 
    539          !                     ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15) 
    540          IF( jpl == 1 ) THEN   ;   zdhbf = MAX( hmelt , dh_i_bott(ji) ) 
    541          ELSE                  ;   zdhbf =              dh_i_bott(ji)  
    542          ENDIF 
    543          zdvres        = zdhbf - dh_i_bott(ji) 
    544          dh_i_bott(ji) = zdhbf 
    545          sfx_thd_1d(ji)  = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice 
    546          !                     ! excessive energy is sent to lateral ablation 
    547          zinda = MAX( 0._wp, SIGN( 1._wp , 1.0 - at_i_b(ji) - epsi10 ) ) 
    548          fsup(ji) =  zinda * rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi10 ) * zdvres * r1_rdtice 
    549       END DO 
    550  
    551       !----------------------------------- 
    552       ! 5.2 More than available ice melts 
    553       !----------------------------------- 
    554       ! then heat applied minus heat content at previous time step should equal heat remaining  
    555       ! 
    556       DO ji = kideb, kiut 
    557          ! Adapt the remaining energy if too much ice melts 
    558          !-------------------------------------------------- 
    559          zdvres     = MAX( 0._wp, - ht_i_b(ji) - dh_i_surf(ji) - dh_i_bott(ji) ) 
    560          zdvsur     = MIN( 0._wp, dh_i_surf(ji) + zdvres ) - dh_i_surf(ji) ! fill the surface first 
    561          zdvbot     = MAX( 0._wp, zdvres - zdvsur ) ! then the bottom 
    562          dh_i_surf (ji) = dh_i_surf(ji) + zdvsur ! clem 
    563          dh_i_bott (ji) = dh_i_bott(ji) + zdvbot ! clem 
    564  
    565          ! new ice thickness (clem) 
    566          zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 
    567          zihgnew    = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 
    568          zhgnew(ji) = zihgnew * zhgnew(ji)      ! ice thickness is put to 0 
    569   
    570          !                     !since ice volume is only used for outputs, we keep it global for all categories 
    571          dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 
    572  
    573         ! remaining heat 
    574          zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) )  
    575  
    576          ! If snow remains, energy is used to melt snow 
    577          zhni =  ht_s_b(ji)      ! snow depth at previous time step 
    578          zihg =  MAX(  zzero , SIGN ( zone , - ht_s_b(ji) )  )   ! =0 if snow  
    579  
    580          ! energy of melting of remaining snow 
    581          zinda = MAX( 0._wp, SIGN( 1._wp , zhni - epsi10 ) ) 
    582          zqt_s(ji) =    ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi10 ) * zinda 
    583          zdhnm     =  - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 ) 
    584          zhnfi     =  zhni + zdhnm 
    585          zfdt_final(ji) =  MAX( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 ) 
    586          ht_s_b(ji)     =  MAX( zzero , zhnfi ) 
    587          zqt_s(ji)      =  zqt_s(ji) * ht_s_b(ji) 
    588          ! we recompute dh_s_tot (clem) 
    589          dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji) 
    590  
    591          ! Mass variations of ice and snow 
    592          !--------------------------------- 
    593          !                                              ! mass variation of the jl category 
    594          zzfmass_s = - a_i_b(ji) * ( zhni       - ht_s_b(ji) ) * rhosn   ! snow 
    595          zzfmass_i =   a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic   ! ice   
    596          ! 
    597          zfmass_i(ji) = zzfmass_i                       ! ice variation saved to compute salt flux (see below) 
    598          ! 
    599          !                                              ! mass variation cumulated over category 
    600          !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s                     ! snow  
    601          !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i                     ! ice  
    602  
    603          ! Remaining heat to the ocean  
    604          !--------------------------------- 
    605          focea(ji)  = - zfdt_final(ji) * r1_rdtice         ! focea is in W.m-2 * dt 
    606  
    607          ! residual salt flux (clem) 
    608          !-------------------------- 
    609          ! surface 
    610          sfx_thd_1d(ji)    = sfx_thd_1d(ji) - sm_i_b(ji)  * a_i_b(ji) * zdvsur * rhoic * r1_rdtice 
    611          ! bottom 
    612          IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) >= 0._wp ) THEN ! melting 
    613             sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji)  * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 
    614          ELSE                                                          ! growth 
    615             sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 
    616          ENDIF 
    617          ! 
    618          ! diagnostic  
    619          ii = MOD( npb(ji) - 1, jpi ) + 1 
    620          ij = ( npb(ji) - 1 ) / jpi + 1 
    621          diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 
    622          diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 
    623          diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 
    624       END DO 
    625  
    626       ftotal_fin (:) = zfdt_final(:)  * r1_rdtice 
    627  
    628       !--------------------------- 
    629       ! heat fluxes                     
    630       !--------------------------- 
    631       DO ji = kideb, kiut 
    632          zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice 
    633          ! 
    634          ! Heat flux 
    635          ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1 
    636          ! excessive total  ablation energy (focea) sent to the ocean 
    637          qfvbq_1d(ji)  = qfvbq_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice 
    638  
    639          zihic   = 1.0 - MAX(  zzero , SIGN( zone , -ht_i_b(ji) )  )      ! equals 0 if ht_i = 0, 1 if ht_i gt 0 
    640          fscbq_1d(ji) =  a_i_b(ji) * fstbif_1d(ji) 
    641          qldif_1d(ji)  = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea   (ji) * a_i_b(ji) * rdt_ice   & 
    642             &                                    + ( 1.0 - zihic   ) * fscbq_1d(ji)             * rdt_ice 
    643       END DO  ! ji 
    644  
    645       !------------------------------------------- 
    646       ! Correct temperature, energy and thickness 
    647       !------------------------------------------- 
    648       DO ji = kideb, kiut 
    649          zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )  
    650          t_su_b(ji) =  zihgnew * t_su_b(ji) + ( 1.0 - zihgnew ) * rtt 
    651       END DO  ! ji 
    652  
    653       DO jk = 1, nlay_i 
    654          DO ji = kideb, kiut 
    655             zihgnew      =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )  
    656             t_i_b(ji,jk) =  zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt 
    657             q_i_b(ji,jk) =  zihgnew * q_i_b(ji,jk) 
    658          END DO 
    659       END DO  ! ji 
    660  
    661       DO ji = kideb, kiut 
    662          ht_i_b(ji) = zhgnew(ji) 
    663       END DO  ! ji 
     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 
     605 
     606      !------------------------------------------- 
     607      ! Update temperature, energy 
     608      !------------------------------------------- 
     609      DO ji = kideb, kiut 
     610         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) ) 
     611      END DO   
     612 
     613      !------------------------------------------- 
     614      ! 5. What to do with remaining energy 
     615      !------------------------------------------- 
     616      ! If heat still available for melting and snow remains, then melt more snow 
     617      !------------------------------------------- 
     618      zdeltah(:,:) = 0._wp ! important 
     619      DO ji = kideb, kiut 
     620         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!     
     635         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
     636         ! 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) 
     640      END DO 
     641       
    664642      ! 
    665643      !------------------------------------------------------------------------------| 
    … …  
    670648      DO ji = kideb, kiut 
    671649         ! 
    672          dh_snowice(ji) = MAX(  zzero , ( rhosn * ht_s_b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic )  ) 
    673          zhgnew(ji)     = MAX(  zhgnew(ji) , zhgnew(ji) + dh_snowice(ji)  ) 
    674          zhnnew         = MIN(  ht_s_b(ji) , ht_s_b(ji) - dh_snowice(ji)  ) 
    675  
    676          !  Changes in ice volume and ice mass. 
    677          dvnbq_1d  (ji) =                a_i_b(ji) * ( zhgnew(ji)-ht_i_b(ji) ) 
    678          dmgwi_1d  (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn 
    679  
    680          !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic  
    681          !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn  
    682  
    683          !        Equivalent salt flux (1) Snow-ice formation component 
    684          !        ----------------------------------------------------- 
    685          ii = MOD( npb(ji) - 1, jpi ) + 1 
    686          ij =    ( npb(ji) - 1 ) / jpi + 1 
    687  
    688          IF( num_sal == 2 ) THEN   ;   zsm_snowice = sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic 
    689          ELSE                      ;   zsm_snowice = sm_i_b(ji)    
    690          ENDIF 
     650         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic )  ) 
     651 
     652         ht_i_1d(ji)     = ht_i_1d(ji) + dh_snowice(ji) 
     653         ht_s_1d(ji)     = ht_s_1d(ji) - dh_snowice(ji) 
     654 
     655         ! Salinity of snow ice 
     656         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) 
     658 
    691659         ! entrapment during snow ice formation 
    692          ! clem: new salinity difference stored (to be used in limthd_ent.F90) 
     660         ! new salinity difference stored (to be used in limthd_ent.F90) 
    693661         IF (  num_sal == 2  ) THEN 
    694             i_ice_switch = MAX( 0._wp , SIGN( 1._wp , zhgnew(ji) - epsi10 ) ) 
     662            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) ) 
    695663            ! salinity dif due to snow-ice formation 
    696             dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch      
     664            dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi10 ) * rswitch      
    697665            ! salinity dif due to bottom growth  
    698             IF (  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  < 0._wp ) THEN 
    699                dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch 
     666            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 
    700668            ENDIF 
    701669         ENDIF 
    702670 
    703          !  Actualize new snow and ice thickness. 
    704          ht_s_b(ji)  = zhnnew 
    705          ht_i_b(ji)  = zhgnew(ji) 
    706  
    707          ! Total ablation ! new lines added to debug 
    708          IF( ht_i_b(ji) <= 0._wp )   a_i_b(ji) = 0._wp 
    709  
    710          ! diagnostic ( snow ice growth ) 
    711          ii = MOD( npb(ji) - 1, jpi ) + 1 
    712          ij =    ( npb(ji) - 1 ) / jpi + 1 
    713          diag_sni_gr(ii,ij)  = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice 
    714          ! 
    715          ! salt flux 
    716          sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
    717          !-------------------------------- 
    718          ! Update mass fluxes (clem) 
    719          !-------------------------------- 
    720          rdm_ice_1d(ji) = rdm_ice_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic  
    721          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn  
     671         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0) 
     672         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
     673         zfmdt          = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0._wp )    ! <0 
     674         zsstK          = sst_m(ii,ij) + rt0                                 
     675         zEw            = rcp * ( zsstK - rt0 ) 
     676         zQm            = zfmdt * zEw  
     677          
     678         ! Contribution to heat flux 
     679         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice  
     680 
     681         ! Contribution to salt flux 
     682         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice  
     683           
     684         ! Contribution to mass flux 
     685         ! All snow is thrown in the ocean, and seawater is taken to replace the volume 
     686         wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
     687         wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice 
     688 
     689         ! update heat content (J.m-2) and layer thickness 
     690         qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw 
     691         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 
     692          
     693         ! Total ablation (to debug) 
     694         IF( ht_i_1d(ji) <= 0._wp )   a_i_1d(ji) = 0._wp 
    722695 
    723696      END DO !ji 
    724       ! 
    725       CALL wrk_dealloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 
    726       CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 
    727       CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i ) 
    728       CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay ) 
    729       ! 
    730       CALL wrk_dealloc( jpij, zviold, zvsold ) ! clem 
     697 
     698      ! 
     699      !------------------------------------------- 
     700      ! Update temperature, energy 
     701      !------------------------------------------- 
     702      !clem bug: we should take snow into account here 
     703      DO ji = kideb, kiut 
     704         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 
     707 
     708      DO jk = 1, nlay_s 
     709         DO ji = kideb,kiut 
     710            ! 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) 
     713            ! 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 ) 
     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 ) 
     722      ! 
    731723      ! 
    732724   END SUBROUTINE lim_thd_dh 
Note: See TracChangeset for help on using the changeset viewer.