After a careful reading of the code, I think that SI3 can handle deposition.
If there is some, a corresponding amount of snow is added to the (0) snow layer where new snow is also added/removed from snowfall/sublimation. This new snow is added at the same temperature as the current snow, or 0degC if no snow. If snow is left in this (0) layer at the end of sbc timestep (after surface melt), the snow is moved to the “normal” snow layers ( layer (1) and above (if more than one layer is used) ).
I will test to see if what is given by our coupler corresponds to what is in the sublimation diagnostics.
Maybe I missed or forgot something, if anyone has any idea, please let me know =)
Summary
My comments after !==>DR:
Set in sbccpl from coupler values
zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) !==>DR: <0 if deposition
zevap_ice_total(:,:) = zevap_ice(:,:,1) !DR: Same
DO jl = 2, jpl
zevap_ice(:,:,jl) = zevap_ice(:,:,1) !DR just the same thing again
END DO
[...]
IF ( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN
! For conservative case zemp_ice has not been defined yet. Do it now.
zemp_ice(:,:) = zevap_ice_total(:,:) * picefr(:,:) - frcv(jpr_snow)%z3(:,:,1) * picefr(:,:) !==>DR: incuded in `emp_ice` but `emp_ice` does not do much in SI3, it seems to be used in CICE
ENDIF
! --- evaporation over ocean (used later for qemp) --- !
zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:) !==>DR: evap_tot = evap_oce * (1-a) + evap_ice *a, so if evap_ice <0 then evap_oce > evap_tot
evap_ice(:,:,:) = zevap_ice(:,:,:) !==>DR: temp. var. to global var.
[...]
IF( iom_use('subl_ai_cea') ) CALL iom_put( 'subl_ai_cea' , zevap_ice_total(:,:) * picefr(:,:) * tmask(:,:,1) ) ! Sublimation over sea-ice (cell average) !==>DR: Also deposition
IF( iom_use('evap_ao_cea') ) CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1) &
& - zevap_ice_total(:,:) * picefr(:,:) ) * tmask(:,:,1) ) ! ice-free oce evap (cell average)
IF ( iom_use('hflx_evap_cea') ) & ! heat flux from evap (cell average)
& CALL iom_put('hflx_evap_cea' , ( frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:) ) &
& * zcptn(:,:) * tmask(:,:,1) )
IF( iom_use('hflx_subl_cea') ) & ! heat flux from sublimation
& CALL iom_put('hflx_subl_cea' , SUM( qevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) * tmask(:,:,1) )
in ice_thd
CALL tab_2d_1d( npti, nptidx(1:npti), evap_ice_1d (1:npti), evap_ice(:,:,kl) ) !==>DR: From 2D to 1D
in ice_thd_dh
! Snow sublimation
!-----------------
! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
! comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
zdeltah (1:npti) = 0._wp ! total snow thickness that sublimates, < 0 !==>DR: > 0 if deposition
zevap_rema(1:npti) = 0._wp
DO ji = 1, npti ! horizontal grid in 1D
zdeltah (ji) = MAX( - evap_ice_1d(ji) * r1_rhos * rDt_ice, - h_s_1d(ji) ) ! amount of snw that sublimates, ==>DR: evap_ice_1d<0 so zdeltah>0
zevap_rema(ji) = evap_ice_1d(ji) * rDt_ice + zdeltah(ji) * rhos ! remaining evap in kg.m-2 (used for ice sublimation later on) ==>DR: because no limit up, zevap_rema = 0
END DO
DO jk = 0, nlay_s ! snow layers just one but 0 for the current snowfall 1 for the layer
DO ji = 1, npti ! horizontal grid in 1D
zdum = MAX( -zh_s(ji,jk), zdeltah(ji) ) ! snow layer thickness that sublimates, < 0 ==> DR: zh_s>0 so zh_s(ji,jk) <0 so zh_s(ji,jk) < zdeltah(ji) always so zdum = zdeltah(ji) > 0
!
hfx_sub_1d (ji) = hfx_sub_1d (ji) + ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! Heat flux of snw that sublimates [W.m-2], < 0 ==>DR: >0 for depostion,
wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux by sublimation ==>DR: <0 for depostion negative substraction means addition
! update thickness
h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdum ) ! ==>DR: Adding total thickness zdum >0
zh_s (ji,jk) = MAX( 0._wp , zh_s (ji,jk) + zdum ) ! ==>DR: Adding layer thickness zdum >0
!!$ IF( zh_s(ji,jk) == 0._wp ) ze_s(ji,jk) = 0._wp
! update sublimation left for the next layer
zdeltah(ji) = MIN( zdeltah(ji) - zdum, 0._wp ) !==>DR: must be zero for the deposition at the first layer, as zdum = zdeltah(ji)
END DO
END DO
[...]
! ! ============ !
! ! Ice !
! ! ============ !
[...]
DO jk = 1, nlay_i ! ice layers
DO ji = 1, npti ! horizontal grid in 1D
[...]
! Ice sublimation
! ---------------
zdum = MAX( - zh_i(ji,jk) , - zevap_rema(ji) * r1_rhoi ) ! ==>DR: is zero for deposition, what comes after is all zero
!
hfx_sub_1d(ji) = hfx_sub_1d(ji) + e_i_1d(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! Heat flux [W.m-2], < 0
wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux > 0
sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux >0
! clem: flux is sent to the ocean for simplicity
! but salt should remain in the ice except
! if all ice is melted. => must be corrected
! update remaining mass flux and thickness
zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi
zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum )
h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum )
dh_i_sub(ji) = dh_i_sub(ji) + zdum
! update heat content (J.m-2) and layer thickness
eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk)
h_i_old (ji,jk) = h_i_old (ji,jk) + zdum
END DO
END DO