MODULE limthd_dh #if defined key_lim3 !!====================================================================== !! *** MODULE limthd_dh *** !! thermodynamic growth and decay of the ice !!====================================================================== !!---------------------------------------------------------------------- !! lim_thd_dh : vertical accr./abl. and lateral ablation of sea ice !! * Modules used USE par_oce ! ocean parameters USE phycst ! physical constants (OCE directory) USE ice_oce ! ice variables USE thd_ice USE iceini USE limistate USE in_out_manager USE ice USE par_ice USE limicepoints IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC lim_thd_dh ! called by lim_thd !! * Module variables REAL(wp) :: & ! constant values epsi20 = 1e-20 , & epsi13 = 1e-13 , & epsi16 = 1e-16 , & zzero = 0.e0 , & zone = 1.e0 !!---------------------------------------------------------------------- !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2005) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_thd_dh(kideb,kiut,jl) !!------------------------------------------------------------------ !! *** ROUTINE lim_thd_dh *** !!------------------------------------------------------------------ !! ** Purpose : !! This routine 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 !! ** Steps !! 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 !! !! ** Arguments !! !! ** Inputs / Outputs !! !! ** External !! !! ** References : Bitz and Lipscomb, JGR 99 !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 !! Vancoppenolle Fichefet and Bitz, GRL 2005 !! !! ** History : !! original code 01-04 (LIM) !! original routine !! (05-2003) Martin Vancoppenolle, Louvain-La-Neuve, Belgium !! (06/07-2005) Martin Vancoppenolle for the 3D version !! !!------------------------------------------------------------------ !! * Arguments INTEGER , INTENT (IN) :: & kideb , & !: Start point on which the the computation is applied kiut , & !: End point on which the the computation is applied jl !: Thickness cateogry number !! * Local variables INTEGER :: & ji , & !: space index jk , & !: ice layer index isnow , & !: switch for presence (1) or absence (0) of snow index , & !: ice point index (limicepoints) zji , & !: 2D corresponding indices to ji zjj , & isnowic , & !: snow ice formation not i_ice_switch , & !: ice thickness above a certain treshold or not iter REAL(wp) :: & zhsnew , & !: new snow thickness zihgnew , & !: switch for total ablation ztmelts , & !: melting point zhn , & zdhcf , & zdhbf , & zhni , & zhnfi , & zihg , & zdhgm , & zdhnm , & zhnnew , & zdfseqv , & zeps = 1.0e-13, & zhisn , & zfracs , & !: fractionation coefficient for bottom salt !: entrapment zds , & !: increment of bottom ice salinity zoldsinew , & !: dummy argument for error diagnosis zcoeff , & !: dummy argument for snowfall partitioning !: over ice and leads zjiindex , & !: zdhs_temp, & zsm_snowice, & !: snow-ice salinity zswi1 , & !: switch for computation of bottom salinity zswi12 , & !: switch for computation of bottom salinity zswi2 , & !: switch for computation of bottom salinity zgrr , & !: bottom growth rate zihic , & !: ztform !: bottom formation temperature REAL(wp) , DIMENSION(jpij) :: & zh_i , & ! ice layer thickness zh_s , & ! snow layer thickness ztfs , & ! melting point zhsold , & ! old snow thickness zqprec , & !: energy of fallen snow zqfont_su , & ! incoming, remaining surface energy zqfont_bo ! incoming, bottom energy REAL(wp) , DIMENSION(jpij) :: & z_f_surf, & ! surface heat for ablation zhgnew ! new ice thickness REAL(wp), DIMENSION(jpij) :: & zdh_s_mel , & ! snow melt zdh_s_pre , & ! snow precipitation zdh_s_sub , & ! snow sublimation zfsalt_melt ! salt flux due to ice melt REAL(wp) , DIMENSION(jpij,jkmax) :: & zdeltah ! Pathological cases REAL(wp), DIMENSION(jpij) :: & zfdt_init , & !: total incoming heat for ice melt zfdt_final , & !: total remaing heat for ice melt zqt_i , & !: total ice heat content zqt_s , & !: total snow heat content zqt_dummy !: dummy heat content REAL(wp), DIMENSION(jpij,jkmax) :: & zqt_i_lay , & !: total ice heat content zqt_s_lay !: total snow heat content ! Heat conservation REAL(wp), DIMENSION(jpij) :: & zfbase, & zdq_i INTEGER, DIMENSION(jpij) :: & innermelt REAL(wp) :: & meance_dh INTEGER :: & num_iter_max, & numce_dh !!----------------------------------------------------------------------------- WRITE(numout,*) 'lim_thd_dh : computation of vertical snow/ice accretion/ablation' WRITE(numout,*) '~~~~~~~~~' zfsalt_melt(:) = 0.0 ftotal_fin(:) = 0.0 zfdt_init(:) = 0.0 zfdt_final(:) = 0.0 DO ji = kideb, kiut old_ht_i_b(ji) = ht_i_b(ji) old_ht_s_b(ji) = ht_s_b(ji) 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.0 zqfont_bo(:) = 0.0 dsm_i_se_1d(:) = 0.0 dsm_i_si_1d(:) = 0.0 ! ! heat conservation ! IF ( con_i ) THEN ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' zf_surf : ', z_f_surf(jiindex_1d) ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' qnsr_ice : ', qnsr_ice_1d(jiindex_1d) ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' qsr_ice : ', qsr_ice_1d(jiindex_1d) ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' fc_su : ', fc_su(jiindex_1d) ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' i0 : ', i0(jiindex_1d) ! ENDIF !------------------------------------------------------------------------------! ! 2) Computing layer thicknesses and snow and sea-ice enthalpies. ! !------------------------------------------------------------------------------! !----------------- ! Layer thickness !----------------- DO ji = kideb,kiut zh_i(ji) = ht_i_b(ji) / nlay_i zh_s(ji) = ht_s_b(ji) / nlay_s END DO !-------------------------------------- ! Snow energy of melting, heat content !-------------------------------------- zqt_s(:) = 0.0 DO jk = 1, nlay_s DO ji = kideb,kiut ! this line could be removed ! q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) ! heat conservation zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s END DO END DO ! ! heat conservation ! qt_s_in(:,jl) = zqt_s(:) ! IF (jiindex_1d.GT.0) & ! WRITE(numout,*) 'jl : ', jl, ' zqts : ', zqt_s(jiindex_1d) / & ! rdt_ice ! IF (jiindex_1d.GT.0) THEN ! WRITE(numout,*) ' Vertical profile : ' ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d) ! WRITE(numout,*) ' qm0 : ', q_s_b(jiindex_1d,1:nlay_s) * ht_s_b(jiindex_1d) / & ! rdt_ice ! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1) ! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1) ! WRITE(numout,*) ' zthick0: ', ht_s_b(jiindex_1d) ! ENDIF !------------------------------------- ! Ice energy of melting, heat content !------------------------------------- zqt_i(:) = 0.0 DO jk = 1, nlay_i DO ji = kideb,kiut ! Melting point in K ! ztmelts = - tmut * s_i_b(ji,jk) + rtt ! this line could be removed ! q_i_b(ji,jk) = rhoic * & ! ( cpic * ( ztmelts - t_i_b(ji,jk) ) & ! + lfus * ( 1.0 - (ztmelts-rtt) / & ! MIN( ( t_i_b(ji,jk) - rtt ) , - zeps ) ) & ! - rcp * ( ztmelts-rtt ) ) zqt_i(ji) = zqt_i(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i zqt_i_lay(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i END DO END DO ! IF (jiindex_1d.GT.0) & ! WRITE(numout,*) 'jl : ', jl, ' zqti : ', zqt_i(jiindex_1d) / & ! rdt_ice ! IF (jiindex_1d.GT.0) THEN ! WRITE(numout,*) ' --- Vertical profile : ' ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d) ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d) ! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d, 1:nlay_i) ! WRITE(numout,*) ' t_i_b : ', t_i_b(jiindex_1d, 1:nlay_i) ! WRITE(numout,*) ' qm0 : ', zqt_i_lay(jiindex_1d,1:nlay_i) / rdt_ice ! WRITE(numout,*) ' zthick0: ', ht_i_b(jiindex_1d) / nlay_i ! ENDIF ! !------------------------------------------------------------------------------| ! 3) Surface ablation and sublimation | !------------------------------------------------------------------------------| ! !--------------------- ! 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.0 ! Melt of fallen snow DO ji = kideb, kiut ! tatm_ice is now in K zqprec(ji) = rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus ) zqfont_su(ji) = z_f_surf(ji) * rdt_ice zdeltah(ji,1) = MIN( 0.0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) ) zqfont_su(ji) = MAX( 0.0 , - 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) END DO ! IF (jiindex_1d .GT. 0) THEN ! WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d) ! WRITE(numout,*) ' zqprec : ', zqprec(jiindex_1d)*( zdh_s_pre(jiindex_1d) + zdh_s_mel(jiindex_1d) ) / rdt_ice ! WRITE(numout,*) ' tatm : ', tatm_ice_1d(jiindex_1d) ! WRITE(numout,*) 'jl : ', jl, ' zqts : ', zqt_s(jiindex_1d) / & ! rdt_ice ! ENDIF ! updated total snow heat content zqt_s(ji) = MAX ( zqt_s(ji) - zqfont_su(ji) , 0.0 ) ! IF (jiindex_1D .GT. 0) & ! WRITE(numout,*) 'jl : ', jl, ' zqts : ', zqt_s(jiindex_1d) / & ! rdt_ice ! 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) ! IF (ji.eq.jiindex_1d) WRITE(numout,*) ' zqfont_su : ', & ! zqfont_su(jiindex_1d) 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 ! !+++++ ! IF (jiindex_1d .GT. 0 ) THEN ! WRITE(numout,*) ' dhs_pre :', zdh_s_pre(jiindex_1d) ! WRITE(numout,*) ' dhs_mel :', zdh_s_mel(jiindex_1d) ! WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d) ! ENDIF ! !+++++ ! 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 ) ! Volume and mass variations of snow dvsbq_1d(ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) & - zdh_s_mel(ji) ) dvsbq_1d(ji) = MIN( zzero, dvsbq_1d(ji) ) rdmsnif_1d(ji) = rhosn*dvsbq_1d(ji) END DO ! ji ! !+++++ ! IF (jiindex_1d .GT. 1) WRITE(numout,*) ' dhs_tot :', dh_s_tot(jiindex_1d) ! IF (jiindex_1d .GT. 1) WRITE(numout,*) ' zqfont_su :', zqfont_su(jiindex_1d) ! !+++++ !---------------------- ! Surface ice ablation !---------------------- DO ji = kideb, kiut ! IF (ji.eq.jiindex_1d) WRITE(numout,*) ' zqfont_su : ', & ! zqfont_su(jiindex_1d) dh_i_surf(ji) = 0.0 ! Heat conservation test z_f_surf(ji) = zqfont_su(ji) / rdt_ice ! heat conservation test zdq_i(ji) = 0.0 END DO ! ji DO jk = 1, nlay_i DO ji = kideb, kiut zdeltah(ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk) zqfont_su(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * & q_i_b(ji,jk) zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) ) dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * & q_i_b(ji,jk) / rdt_ice ! Contribution to salt flux -> to change zji = MOD( npb(ji) - 1, jpi ) + 1 zjj = ( npb(ji) - 1 ) / jpi + 1 zfsalt_melt(ji) = zfsalt_melt(ji) + & ! ( sss_io(zji,zjj) - s_i_b(ji,jk) ) * & ( sss_io(zji,zjj) - sm_i_b(ji) ) * & a_i_b(ji) * & MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice !+++++++++++ ! IF ( (zji.eq.jiindex) .AND. (zjj.eq.jjindex) ) THEN ! WRITE(numout,*) ' surf abl ' ! WRITE(numout,*) ' zfsalt_melt : ', zfsalt_melt(ji) / ( sss_io(zji,zjj)+epsi16) ! WRITE(numout,*) ' sss_io : ', sss_io(zji,zjj) ! ENDIF !+++++++++++ END DO ! ji END DO ! jk !------------------- ! Conservation test !------------------- IF ( con_i ) THEN numce_dh = 0 meance_dh = 0.0 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!.GE. 1.0e-3 ) .AND. & ! ( ( ht_i_b(ji) + dh_i_bott(ji) ) .GE. 1.0e-6 ) ) THEN WRITE(numout,*) ' ALERTE heat loss for surface melt ' WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, 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_io : ', sss_io(zji,zjj) ENDIF END DO ! ji IF ( numce_dh .GT. 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 ! con_i !------------------ ! Snow sublimation !------------------ ! !++++++ ! zqt_dummy(:) = 0.0 ! DO jk = 1, nlay_s ! DO ji = kideb,kiut ! q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) ! ! heat conservation ! zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * old_ht_s_b(ji) / nlay_s ! END DO ! END DO ! IF (jiindex_1d .GT. 0) & ! WRITE(numout,*) 'jl : ', jl, ' Old zqts : ', zqt_dummy(jiindex_1d) / & ! rdt_ice !++++++ ! !++++++ ! zqt_dummy(:) = 0.0 ! DO jk = 1, nlay_s ! DO ji = kideb,kiut ! q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) ! ! heat conservation ! zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s ! END DO ! END DO ! IF (jiindex_1d .GT. 0) & ! WRITE(numout,*) 'jl : ', jl, ' New zqts : ', zqt_dummy(jiindex_1d) / & ! rdt_ice ! !++++++ DO ji = kideb, kiut ! if qla is positive (upwards), heat goes to the atmosphere, therefore ! snow sublimates, if qla is negative (downwards), snow condensates zdh_s_sub(ji) = - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice 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 ! which may be different in case of total snow melt dh_s_tot(ji) = ht_s_b(ji) - zhsold(ji) ! IF ( ht_s_b(ji) .LE. 0.0 ) THEN ! dh_s_tot(ji) = MAX( 0.0 , dh_s_tot(ji) ) ! ENDIF qt_s_in(ji,jl) = qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1) END DO !ji ! IF ( jiindex_1d .GT. 0 ) THEN ! WRITE(numout,*) ' dh_s_sub : ', zdh_s_sub(jiindex_1d) ! WRITE(numout,*) ' dhs_tot : ', dh_s_tot(jiindex_1d) ! ENDIF !++++++ zqt_dummy(:) = 0.0 DO jk = 1, nlay_s DO ji = kideb,kiut q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) ! heat conservation zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s END DO END DO ! IF (jiindex_1d .GT. 0) & ! WRITE(numout,*) 'jl : ', jl, ' Old zqts : ', zqt_dummy(jiindex_1d) / & ! rdt_ice ! !++++++ ! IF (jiindex_1d .GT. 0) THEN ! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1) ! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1) ! ENDIF DO jk = 1, nlay_s !n DO ji = kideb, kiut !n ! 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 ! IF (jiindex_1d .GT. 0) THEN ! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1) ! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1) ! ENDIF ! !------------------------------------------------------------------------------! ! 4) Basal growth / melt ! !------------------------------------------------------------------------------! ! ! Ice basal growth / melt is given by the ratio of heat budget over basal ! ice heat content. Basal heat budget is given by the difference 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 !--------------------------------------------- ! Basal growth - salinity not varying in time !--------------------------------------------- IF ( num_sal .NE. 2 ) THEN DO ji = kideb, kiut IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN s_i_new(ji) = sm_i_b(ji) ! Melting point in K ztmelts = - tmut * s_i_new(ji) + rtt ! New ice heat content (Bitz and Lipscomb, 1999) ztform = t_i_b(ji,nlay_i) ! t_bo_b crashes in the ! Baltic q_i_b(ji,nlay_i+1) = rhoic * & ( cpic * ( ztmelts - ztform ) & + lfus * ( 1.0 - ( ztmelts - rtt ) / & ( ztform - rtt ) ) & - rcp * ( ztmelts-rtt ) ) ! Basal growth rate = - F*dt / q dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + & qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) ! IF ( ji .EQ. jiindex_1d ) THEN ! WRITE(numout,*) ' s_i_new : ', s_i_new(ji) ! WRITE(numout,*) ' ztmelts : ', ztmelts ! WRITE(numout,*) ' t_bo : ', t_bo_b(ji) ! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,nlay_i+1) ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) ! ENDIF ENDIF ! heat budget END DO ! ji ENDIF ! num_sal ! IF ( jiindex_1d .GT. 0 ) THEN ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d) ! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d,nlay_i+1) ! WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(jiindex_1d) ! WRITE(numout,*) ' fbif_1d : ', fbif_1d(jiindex_1d) ! WRITE(numout,*) ' qlbbq_1d : ', qlbbq_1d(jiindex_1d) ! ENDIF !----------------------------------------- ! Basal growth - salinity varying in time !----------------------------------------- IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! the growth rate (dh_i_bott) is function of the new ice ! heat content (q_i_b(nlay_i+1)). q_i_b depends on the new ice ! salinity (snewice). snewice depends on dh_i_bott ! it converges quickly, so, no problem ! Initial value (tested 1D, can be anything between 1 and 20) num_iter_max = 4 s_i_new(:) = 4.0 ! Iterative procedure DO iter = 1, num_iter_max DO ji = kideb, kiut IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN zji = MOD( npb(ji) - 1, jpi ) + 1 zjj = ( npb(ji) - 1 ) / jpi + 1 ! Melting point in K ztmelts = - tmut * s_i_new(ji) + rtt ! New ice heat content (Bitz and Lipscomb, 1999) q_i_b(ji,nlay_i+1) = rhoic * & ( cpic * ( ztmelts - t_bo_b(ji) ) & + lfus * ( 1.0 - ( ztmelts - rtt ) / & ( t_bo_b(ji) - rtt ) ) & - rcp * ( ztmelts-rtt ) ) ! !+++++++ ! IF (ji.EQ.jiindex_1d) THEN ! WRITE(numout,*) ' ztmelts : ', ztmelts ! WRITE(numout,*) ' t_bo_b : ', t_bo_b(ji) ! WRITE(numout,*) ' q_i_b(ji,nlay_i+1) : ', q_i_b(ji,nlay_i+1) ! ENDIF ! !+++++++ ! Bottom growth rate = - F*dt / q dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) & + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) ! New ice salinity ( Cox and Weeks, JGR, 1988 ) ! zswi2 (1) if dh_i_bott/rdt .GT. 3.6e-7 ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8 ! zswi1 (1) if dh_i_bott/rdt .LT. 2.0e-8 ! zgrr = MAX( dh_i_bott(ji) / rdt_ice, 0.0 ) zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , zeps ) ) 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 ! IF (ji.EQ.jiindex_1d) THEN ! WRITE(numout,*) ' zgrr : ', zgrr ! WRITE(numout,*) ' zswi2 : ', zswi2, ' zswi12 :', zswi12, ' zswi1 :', zswi1 ! ENDIF zfracs = zswi1 * 0.12 + & zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) + & zswi2 * 0.26 / & ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) zds = zfracs*sss_io(zji,zjj) - s_i_new(ji) s_i_new(ji) = zfracs * sss_io(zji,zjj) ENDIF ! fc_bo_i END DO ! ji END DO ! iter ! Final values DO ji = kideb, kiut IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN ! New ice salinity must not exceed 15 psu zoldsinew = s_i_new(ji) s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) ! Metling point in K ztmelts = - tmut * s_i_new(ji) + rtt ! New ice heat content (Bitz and Lipscomb, 1999) q_i_b(ji,nlay_i+1) = rhoic * & ( cpic * ( ztmelts - t_bo_b(ji) ) & + lfus * ( 1.0 - ( ztmelts - rtt ) / & ( t_bo_b(ji) - rtt ) ) & - rcp * ( ztmelts-rtt ) ) ! Basal growth rate = - F*dt / q ! sometimes, growth is very high because of very high bottom cond ! flux... this is dangerous dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + & qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) ! new lines !----------------- ! Salinity update !----------------- ! entrapment during bottom growth dsm_i_se_1d(ji) = ( s_i_new(ji)*dh_i_bott(ji) + & sm_i_b(ji)*ht_i_b(ji) ) / & MAX( ht_i_b(ji) + dh_i_bott(ji) ,zeps ) & - sm_i_b(ji) ENDIF ! heat budget END DO ! ji ENDIF ! num_sal ! IF ( jiindex_1d .GT. 0 ) THEN ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d) ! WRITE(numout,*) ' s_i_new : ', s_i_new(jiindex_1d) ! ENDIF !------------ ! Basal melt !------------ !++++++ meance_dh = 0.0 numce_dh = 0 innermelt(:) = 0 !++++++ DO ji = kideb, kiut ! heat convergence at the surface > 0 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN s_i_new(ji) = s_i_b(ji,nlay_i) zqfont_bo(ji) = rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) zfbase(ji) = zqfont_bo(ji) / rdt_ice ! heat conservation test zdq_i(ji) = 0.0 dh_i_bott(ji) = 0.0 ENDIF END DO ! IF ( jiindex_1d .GT. 0 ) THEN ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice ! WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(jiindex_1d) ! WRITE(numout,*) ' fbif : ', fbif_1d(jiindex_1d) ! WRITE(numout,*) ' qlbbq : ', qlbbq_1d(jiindex_1d) ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d) ! ENDIF DO jk = nlay_i, 1, -1 DO ji = kideb, kiut IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN ztmelts = - tmut * s_i_b(ji,jk) + rtt ! sometimes internal temperature gt melting point ! the whole layer is removed ! IF (ji.EQ.jiindex_1d) THEN ! WRITE(numout,*) ' jk : ', jk ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(ji) ! WRITE(numout,*) ' t_i_b : ', t_i_b(ji,jk) ! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,jk) ! ENDIF IF ( t_i_b(ji,jk) .GE. ztmelts ) THEN zdeltah(ji,jk) = - zh_i(ji) ! zqfont_bo(ji) = zqfont_bo(ji) dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) !++++++ innermelt(ji) = 1 !++++++ ! IF (ji.EQ.jiindex_1d) THEN ! WRITE(numout,*) ' Inner melt: ', innermelt(ji) ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) ! ENDIF ELSE ! normal ablation zdeltah(ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk) zqfont_bo(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * & q_i_b(ji,jk) zdeltah(ji,jk) = MAX(zdeltah(ji,jk), - zh_i(ji) ) dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * & q_i_b(ji,jk) / rdt_ice ! contribution to salt flux zji = MOD( npb(ji) - 1, jpi ) + 1 zjj = ( npb(ji) - 1 ) / jpi + 1 zfsalt_melt(ji) = zfsalt_melt(ji) + & ! ( sss_io(zji,zjj) - s_i_b(ji,jk) ) * & ( sss_io(zji,zjj) - sm_i_b(ji) ) * & a_i_b(ji) * & MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice ! IF (ji.EQ.jiindex_1d) THEN ! WRITE(numout,*) ' No inner melt ' ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(ji) ! WRITE(numout,*) ' t_i_b : ', t_i_b(ji,jk) ! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,jk) ! WRITE(numout,*) ! ENDIF ENDIF ! IF ( ji .EQ. jiindex_1d ) THEN ! WRITE(numout,*) ' basal melt ' ! WRITE(numout,*) ' sss_io : ', sss_io(zji,zjj) ! WRITE(numout,*) ' zfsalt_melt : ', zfsalt_melt(ji) ! WRITE(numout,*) ' zfsalt_melt good units : ', zfsalt_melt(ji) / ( sss_io(zji,zjj) + epsi16 ) ! ENDIF !+++++ ENDIF END DO ! ji END DO ! jk IF ( con_i ) THEN DO ji = kideb, kiut IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN !------------------- ! Conservation test !------------------- IF ( ( zfbase(ji) + zdq_i(ji) ) .GE. 1.0e-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!.GE. 1.0e-3 ) .AND. & ! ( ( ht_i_b(ji) + dh_i_bott(ji) ) .GE. 1.0e-6 ) ) THEN WRITE(numout,*) ' ALERTE heat loss for basal melt ' WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, 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_io : ', sss_io(zji,zjj) WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) WRITE(numout,*) ' innermelt : ', innermelt(ji) ENDIF ENDIF ! heat convergence at the surface END DO ! ji IF ( numce_dh .GT. 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 ! con_i ! IF ( jiindex_1d .GT. 0 ) THEN ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d) ! ENDIF ! !------------------------------------------------------------------------------! ! 5) Pathological cases ! !------------------------------------------------------------------------------! ! !------------------------------------------ ! 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) zdhbf = dh_i_bott(ji) IF (jpl.EQ.1) zdhbf = MAX( hmelt , dh_i_bott(ji) ) ! excessive energy is sent to lateral ablation fsup(ji) = rhoic*lfus * at_i_b(ji) / MAX( ( 1.0 - at_i_b(ji) ),epsi13) & * ( zdhbf - dh_i_bott(ji) ) / rdt_ice dh_i_bott(ji) = zdhbf !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) !new ice thickness zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) ! diagnostic ( bottom ice growth ) zji = MOD( npb(ji) - 1, jpi ) + 1 zjj = ( npb(ji) - 1 ) / jpi + 1 diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) & / rdt_ice diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) & / rdt_ice diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) & / rdt_ice END DO ! IF (jiindex_1d .GT. 0 ) THEN ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d) ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d) ! ENDIF !-------------------------- ! If too much 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 !-------------------------------------------------- zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice ! 0 if no more ice zhgnew(ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 ! remaining heat zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) + zqfont_bo(ji) ) ! & ! + zihgnew * MAX ( zfdt_init(ji) - zqt_i(ji) - zqt_s(ji), 0.0 ) ! !++++++++++ ! IF (ji.EQ.jiindex_1d) THEN ! WRITE(numout,*) ' zfdt_init : ', zfdt_init(ji) / rdt_ice ! WRITE(numout,*) ' zqt_i : ', zqt_i(ji) / rdt_ice ! WRITE(numout,*) ' 1. zqt_s : ', zqt_s(ji) / rdt_ice ! WRITE(numout,*) ' zqt : ', ( zqt_i(ji) + zqt_s(ji) ) / rdt_ice ! WRITE(numout,*) ' zhgnew : ', zhgnew(ji) ! WRITE(numout,*) ' zihgnew : ', zihgnew ! WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji) ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) ! WRITE(numout,*) ' 1. zfdt_final: ', zfdt_final(ji) / rdt_ice ! WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d) / rdt_ice ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice ! ENDIF ! !++++++++++ ! 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 zqt_s(ji) = ( 1. - zihg) * zqt_s(ji) / MAX( zhni, zeps ) zdhnm = - ( 1. - zihg ) * ( 1. - zihgnew ) * ( zfdt_final(ji) / & MAX( zqt_s(ji) , zeps ) ) 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) ! IF (ji.EQ.jiindex_1d) THEN ! WRITE(numout,*) ' 2. zfdt_final : ', zfdt_final(ji) / rdt_ice ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) ! WRITE(numout,*) ' zhgnew : ', zhgnew(ji) ! WRITE(numout,*) ' 2. zqt_s : ', zqt_s(ji) / rdt_ice ! WRITE(numout,*) ' zdhnm : ', zdhnm ! WRITE(numout,*) ! ENDIF ! Mass variations of ice and snow !--------------------------------- rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) * & (zhgnew(ji)-ht_i_b(ji))*rhoic ! good rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) * & (ht_s_b(ji)-zhni)*rhosn ! good too ! Remaining heat to the ocean !--------------------------------- ! check the switches here ! focea is in W.m-2 * dt focea(ji) = - zfdt_final(ji) / rdt_ice ! IF ( (zji.eq.jiindex) .AND. (zjj.eq.jjindex) ) THEN ! WRITE(numout,*) ' focea : ', focea(ji) ! ENDIF END DO ! IF (jiindex_1d .GT. 0 ) THEN ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d) ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d) ! ENDIF ftotal_fin (:) = zfdt_final(:) / rdt_ice ! IF (jiindex_1d .GT. 1 ) WRITE(numout,*) ' ftotal_fin : ', ftotal_fin(jiindex_1d) !--------------------------- ! Salt flux and heat fluxes !--------------------------- DO ji = kideb, kiut zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice ! Salt flux zji = MOD( npb(ji) - 1, jpi ) + 1 zjj = ( npb(ji) - 1 ) / jpi + 1 IF ( num_sal .NE. 4 ) & fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) + & (1.0 - zihgnew) * rdmicif_1d(ji) * & ( sss_io(zji,zjj) - sm_i_b(ji) ) / rdt_ice ! new lines IF ( num_sal .EQ. 4 ) & fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) + & (1.0 - zihgnew) * rdmicif_1d(ji) * & ( sss_io(zji,zjj) - bulk_sal ) / rdt_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 ! ! astarojna zihic = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) ) ! equals 0 if ht_i = 0, 1 if ht_i gt 0 ! fstbif_1d est cumulatif merde fscbq_1d(ji) = a_i_b(ji) * fstbif_1d(ji) ! here there is a bug 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) ! IF (ji.eq.jiindex_1d) THEN ! WRITE(numout,*) ' jk : ', jk ! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,jk) ! WRITE(numout,*) ' t_i_b : ', t_i_b(jiindex_1d, 1:nlay_i) ! WRITE(numout,*) ' zihgnew : ', zihgnew, ' zhgnew : ', & ! zhgnew(ji) ! ENDIF END DO END DO ! ji ! IF (jiindex_1d.GT.0) THEN ! WRITE(numout,*) ' --- Vertical profile : ' ! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d, 1:nlay_i) ! ENDIF DO ji = kideb, kiut ht_i_b(ji) = zhgnew(ji) END DO ! ji ! IF (jiindex_1d .GT. 0 ) THEN ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d) ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d) ! ENDIF ! !------------------------------------------------------------------------------| ! 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 #if defined key_lim_fdd !(presently Activated) rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) & * ( zhgnew(ji) - ht_i_b(ji) )*rhoic rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) & * ( zhnnew - ht_s_b(ji) )*rhosn ! Equivalent salt flux (1) Snow-ice formation component ! ----------------------------------------------------- zji = MOD( npb(ji) - 1, jpi ) + 1 zjj = ( npb(ji) - 1 ) / jpi + 1 ! zsm_snowice = MIN ( s_i_max, ( rhoic - rhosn ) / rhoic * & ! sss_io(zji,zjj) ) zsm_snowice = ( rhoic - rhosn ) / rhoic * & sss_io(zji,zjj) IF ( num_sal .NE. 2 ) zsm_snowice = sm_i_b(ji) IF ( num_sal .NE. 4 ) & fseqv_1d(ji) = fseqv_1d(ji) + & ( sss_io(zji,zjj) - zsm_snowice ) * & a_i_b(ji) * & ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice ! new lines IF ( num_sal .EQ. 4 ) & fseqv_1d(ji) = fseqv_1d(ji) + & ( sss_io(zji,zjj) - bulk_sal ) * & a_i_b(ji) * & ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice ! entrapment during snow ice formation i_ice_switch = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i_b(ji) + 1.0d-6 ) ) isnowic = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - dh_snowice(ji) ) ) * & i_ice_switch IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & dsm_i_si_1d(ji) = ( zsm_snowice*dh_snowice(ji) & + sm_i_b(ji) * ht_i_b(ji) & / MAX( ht_i_b(ji) + dh_snowice(ji), zeps) & - sm_i_b(ji) ) * isnowic #else rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) & * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic & + ( zhnnew - ht_s_b(ji) ) * rhosn ) #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).LE.0.0 ) a_i_b(ji) = 0.0 ! diagnostic ( snow ice growth ) zji = MOD( npb(ji) - 1, jpi ) + 1 zjj = ( npb(ji) - 1 ) / jpi + 1 diag_sni_gr(zji,zjj) = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / & rdt_ice END DO !ji ! IF (jiindex_1d .GT. 1 ) THEN ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d) ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d) ! ENDIF END SUBROUTINE lim_thd_dh #else !!====================================================================== !! *** MODULE limthd_dh *** !! no sea ice model !!====================================================================== CONTAINS SUBROUTINE lim_thd_dh ! Empty routine END SUBROUTINE lim_thd_dh #endif END MODULE limthd_dh