MODULE limthd_dh !!====================================================================== !! *** MODULE limthd_dh *** !! LIM-3 : thermodynamic growth and decay of the ice !!====================================================================== !! History : LIM ! 2003-05 (M. Vancoppenolle) Original code in 1D !! ! 2005-06 (M. Vancoppenolle) 3D version !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw & wfx_ice !! 3.4 ! 2011-02 (G. Madec) dynamical allocation !! 3.5 ! 2012-10 (G. Madec & co) salt flux + bug fixes !!---------------------------------------------------------------------- #if defined key_lim3 !!---------------------------------------------------------------------- !! 'key_lim3' LIM3 sea-ice model !!---------------------------------------------------------------------- !! lim_thd_dh : vertical accr./abl. and lateral ablation of sea ice !!---------------------------------------------------------------------- USE par_oce ! ocean parameters USE phycst ! physical constants (OCE directory) USE sbc_oce ! Surface boundary condition: ocean fields USE ice ! LIM variables USE par_ice ! LIM parameters USE thd_ice ! LIM thermodynamics USE in_out_manager ! I/O manager USE lib_mpp ! MPP library USE wrk_nemo ! work arrays USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) IMPLICIT NONE PRIVATE PUBLIC lim_thd_dh ! called by lim_thd !!---------------------------------------------------------------------- !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_thd_dh( kideb, kiut ) !!------------------------------------------------------------------ !! *** ROUTINE lim_thd_dh *** !! !! ** Purpose : determines variations of ice and snow thicknesses. !! !! ** Method : Ice/Snow surface melting arises from imbalance in surface fluxes !! Bottom accretion/ablation arises from flux budget !! Snow thickness can increase by precipitation and decrease by sublimation !! If snow load excesses Archmiede limit, snow-ice is formed by !! the flooding of sea-water in the snow !! !! 1) Compute available flux of heat for surface ablation !! 2) Compute snow and sea ice enthalpies !! 3) Surface ablation and sublimation !! 4) Bottom accretion/ablation !! 5) Case of Total ablation !! 6) Snow ice formation !! !! References : Bitz and Lipscomb, 1999, J. Geophys. Res. !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 !! Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let. !! Vancoppenolle et al.,2009, Ocean Modelling !!------------------------------------------------------------------ INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied !! INTEGER :: ji , jk ! dummy loop indices INTEGER :: ii, ij ! 2D corresponding indices to ji INTEGER :: iter REAL(wp) :: ztmelts ! local scalar REAL(wp) :: zdh, zfdum ! REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads REAL(wp) :: zs_snic ! snow-ice salinity REAL(wp) :: zswi1 ! switch for computation of bottom salinity REAL(wp) :: zswi12 ! switch for computation of bottom salinity REAL(wp) :: zswi2 ! switch for computation of bottom salinity REAL(wp) :: zgrr ! bottom growth rate REAL(wp) :: zt_i_new ! bottom formation temperature REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean REAL(wp) :: zEi ! specific enthalpy of sea ice (J/kg) REAL(wp) :: zEw ! specific enthalpy of exchanged water (J/kg) REAL(wp) :: zdE ! specific enthalpy difference (J/kg) REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean REAL(wp) :: zsstK ! SST in Kelvin REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness REAL(wp), POINTER, DIMENSION(:) :: zqprec ! energy of fallen snow (J.m-3) REAL(wp), POINTER, DIMENSION(:) :: zq_su ! heat for surface ablation (J.m-2) REAL(wp), POINTER, DIMENSION(:) :: zq_bo ! heat for bottom ablation (J.m-2) REAL(wp), POINTER, DIMENSION(:) :: zq_1cat ! corrected heat in case 1-cat and hmelt>15cm (J.m-2) REAL(wp), POINTER, DIMENSION(:) :: zq_rema ! remaining heat at the end of the routine (J.m-2) REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2) INTEGER , POINTER, DIMENSION(:) :: icount ! number of layers vanished by melting REAL(wp), POINTER, DIMENSION(:) :: zdh_s_mel ! snow melt REAL(wp), POINTER, DIMENSION(:) :: zdh_s_pre ! snow precipitation REAL(wp), POINTER, DIMENSION(:) :: zdh_s_sub ! snow sublimation REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah REAL(wp), POINTER, DIMENSION(:,:) :: zh_i ! ice layer thickness REAL(wp), POINTER, DIMENSION(:) :: zqh_i ! total ice heat content (J.m-2) REAL(wp), POINTER, DIMENSION(:) :: zqh_s ! total snow heat content (J.m-2) REAL(wp), POINTER, DIMENSION(:) :: zq_s ! total snow enthalpy (J.m-3) ! mass and salt flux (clem) REAL(wp) :: zdvres, zswitch_sal ! Heat conservation INTEGER :: num_iter_max !!------------------------------------------------------------------ ! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values) SELECT CASE( num_sal ) ! varying salinity or not CASE( 1, 3, 4 ) ; zswitch_sal = 0 ! prescribed salinity profile CASE( 2 ) ; zswitch_sal = 1 ! varying salinity profile END SELECT CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) CALL wrk_alloc( jpij, nlay_i+1, zdeltah, zh_i ) CALL wrk_alloc( jpij, icount ) dh_i_surf (:) = 0._wp ; dh_i_bott (:) = 0._wp ; dh_snowice(:) = 0._wp dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp zqprec (:) = 0._wp ; zq_su (:) = 0._wp ; zq_bo (:) = 0._wp ; zf_tt (:) = 0._wp zq_1cat(:) = 0._wp ; zq_rema(:) = 0._wp zh_s (:) = 0._wp zdh_s_pre(:) = 0._wp zdh_s_mel(:) = 0._wp zdh_s_sub(:) = 0._wp zqh_s (:) = 0._wp zqh_i (:) = 0._wp zh_i (:,:) = 0._wp zdeltah (:,:) = 0._wp icount (:) = 0 ! initialize layer thicknesses and enthalpies h_i_old (:,0:nlay_i+1) = 0._wp qh_i_old(:,0:nlay_i+1) = 0._wp DO jk = 1, nlay_i DO ji = kideb, kiut h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) ENDDO ENDDO ! !------------------------------------------------------------------------------! ! 1) Calculate available heat for surface and bottom ablation ! !------------------------------------------------------------------------------! ! DO ji = kideb, kiut rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) ztmelts = rswitch * rtt + ( 1._wp - rswitch ) * rtt zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts ) ) zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) END DO ! !------------------------------------------------------------------------------! ! If snow temperature is above freezing point, then snow melts ! (should not happen but sometimes it does) !------------------------------------------------------------------------------! DO ji = kideb, kiut IF( t_s_1d(ji,1) > rtt ) THEN !!! Internal melting ! Contribution to heat flux to the ocean [W.m-2], < 0 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice ! Contribution to mass flux wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice ! updates ht_s_1d(ji) = 0._wp q_s_1d (ji,1) = 0._wp t_s_1d (ji,1) = rtt END IF END DO !------------------------------------------------------------! ! 2) Computing layer thicknesses and enthalpies. ! !------------------------------------------------------------! ! DO ji = kideb, kiut zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) END DO ! DO jk = 1, nlay_s DO ji = kideb, kiut zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji) END DO END DO ! DO jk = 1, nlay_i DO ji = kideb, kiut zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) zqh_i(ji) = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) END DO END DO ! !------------------------------------------------------------------------------| ! 3) Surface ablation and sublimation | !------------------------------------------------------------------------------| ! !------------------------- ! 3.1 Snow precips / melt !------------------------- ! Snow accumulation in one thermodynamic time step ! snowfall is partitionned between leads and ice ! if snow fall was uniform, a fraction (1-at_i) would fall into leads ! but because of the winds, more snow falls on leads than on sea ice ! and a greater fraction (1-at_i)^beta of the total mass of snow ! (beta < 1) falls in leads. ! In reality, beta depends on wind speed, ! and should decrease with increasing wind speed but here, it is ! considered as a constant. an average value is 0.66 ! Martin Vancoppenolle, December 2006 DO ji = kideb, kiut !----------- ! Snow fall !----------- ! thickness change zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**betas ) / at_i_1d(ji) zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) zqprec (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus ) IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp ! heat flux from snow precip (>0, W.m-2) hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice ! mass flux, <0 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice ! update thickness ht_s_1d (ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) !--------------------- ! Melt of falling snow !--------------------- ! thickness change IF( zdh_s_pre(ji) > 0._wp ) THEN rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) ) zdh_s_mel (ji) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting ! heat used to melt snow (W.m-2, >0) hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice ! snow melting only = water into the ocean (then without snow precip), >0 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_mel(ji) * r1_rdtice ! updates available heat + thickness zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) ) ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) ) zh_s (ji) = ht_s_1d(ji) / REAL( nlay_s ) ENDIF END DO ! If heat still available, then melt more snow zdeltah(:,:) = 0._wp ! important DO jk = 1, nlay_s DO ji = kideb, kiut ! thickness change rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) rswitch = rswitch * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) + epsi20 ) ) ) zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) ! heat used to melt snow(W.m-2, >0) hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice ! snow melting only = water into the ocean (then without snow precip) wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! updates available heat + thickness zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) END DO END DO !---------------------- ! 3.2 Snow sublimation !---------------------- ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates ! clem comment: not counted in mass exchange in limsbc since this is an exchange with atm. (not ocean) ! clem comment: ice should also sublimate IF( lk_cpl ) THEN ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) zdh_s_sub(:) = 0._wp ELSE ! forced mode: snow thickness change due to sublimation DO ji = kideb, kiut zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) ! Heat flux by sublimation [W.m-2], < 0 ! sublimate first snow that had fallen, then pre-existing snow zcoeff = ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) * zqprec(ji) + & & ( 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) ) & & * a_i_1d(ji) * r1_rdtice hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff ! Mass flux by sublimation wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice ! new snow thickness ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) END DO ENDIF ! --- Update snow diags --- ! DO ji = kideb, kiut dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) END DO ! ji !------------------------------------------- ! 3.3 Update temperature, energy !------------------------------------------- ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow) zq_s(:) = 0._wp DO jk = 1, nlay_s DO ji = kideb,kiut rswitch = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) + epsi20 ) ) q_s_1d(ji,jk) = ( 1._wp - rswitch ) / MAX( ht_s_1d(ji), epsi20 ) * & & ( ( MAX( 0._wp, dh_s_tot(ji) ) ) * zqprec(ji) + & & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) ) zq_s(ji) = zq_s(ji) + q_s_1d(ji,jk) END DO END DO !-------------------------- ! 3.4 Surface ice ablation !-------------------------- zdeltah(:,:) = 0._wp ! important DO jk = 1, nlay_i DO ji = kideb, kiut zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0] ztmelts = - tmut * s_i_1d(ji,jk) + rtt ! Melting point of layer k [K] zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] zdE = zEi - zEw ! Specific enthalpy difference < 0 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] zdeltah(ji,jk) = - zfmdt / rhoic ! Melt of layer jk [m, <0] 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] zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice ! Contribution to heat flux [W.m-2], < 0 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Total heat flux used in this process [W.m-2], > 0 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Contribution to mass flux wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! record which layers have disappeared (for bottom melting) ! => icount=0 : no layer has vanished ! => icount=5 : 5 layers have vanished rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) icount(ji) = icount(ji) + NINT( rswitch ) zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) ! update heat content (J.m-2) and layer thickness qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) END DO END DO ! update ice thickness DO ji = kideb, kiut ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) ) END DO ! !------------------------------------------------------------------------------! ! 4) Basal growth / melt ! !------------------------------------------------------------------------------! ! !------------------ ! 4.1 Basal growth !------------------ ! Basal growth is driven by heat imbalance at the ice-ocean interface, ! between the inner conductive flux (fc_bo_i), from the open water heat flux ! (fhld) and the turbulent ocean flux (fhtur). ! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice ! If salinity varies in time, an iterative procedure is required, because ! the involved quantities are inter-dependent. ! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi), ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott) ! -> need for an iterative procedure, which converges quickly IF ( num_sal == 2 ) THEN num_iter_max = 5 ELSE num_iter_max = 1 ENDIF !clem debug. Just to be sure that enthalpy at nlay_i+1 is null DO ji = kideb, kiut q_i_1d(ji,nlay_i+1) = 0._wp END DO ! Iterative procedure DO iter = 1, num_iter_max DO ji = kideb, kiut IF( zf_tt(ji) < 0._wp ) THEN ! New bottom ice salinity (Cox & Weeks, JGR88 ) !--- zswi1 if dh/dt < 2.0e-8 !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7 !--- zswi2 if dh/dt > 3.6e-7 zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) ) zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) ) zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) zswi1 = 1. - zswi2 * zswi12 zfracs = MIN ( zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) & & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 ) ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 s_i_new(ji) = zswitch_sal * zfracs * sss_m(ii,ij) & ! New ice salinity + ( 1. - zswitch_sal ) * sm_i_1d(ji) ! New ice growth ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K) zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) & - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) ) & & + rcp * ( ztmelts-rtt ) zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) q_i_1d(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0) ENDIF ! fc_bo_i END DO ! ji END DO ! iter ! Contribution to Energy and Salt Fluxes DO ji = kideb, kiut IF( zf_tt(ji) < 0._wp ) THEN ! New ice growth zfmdt = - rhoic * dh_i_bott(ji) ! Mass flux x time step (kg/m2, < 0) ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K) zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) & - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) ) & & + rcp * ( ztmelts-rtt ) zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) ! Contribution to heat flux to the ocean [W.m-2], >0 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Total heat flux used in this process [W.m-2], <0 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Contribution to salt flux, <0 sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_1d(ji) * zfmdt * r1_rdtice ! Contribution to mass flux, <0 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice ! update heat content (J.m-2) and layer thickness 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) h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji) ENDIF END DO !---------------- ! 4.2 Basal melt !---------------- zdeltah(:,:) = 0._wp ! important DO jk = nlay_i, 1, -1 DO ji = kideb, kiut IF( zf_tt(ji) >= 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared from surface melting ztmelts = - tmut * s_i_1d(ji,jk) + rtt ! Melting point of layer jk (K) IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) !!zEw = rcp * ( t_i_1d(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) ! set up at 0 since no energy is needed to melt water...(it is already melted) zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing ! this should normally not happen, but sometimes, heat diffusion leads to this dh_i_bott (ji) = dh_i_bott(ji) + zdeltah(ji,jk) zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice ! Contribution to mass flux wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! update heat content (J.m-2) and layer thickness qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) ELSE !!! Basal melting zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) zEw = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0) zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0) zdeltah(ji,jk) = - zfmdt / rhoic ! Gross thickness change zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change zq_bo(ji) = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) ! Update basal melt zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 zQm = zfmdt * zEw ! Heat exchanged with ocean ! Contribution to heat flux to the ocean [W.m-2], <0 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice ! Total heat flux used in this process [W.m-2], >0 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Contribution to mass flux wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! update heat content (J.m-2) and layer thickness qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) ENDIF ENDIF END DO ! ji END DO ! jk !------------------------------------------------------------------------------! ! Excessive ablation in a 1-category model ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15) !------------------------------------------------------------------------------! ! ??? keep ??? ! clem bug: I think this should be included above, so we would not have to ! track heat/salt/mass fluxes backwards ! IF( jpl == 1 ) THEN ! DO ji = kideb, kiut ! IF( zf_tt(ji) >= 0._wp ) THEN ! zdh = MAX( hmelt , dh_i_bott(ji) ) ! zdvres = zdh - dh_i_bott(ji) ! >=0 ! dh_i_bott(ji) = zdh ! ! ! excessive energy is sent to lateral ablation ! rswitch = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_1d(ji) - epsi20 ) ) ! zq_1cat(ji) = rswitch * rhoic * lfus * at_i_1d(ji) / MAX( 1._wp - at_i_1d(ji) , epsi20 ) * zdvres ! J.m-2 >=0 ! ! ! correct salt and mass fluxes ! 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 ! wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdvres * r1_rdtice ! ENDIF ! END DO ! ENDIF !------------------------------------------- ! Update temperature, energy !------------------------------------------- DO ji = kideb, kiut ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) ) END DO !------------------------------------------- ! 5. What to do with remaining energy !------------------------------------------- ! If heat still available for melting and snow remains, then melt more snow !------------------------------------------- zdeltah(:,:) = 0._wp ! important DO ji = kideb, kiut zq_rema(ji) = zq_su(ji) + zq_bo(ji) ! zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) ! =1 if snow ! zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zq_s(ji) + epsi20 ) ) ! zdeltah (ji,1) = - zindh * zindq * zq_rema(ji) / MAX( zq_s(ji), epsi20 ) ! zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting ! zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) ! dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) ! ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1) ! ! zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji) ! update available heat (J.m-2) ! ! heat used to melt snow ! hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0) ! ! Contribution to mass flux ! wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice ! ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice IF( ln_nicep .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) END DO ! !------------------------------------------------------------------------------| ! 6) Snow-Ice formation | !------------------------------------------------------------------------------| ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level, ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice DO ji = kideb, kiut ! dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic ) ) ht_i_1d(ji) = ht_i_1d(ji) + dh_snowice(ji) ht_s_1d(ji) = ht_s_1d(ji) - dh_snowice(ji) ! Salinity of snow ice ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) ! entrapment during snow ice formation ! new salinity difference stored (to be used in limthd_ent.F90) IF ( num_sal == 2 ) THEN rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) ) ! salinity dif due to snow-ice formation dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi10 ) * rswitch ! salinity dif due to bottom growth IF ( zf_tt(ji) < 0._wp ) THEN dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi10 ) * rswitch ENDIF ENDIF ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0) ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 zfmdt = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0._wp ) ! <0 zsstK = sst_m(ii,ij) + rt0 zEw = rcp * ( zsstK - rt0 ) zQm = zfmdt * zEw ! Contribution to heat flux hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Contribution to salt flux sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice ! Contribution to mass flux ! All snow is thrown in the ocean, and seawater is taken to replace the volume wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice ! update heat content (J.m-2) and layer thickness qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) ! Total ablation (to debug) IF( ht_i_1d(ji) <= 0._wp ) a_i_1d(ji) = 0._wp END DO !ji ! !------------------------------------------- ! Update temperature, energy !------------------------------------------- !clem bug: we should take snow into account here DO ji = kideb, kiut rswitch = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rtt END DO ! ji DO jk = 1, nlay_s DO ji = kideb,kiut ! mask enthalpy rswitch = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) ) q_s_1d(ji,jk) = ( 1.0 - rswitch ) * q_s_1d(ji,jk) ! recalculate t_s_1d from q_s_1d t_s_1d(ji,jk) = rtt + ( 1._wp - rswitch ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) END DO END DO CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i ) CALL wrk_dealloc( jpij, icount ) ! ! END SUBROUTINE lim_thd_dh #else !!---------------------------------------------------------------------- !! Default option NO LIM3 sea-ice model !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_thd_dh ! Empty routine END SUBROUTINE lim_thd_dh #endif !!====================================================================== END MODULE limthd_dh