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 rdm_snw & rdm_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 REAL(wp) :: epsi20 = 1.e-20 ! constant values REAL(wp) :: epsi10 = 1.e-10 ! REAL(wp) :: epsi13 = 1.e-13 ! REAL(wp) :: zzero = 0._wp ! REAL(wp) :: zone = 1._wp ! !!---------------------------------------------------------------------- !! 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, jl ) !!------------------------------------------------------------------ !! *** 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 , INTENT(in) :: jl ! Thickness cateogry number !! INTEGER :: ji , jk ! dummy loop indices INTEGER :: ii, ij ! 2D corresponding indices to ji INTEGER :: i_use ! debugging flag INTEGER :: isnow ! switch for presence (1) or absence (0) of snow INTEGER :: isnowic ! snow ice formation not INTEGER :: i_ice_switch ! ice thickness above a certain treshold or not INTEGER :: iter REAL(wp) :: zzfmass_i, zihgnew ! local scalar REAL(wp) :: zzfmass_s, zhsnew, ztmelts ! local scalar REAL(wp) :: zhn, zdhcf, zdhbf, zhni, zhnfi, zihg ! REAL(wp) :: zdhnm, zhnnew, zhisn, zihic, zzc ! 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_i ! ice layer thickness REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness REAL(wp), POINTER, DIMENSION(:) :: ztfs ! melting point REAL(wp), POINTER, DIMENSION(:) :: zhsold ! old snow thickness REAL(wp), POINTER, DIMENSION(:) :: zqprec ! energy of fallen snow REAL(wp), POINTER, DIMENSION(:) :: zqfont_su ! incoming, remaining surface energy REAL(wp), POINTER, DIMENSION(:) :: zqfont_bo ! incoming, bottom energy REAL(wp), POINTER, DIMENSION(:) :: z_f_surf ! surface heat for ablation REAL(wp), POINTER, DIMENSION(:) :: zhgnew ! new ice thickness REAL(wp), POINTER, DIMENSION(:) :: zfmass_i ! 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 ! Pathological cases REAL(wp), POINTER, DIMENSION(:) :: zfdt_init ! total incoming heat for ice melt REAL(wp), POINTER, DIMENSION(:) :: zfdt_final ! total remaing heat for ice melt REAL(wp), POINTER, DIMENSION(:) :: zqt_i ! total ice heat content REAL(wp), POINTER, DIMENSION(:) :: zqt_s ! total snow heat content REAL(wp), POINTER, DIMENSION(:) :: zqt_dummy ! dummy heat content REAL(wp), POINTER, DIMENSION(:,:) :: zqt_i_lay ! total ice heat content ! mass and salt flux (clem) REAL(wp) :: zdvres, zdvsur, zdvbot, zswitch_sal REAL(wp), POINTER, DIMENSION(:) :: zviold, zvsold ! old ice volume... ! Heat conservation INTEGER :: num_iter_max, numce_dh REAL(wp) :: meance_dh REAL(wp) :: zinda REAL(wp), POINTER, DIMENSION(:) :: zinnermelt REAL(wp), POINTER, DIMENSION(:) :: zfbase, zdq_i !!------------------------------------------------------------------ ! 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_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) CALL wrk_alloc( jpij, zinnermelt, zfbase, zdq_i ) CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay ) CALL wrk_alloc( jpij, zviold, zvsold ) ! clem ftotal_fin(:) = 0._wp zfdt_init (:) = 0._wp zfdt_final(:) = 0._wp dh_i_surf (:) = 0._wp dh_i_bott (:) = 0._wp dh_snowice(:) = 0._wp DO ji = kideb, kiut old_ht_i_b(ji) = ht_i_b(ji) old_ht_s_b(ji) = ht_s_b(ji) zviold(ji) = a_i_b(ji) * ht_i_b(ji) ! clem zvsold(ji) = a_i_b(ji) * ht_s_b(ji) ! clem END DO ! !------------------------------------------------------------------------------! ! 1) Calculate available heat for surface ablation ! !------------------------------------------------------------------------------! ! DO ji = kideb, kiut isnow = INT( 1.0 - MAX( 0.0 , SIGN( 1.0 , - ht_s_b(ji) ) ) ) ztfs (ji) = isnow * rtt + ( 1.0 - isnow ) * rtt z_f_surf (ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) z_f_surf (ji) = MAX( zzero , z_f_surf(ji) ) * MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) ) zfdt_init(ji) = ( z_f_surf(ji) + MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) * rdt_ice END DO ! ji zqfont_su (:) = 0._wp zqfont_bo (:) = 0._wp dsm_i_se_1d(:) = 0._wp dsm_i_si_1d(:) = 0._wp ! !------------------------------------------------------------------------------! ! 2) Computing layer thicknesses and snow and sea-ice enthalpies. ! !------------------------------------------------------------------------------! ! DO ji = kideb, kiut ! Layer thickness zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) END DO ! zqt_s(:) = 0._wp ! Total enthalpy of the snow DO jk = 1, nlay_s DO ji = kideb, kiut zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) END DO END DO ! zqt_i(:) = 0._wp ! Total enthalpy of the ice DO jk = 1, nlay_i DO ji = kideb, kiut zzc = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) zqt_i(ji) = zqt_i(ji) + zzc zqt_i_lay(ji,jk) = zzc 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 ! Snow fall DO ji = kideb, kiut zcoeff = ( 1.0 - ( 1.0 - at_i_b(ji) )**betas ) / at_i_b(ji) zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn END DO zdh_s_mel(:) = 0._wp ! Melt of fallen snow DO ji = kideb, kiut ! tatm_ice is now in K zqprec (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus ) zqfont_su(ji) = z_f_surf(ji) * rdt_ice zdeltah (ji,1) = MIN( 0.e0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) ) zqfont_su(ji) = MAX( 0.e0 , - zdh_s_pre(ji) - zdeltah(ji,1) ) * zqprec(ji) zdeltah (ji,1) = MAX( - zdh_s_pre(ji) , zdeltah(ji,1) ) zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) ! heat conservation qt_s_in(ji,jl) = qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji) zqt_s (ji) = zqt_s (ji) + zqprec(ji) * zdh_s_pre(ji) zqt_s (ji) = MAX( zqt_s(ji) - zqfont_su(ji) , 0.e0 ) END DO ! Snow melt due to surface heat imbalance DO jk = 1, nlay_s DO ji = kideb, kiut zdeltah (ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk) zqfont_su(ji) = MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * q_s_b(ji,jk) zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) ! resulting melt of snow END DO END DO ! Apply snow melt to snow depth DO ji = kideb, kiut dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) ! Old and new snow depths zhsold(ji) = ht_s_b(ji) zhsnew = ht_s_b(ji) + dh_s_tot(ji) ! If snow is still present zhn = 1, else zhn = 0 zhn = 1.0 - MAX( zzero , SIGN( zone , - zhsnew ) ) ht_s_b(ji) = MAX( zzero , zhsnew ) ! we recompute dh_s_tot (clem) dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) ! Volume and mass variations of snow dvsbq_1d (ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) ) dvsbq_1d (ji) = MIN( zzero, dvsbq_1d(ji) ) !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji) END DO ! ji !-------------------------- ! 3.2 Surface ice ablation !-------------------------- DO ji = kideb, kiut z_f_surf (ji) = zqfont_su(ji) * r1_rdtice ! heat conservation test zdq_i (ji) = 0._wp END DO ! ji DO jk = 1, nlay_i DO ji = kideb, kiut zEi = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0] ztmelts = - tmut * s_i_b(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 = - zqfont_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] zdeltah(ji,jk) = - zfmdt / rhoic ! Melt of layer jk [m, <0] zqfont_su(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * rhoic * ( - zdE ) ! Energy remaining in case of melting of the full layer [J/m2, >0] zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0] 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 = MAX ( zfmdt, 0._wp ) * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] ! Contribution to salt flux sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) & & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice ! Contribution to heat flux rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji) ! Conservation test zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * rhoic * zdE * r1_rdtice ! ? still don't know END DO END DO ! !------------------- IF( con_i .AND. jiindex_1d > 0 ) THEN ! Conservation test ! !------------------- numce_dh = 0 meance_dh = 0._wp DO ji = kideb, kiut IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN numce_dh = numce_dh + 1 meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji) ENDIF IF( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN! WRITE(numout,*) ' ALERTE heat loss for surface melt ' WRITE(numout,*) ' ii, ij, jl :', ii, ij, jl WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) WRITE(numout,*) ' z_f_surf : ', z_f_surf(ji) WRITE(numout,*) ' zdq_i : ', zdq_i(ji) WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji) WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji) WRITE(numout,*) ' qlbbq_1d : ', qlbbq_1d(ji) WRITE(numout,*) ' s_i_new : ', s_i_new(ji) WRITE(numout,*) ' sss_m : ', sss_m(ii,ij) ENDIF END DO ! IF( numce_dh > 0 ) meance_dh = meance_dh / numce_dh WRITE(numout,*) ' Error report - Category : ', jl WRITE(numout,*) ' ~~~~~~~~~~~~ ' WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh ! ENDIF !---------------------- ! 3.3 Snow sublimation !---------------------- DO ji = kideb, kiut ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates #if defined key_coupled zdh_s_sub(ji) = 0._wp ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) #else ! ! forced mode: snow thickness change due to sublimation zdh_s_sub(ji) = - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice #endif dh_s_tot (ji) = dh_s_tot(ji) + zdh_s_sub(ji) zdhcf = ht_s_b(ji) + zdh_s_sub(ji) ht_s_b (ji) = MAX( zzero , zdhcf ) ! we recompute dh_s_tot dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) qt_s_in (ji,jl) = qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1) END DO zqt_dummy(:) = 0.e0 DO jk = 1, nlay_s DO ji = kideb,kiut q_s_b (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) ! heat conservation END DO END DO DO jk = 1, nlay_s DO ji = kideb, kiut ! In case of disparition of the snow, we have to update the snow temperatures zhisn = MAX( zzero , SIGN( zone, - ht_s_b(ji) ) ) t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk) END DO 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 ! (qlbbqb) and the turbulent ocean flux (fbif). ! fc_bo_i is positive downwards. fbif and qlbbq 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 ! Initialize dh_i_bott dh_i_bott(:) = 0.0e0 ! Iterative procedure DO iter = 1, num_iter_max DO ji = kideb, kiut IF( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0 ) 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 , epsi13 ) ) zswi2 = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) zswi12 = MAX( zzero , SIGN( zone , 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_b(ji) ! New ice growth ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K) zt_i_new = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(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_b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) dh_i_bott(ji) = rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / ( zdE * rhoic ) !!! not sure we still need the next line... useful to keep this in memory ? check limthd_ent... q_i_b(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 zEw = rcp * ( t_bo_b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) zfmdt = - rhoic * dh_i_bott(ji) ! Mass flux x time step (kg/m2, < 0) ! Contribution to energy flux to the ocean (J/m2) rdq_ice_1d(ji) = rdq_ice_1d(ji) + zEw * a_i_b(ji) * zfmdt ! Contribution to salt flux () sfx_thd_1d(ji) = sfx_thd_1d(ji) + s_i_new(ji) * a_i_b(ji) * zfmdt * r1_rdtice END DO !---------------- ! 4.2 Basal melt !---------------- meance_dh = 0._wp numce_dh = 0 zinnermelt(:) = 0._wp DO ji = kideb, kiut ! heat convergence at the surface > 0 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp ) THEN s_i_new(ji) = s_i_b(ji,nlay_i) ! is this line still useful ? zqfont_bo(ji) = rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) zfbase(ji) = zqfont_bo(ji) * r1_rdtice ! heat conservation test zdq_i(ji) = 0._wp dh_i_bott(ji) = 0._wp ENDIF END DO DO jk = nlay_i, 1, -1 DO ji = kideb, kiut IF( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) >= 0._wp ) THEN ztmelts = - tmut * s_i_b(ji,jk) + rtt ! Melting point of layer jk (K) IF( t_i_b(ji,jk) >= ztmelts ) THEN !!! Internal melting zdeltah (ji,jk) = - zh_i(ji) ! 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) zinnermelt(ji) = 1._wp zQm = 0. ! Not sure which specific enthalpy we should use here (MV HC 2014) ! If that happens, heat is probably not well counted ! Put zero by default, but more bug analysis should be done to investigate this case ELSE !!! Basal melting zEi = - q_i_b(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 = - zqfont_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0) zdeltah(ji,jk) = - zfmdt / rhoic ! Gross thickness change zqfont_bo(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * rhoic * ( - zdE ) ! Update heat available zdeltah(ji,jk) = MAX( zdeltah(ji,jk), - zh_i(ji) ) ! Update thickness change dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) ! Update basal melt zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step zQm = MAX ( zfmdt , 0.0 ) * zEw ! Heat exchanged with ocean zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * rhoic * zdE * r1_rdtice ! for heat conservation ENDIF ! Contribution to salt flux sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) & & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice ! Contribution to energy flux to the ocean (J/m2) rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji) ENDIF END DO ! ji END DO ! jk ! !------------------- IF( con_i .AND. jiindex_1d > 0 ) THEN ! Conservation test ! !------------------- DO ji = kideb, kiut IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0 ) THEN IF( ( zfbase(ji) + zdq_i(ji) ) >= 1.e-3 ) THEN numce_dh = numce_dh + 1 meance_dh = meance_dh + zfbase(ji) + zdq_i(ji) ENDIF IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN WRITE(numout,*) ' ALERTE heat loss for basal melt : ii, ij, jl :', ii, ij, jl WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) WRITE(numout,*) ' zfbase : ', zfbase(ji) WRITE(numout,*) ' zdq_i : ', zdq_i(ji) WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji) WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji) WRITE(numout,*) ' qlbbq_1d : ', qlbbq_1d(ji) WRITE(numout,*) ' s_i_new : ', s_i_new(ji) WRITE(numout,*) ' sss_m : ', sss_m(ii,ij) WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) WRITE(numout,*) ' innermelt : ', INT( zinnermelt(ji) ) ENDIF ENDIF END DO IF( numce_dh > 0 ) meance_dh = meance_dh / numce_dh WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d) ! ENDIF ! !------------------------------------------------------------------------------! ! 5) Pathological cases ! !------------------------------------------------------------------------------! ! !---------------------------------------------- ! 5.1 Excessive ablation in a 1-category model !---------------------------------------------- DO ji = kideb, kiut ! ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15) IF( jpl == 1 ) THEN ; zdhbf = MAX( hmelt , dh_i_bott(ji) ) ELSE ; zdhbf = dh_i_bott(ji) ENDIF zdvres = zdhbf - dh_i_bott(ji) dh_i_bott(ji) = zdhbf sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice ! ! excessive energy is sent to lateral ablation zinda = MAX( 0._wp, SIGN( 1._wp , 1.0 - at_i_b(ji) - epsi10 ) ) fsup(ji) = zinda * rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi10 ) * zdvres * r1_rdtice END DO !----------------------------------- ! 5.2 More than available ice melts !----------------------------------- ! then heat applied minus heat content at previous time step should equal heat remaining ! DO ji = kideb, kiut ! Adapt the remaining energy if too much ice melts !-------------------------------------------------- zdvres = MAX( 0._wp, - ht_i_b(ji) - dh_i_surf(ji) - dh_i_bott(ji) ) zdvsur = MIN( 0._wp, dh_i_surf(ji) + zdvres ) - dh_i_surf(ji) ! fill the surface first zdvbot = MAX( 0._wp, zdvres - zdvsur ) ! then the bottom dh_i_surf (ji) = dh_i_surf(ji) + zdvsur ! clem dh_i_bott (ji) = dh_i_bott(ji) + zdvbot ! clem ! new ice thickness (clem) zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice zhgnew(ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 ! !since ice volume is only used for outputs, we keep it global for all categories dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) ! remaining heat zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) + zqfont_bo(ji) ) ! If snow remains, energy is used to melt snow zhni = ht_s_b(ji) ! snow depth at previous time step zihg = MAX( zzero , SIGN ( zone , - ht_s_b(ji) ) ) ! =0 if snow ! energy of melting of remaining snow zinda = MAX( 0._wp, SIGN( 1._wp , zhni - epsi10 ) ) zqt_s(ji) = ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi10 ) * zinda zdhnm = - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 ) zhnfi = zhni + zdhnm zfdt_final(ji) = MAX( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 ) ht_s_b(ji) = MAX( zzero , zhnfi ) zqt_s(ji) = zqt_s(ji) * ht_s_b(ji) ! we recompute dh_s_tot (clem) dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) ! Mass variations of ice and snow !--------------------------------- ! ! mass variation of the jl category zzfmass_s = - a_i_b(ji) * ( zhni - ht_s_b(ji) ) * rhosn ! snow zzfmass_i = a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic ! ice ! zfmass_i(ji) = zzfmass_i ! ice variation saved to compute salt flux (see below) ! ! ! mass variation cumulated over category !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s ! snow !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i ! ice ! Remaining heat to the ocean !--------------------------------- focea(ji) = - zfdt_final(ji) * r1_rdtice ! focea is in W.m-2 * dt ! residual salt flux (clem) !-------------------------- ! surface sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvsur * rhoic * r1_rdtice ! bottom IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) >= 0._wp ) THEN ! melting sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice ELSE ! growth sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice ENDIF ! ! diagnostic ii = MOD( npb(ji) - 1, jpi ) + 1 ij = ( npb(ji) - 1 ) / jpi + 1 diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice END DO ftotal_fin (:) = zfdt_final(:) * r1_rdtice !--------------------------- ! heat fluxes !--------------------------- DO ji = kideb, kiut zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) ! =1 if ice ! ! Heat flux ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1 ! excessive total ablation energy (focea) sent to the ocean qfvbq_1d(ji) = qfvbq_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice zihic = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) ) ! equals 0 if ht_i = 0, 1 if ht_i gt 0 fscbq_1d(ji) = a_i_b(ji) * fstbif_1d(ji) qldif_1d(ji) = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea (ji) * a_i_b(ji) * rdt_ice & & + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice END DO ! ji !------------------------------------------- ! Correct temperature, energy and thickness !------------------------------------------- DO ji = kideb, kiut zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) t_su_b(ji) = zihgnew * t_su_b(ji) + ( 1.0 - zihgnew ) * rtt END DO ! ji DO jk = 1, nlay_i DO ji = kideb, kiut zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) t_i_b(ji,jk) = zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt q_i_b(ji,jk) = zihgnew * q_i_b(ji,jk) END DO END DO ! ji DO ji = kideb, kiut ht_i_b(ji) = zhgnew(ji) END DO ! ji ! !------------------------------------------------------------------------------| ! 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( zzero , ( rhosn * ht_s_b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic ) ) zhgnew(ji) = MAX( zhgnew(ji) , zhgnew(ji) + dh_snowice(ji) ) zhnnew = MIN( ht_s_b(ji) , ht_s_b(ji) - dh_snowice(ji) ) ! Changes in ice volume and ice mass. dvnbq_1d (ji) = a_i_b(ji) * ( zhgnew(ji)-ht_i_b(ji) ) dmgwi_1d (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew - ht_s_b(ji) ) * rhosn ! 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_b(ji) ! entrapment during snow ice formation ! clem: new salinity difference stored (to be used in limthd_ent.F90) IF ( num_sal == 2 ) THEN i_ice_switch = MAX( 0._wp , SIGN( 1._wp , zhgnew(ji) - epsi10 ) ) ! salinity dif due to snow-ice formation dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch ! salinity dif due to bottom growth IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0._wp ) THEN dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch ENDIF ENDIF ! Actualize new snow and ice thickness. ht_s_b(ji) = zhnnew ht_i_b(ji) = zhgnew(ji) ! Total ablation ! new lines added to debug IF( ht_i_b(ji) <= 0._wp ) a_i_b(ji) = 0._wp ! diagnostic ( snow ice growth ) ii = MOD( npb(ji) - 1, jpi ) + 1 !MV HC 2014 useless ij = ( npb(ji) - 1 ) / jpi + 1 ! MV HC 2014 useless diag_sni_gr(ii,ij) = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice ! Contribution to energy flux to the ocean [J/m2], >0 zfmdt = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0.0 ) ! <0 zsstK = sst_m(ii,ij) + rt0 zEw = rcp * ( zsstK - rt0 ) zQm = zfmdt * zEw rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji) ! Contribution to salt flux sfx_thd_1d(ji) = sfx_thd_1d(ji) + sss_m(ii,ij) * a_i_b(ji) * zfmdt * r1_rdtice ! Contribution to mass fluxes rdm_ice_1d(ji) = rdm_ice_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic rdm_snw_1d(ji) = rdm_snw_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn END DO !ji ! CALL wrk_dealloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i ) CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay ) ! CALL wrk_dealloc( jpij, zviold, zvsold ) ! clem ! 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