Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
MODULE icethd_dh
!!======================================================================
!! *** MODULE icethd_dh ***
!! seaice : thermodynamic growth and melt
!!======================================================================
!! History : ! 2003-05 (M. Vancoppenolle) Original code in 1D
!! ! 2005-06 (M. Vancoppenolle) 3D version
!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
!!----------------------------------------------------------------------
#if defined key_si3
!!----------------------------------------------------------------------
!! 'key_si3' SI3 sea-ice model
!!----------------------------------------------------------------------
!! ice_thd_dh : vertical sea-ice growth and melt
!!----------------------------------------------------------------------
USE dom_oce ! ocean space and time domain
USE phycst ! physical constants
USE ice ! sea-ice: variables
USE ice1D ! sea-ice: thermodynamics variables
USE icethd_sal ! sea-ice: salinity profiles
USE icevar ! for CALL ice_var_snwblow
!
USE in_out_manager ! I/O manager
USE lib_mpp ! MPP library
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
IMPLICIT NONE
PRIVATE
PUBLIC ice_thd_dh ! called by ice_thd
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: icethd_dh.F90 14686 2021-04-08 15:36:01Z clem $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE ice_thd_dh
!!------------------------------------------------------------------
!! *** ROUTINE ice_thd_dh ***
!!
!! ** Purpose : compute ice and snow thickness changes due to growth/melting
!!
!! ** 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
!!
!! - Compute available flux of heat for surface ablation
!! - Compute snow and sea ice enthalpies
!! - Surface ablation and sublimation
!! - Bottom accretion/ablation
!! - Snow ice formation
!!
!! ** Note : h=max(0,h+dh) are often used to ensure positivity of h.
!! very small negative values can occur otherwise (e.g. -1.e-20)
!!
!! 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 :: ji, jk ! dummy loop indices
INTEGER :: iter ! local integer
REAL(wp) :: ztmelts ! local scalar
REAL(wp) :: zdum
REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment
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) :: z1_rho ! 1/(rhos+rho0-rhoi)
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), DIMENSION(jpij) :: zq_top ! heat for surface ablation (J.m-2)
REAL(wp), DIMENSION(jpij) :: zq_bot ! heat for bottom ablation (J.m-2)
REAL(wp), DIMENSION(jpij) :: zq_rema ! remaining heat at the end of the routine (J.m-2)
REAL(wp), DIMENSION(jpij) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2)
REAL(wp) :: zevap_rema ! remaining mass flux from sublimation (kg.m-2)
REAL(wp) :: zdeltah
REAL(wp), DIMENSION(jpij) :: zsnw ! distribution of snow after wind blowing
INTEGER , DIMENSION(nlay_i) :: icount ! number of layers vanishing by melting
REAL(wp), DIMENSION(0:nlay_i+1) :: zh_i ! ice layer thickness (m)
REAL(wp), DIMENSION(0:nlay_s ) :: zh_s ! snw layer thickness (m)
REAL(wp), DIMENSION(0:nlay_s ) :: ze_s ! snw layer enthalpy (J.m-3)
REAL(wp), DIMENSION(0:nlay_i+1) :: zh_i_old ! old thickness
REAL(wp), DIMENSION(0:nlay_i+1) :: ze_i_old ! old enthalpy
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
REAL(wp) :: zswitch_sal
INTEGER :: num_iter_max ! Heat conservation
!!------------------------------------------------------------------
! Discriminate between time varying salinity and constant
SELECT CASE( nn_icesal ) ! varying salinity or not
CASE( 1 , 3 ) ; zswitch_sal = 0._wp ! prescribed salinity profile
CASE( 2 ) ; zswitch_sal = 1._wp ! varying salinity profile
END SELECT
!
! ! ============================================== !
! ! Available heat for surface and bottom ablation !
! ! ============================================== !
IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN
DO ji = 1, npti
zq_top(ji) = MAX( 0._wp, qml_ice_1d(ji) * rDt_ice )
END DO
ELSE
DO ji = 1, npti
zdum = qns_ice_1d(ji) + qsr_ice_1d(ji) - qtr_ice_top_1d(ji) - qcn_ice_top_1d(ji)
qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
zq_top(ji) = MAX( 0._wp, qml_ice_1d(ji) * rDt_ice )
END DO
ENDIF
!
DO ji = 1, npti
zf_tt(ji) = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) + qtr_ice_bot_1d(ji) * frq_m_1d(ji)
zq_bot(ji) = MAX( 0._wp, zf_tt(ji) * rDt_ice )
END DO
! ! ========== !
! ! Other init !
! ! ========== !
!
! snow distribution over ice after wind blowing
CALL ice_var_snwblow( 1._wp - at_i_1d(1:npti), zsnw(1:npti) )
!
! for snw-ice formation
IF( nn_icesal == 2 ) THEN ; num_iter_max = 5 ! salinity varying in time
ELSE ; num_iter_max = 1
ENDIF
! ! ==================== !
! ! Start main loop here !
! ! ==================== !
! initialize ice layer thicknesses and enthalpies
ze_i_old(0:nlay_i+1) = 0._wp
zh_i_old(0:nlay_i+1) = 0._wp
DO jk = 1, nlay_i
ze_i_old(jk) = h_i_1d(ji) * r1_nlay_i * e_i_1d(ji,jk)
zh_i_old(jk) = h_i_1d(ji) * r1_nlay_i
END DO
!
! initialize snw layer thicknesses and enthalpies
zh_s(0) = 0._wp
ze_s(0) = 0._wp
DO jk = 1, nlay_s
zh_s(jk) = h_s_1d(ji) * r1_nlay_s
ze_s(jk) = e_s_1d(ji,jk)
END DO
!
! ! ============ !
! ! Snow !
! ! ============ !
!
! Internal melting
! ----------------
! IF snow temperature is above freezing point, THEN snow melts (should not happen but sometimes it does)
DO jk = 1, nlay_s
hfx_res_1d (ji) = hfx_res_1d (ji) - ze_s(jk) * zh_s(jk) * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0
wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos * zh_s(jk) * a_i_1d(ji) * r1_Dt_ice ! mass flux
dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(jk)
h_s_1d (ji) = MAX( 0._wp, h_s_1d (ji) - zh_s(jk) )
zh_s(0) = zsnw(ji) * sprecip_1d(ji) * rDt_ice * r1_rhos / at_i_1d(ji) ! thickness of precip
ze_s(0) = MAX( 0._wp, - qprec_ice_1d(ji) ) ! enthalpy of the precip (>0, J.m-3)
hfx_spr_1d(ji) = hfx_spr_1d(ji) + ze_s(0) * zh_s(0) * a_i_1d(ji) * r1_Dt_ice ! heat flux from snow precip (>0, W.m-2)
wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos * zh_s(0) * a_i_1d(ji) * r1_Dt_ice ! mass flux, <0
! Snow melting
! ------------
! If heat still available (zq_top > 0)
! then all snw precip has been melted and we need to melt more snow
DO jk = 0, nlay_s
IF( zh_s(jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN
rswitch = MAX( 0._wp , SIGN( 1._wp , ze_s(jk) - epsi20 ) )
zdum = - rswitch * zq_top(ji) / MAX( ze_s(jk), epsi20 ) ! thickness change
zdum = MAX( zdum , - zh_s(jk) ) ! bound melting
hfx_snw_1d (ji) = hfx_snw_1d (ji) - ze_s(jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! heat used to melt snow(W.m-2, >0)
wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * zdum * a_i_1d(ji) * r1_Dt_ice ! snow melting only = water into the ocean
dh_s_mlt(ji) = dh_s_mlt(ji) + zdum
zq_top (ji) = MAX( 0._wp , zq_top (ji) + zdum * ze_s(jk) )
h_s_1d (ji) = MAX( 0._wp , h_s_1d (ji) + zdum )
zh_s (jk) = MAX( 0._wp , zh_s (jk) + zdum )
!!$ IF( zh_s(jk) == 0._wp ) ze_s(jk) = 0._wp
! 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 = MAX( - evap_ice_1d(ji) * r1_rhos * rDt_ice, - h_s_1d(ji) ) ! amount of snw that sublimates, < 0
zevap_rema = evap_ice_1d(ji) * rDt_ice + zdeltah * rhos ! remaining evap in kg.m-2 (used for ice sublimation later on)
zdum = MAX( -zh_s(jk), zdeltah ) ! snow layer thickness that sublimates, < 0
hfx_sub_1d (ji) = hfx_sub_1d (ji) + ze_s(jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! Heat flux of snw that sublimates [W.m-2], < 0
wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux by sublimation
zh_s (jk) = MAX( 0._wp , zh_s (jk) + zdum )
!!$ IF( zh_s(jk) == 0._wp ) ze_s(jk) = 0._wp
!
! ! ============ !
! ! Ice !
! ! ============ !
! Surface ice melting
!--------------------
DO jk = 1, nlay_i
ztmelts = - rTmlt * sz_i_1d(ji,jk) ! Melting point of layer k [C]
IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting
zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of layer k [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)
zdum = MIN( 0._wp , - zh_i(jk) ) ! internal melting occurs when the internal temperature is above freezing
! this should normally not happen, but sometimes, heat diffusion leads to this
zfmdt = - zdum * rhoi ! Recompute mass flux [kg/m2, >0]
!
dh_i_itm(ji) = dh_i_itm(ji) + zdum ! Cumulate internal melting
!
hfx_res_1d(ji) = hfx_res_1d(ji) + zEi * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0
! ice enthalpy zEi is "sent" to the ocean
wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux
sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux
! using s_i_1d and not sz_i_1d(jk) is ok
ELSE !-- Surface melting
zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of layer k [J/kg, <0]
zEw = rcp * ztmelts ! Specific enthalpy of resulting meltwater [J/kg, <0]
zdE = zEi - zEw ! Specific enthalpy difference < 0
zfmdt = - zq_top(ji) / zdE ! Mass flux to the ocean [kg/m2, >0]
zdum = - zfmdt * r1_rhoi ! Melt of layer jk [m, <0]
zdum = MIN( 0._wp , MAX( zdum , - zh_i(jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0]
zq_top(ji) = MAX( 0._wp , zq_top(ji) - zdum * rhoi * zdE ) ! update available heat
dh_i_sum(ji) = dh_i_sum(ji) + zdum ! Cumulate surface melt
zfmdt = - rhoi * zdum ! Recompute mass flux [kg/m2, >0]
zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0]
hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux [W.m-2], < 0
hfx_sum_1d(ji) = hfx_sum_1d(ji) - zdE * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux used in this process [W.m-2], > 0
wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux
sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux >0
! using s_i_1d and not sz_i_1d(jk) is ok)
END IF
! update thickness
zh_i (jk) = MAX( 0._wp, zh_i (jk) + zdum )
h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum )
!
! update heat content (J.m-2) and layer thickness
ze_i_old(jk) = ze_i_old(jk) + zdum * e_i_1d(ji,jk)
zh_i_old(jk) = zh_i_old(jk) + zdum
zdum = MAX( - zh_i(jk) , - zevap_rema * r1_rhoi )
!
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 = zevap_rema + zdum * rhoi
zh_i (jk) = MAX( 0._wp, zh_i (jk) + zdum )
h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum )
dh_i_sub(ji) = dh_i_sub(ji) + zdum
ze_i_old(jk) = ze_i_old(jk) + zdum * e_i_1d(ji,jk)
zh_i_old(jk) = zh_i_old(jk) + zdum
! 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(jk) ) )
! remaining "potential" evap is sent to ocean
wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema * a_i_1d(ji) * r1_Dt_ice ! <=0 (net evap for the ocean in kg.m-2.s-1)
! Ice Basal growth
!------------------
! Basal growth is driven by heat imbalance at the ice-ocean interface,
! between the inner conductive flux (qcn_ice_bot), from the open water heat flux
! (fhld) and the sensible ice-ocean flux (qsb_ice_bot).
! qcn_ice_bot is positive downwards. qsb_ice_bot 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_bog) depends upon new ice specific enthalpy (zEi),
! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bog)
! -> need for an iterative procedure, which converges quickly
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
IF( zf_tt(ji) < 0._wp ) THEN
DO iter = 1, num_iter_max ! iterations
! 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_bog(ji) * r1_Dt_ice , 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 )
s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji) ! New ice salinity
ztmelts = - rTmlt * s_i_new(ji) ! New ice melting point (C)
zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
zEi = rcpi * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0)
& - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp * ztmelts
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_bog(ji) = rDt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) )
END DO
! Contribution to Energy and Salt Fluxes
zfmdt = - rhoi * dh_i_bog(ji) ! Mass flux x time step (kg/m2, < 0)
hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux to the ocean [W.m-2], >0
hfx_bog_1d(ji) = hfx_bog_1d(ji) - zdE * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux used in this process [W.m-2], <0
wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * dh_i_bog(ji) * a_i_1d(ji) * r1_Dt_ice ! Mass flux, <0
sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * dh_i_bog(ji) * s_i_new(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux, <0
! update thickness
ze_i_old(nlay_i+1) = ze_i_old(nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi)
zh_i_old(nlay_i+1) = zh_i_old(nlay_i+1) + dh_i_bog(ji)
! Ice Basal melt
!---------------
DO jk = nlay_i, 1, -1
IF( zf_tt(ji) > 0._wp .AND. jk > icount(jk) ) THEN ! do not calculate where layer has already disappeared by surface melting
ztmelts = - rTmlt * sz_i_1d(ji,jk) ! Melting point of layer jk (C)
IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting
zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (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)
zdum = MIN( 0._wp , - zh_i(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_itm (ji) = dh_i_itm(ji) + zdum
!
zfmdt = - zdum * rhoi ! Mass flux x time step > 0
!
hfx_res_1d(ji) = hfx_res_1d(ji) + zEi * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0
! ice enthalpy zEi is "sent" to the ocean
wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux
sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux
! using s_i_1d and not sz_i_1d(jk) is ok
ELSE !-- Basal melting
zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0)
zEw = rcp * ztmelts ! Specific enthalpy of meltwater (J/kg, <0)
zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
zfmdt = - zq_bot(ji) / zdE ! Mass flux x time step (kg/m2, >0)
zdum = - zfmdt * r1_rhoi ! Gross thickness change
zdum = MIN( 0._wp , MAX( zdum, - zh_i(jk) ) ) ! bound thickness change
zq_bot(ji) = MAX( 0._wp , zq_bot(ji) - zdum * rhoi * zdE ) ! update available heat. MAX is necessary for roundup errors
dh_i_bom(ji) = dh_i_bom(ji) + zdum ! Update basal melt
zfmdt = - zdum * rhoi ! Mass flux x time step > 0
zQm = zfmdt * zEw ! Heat exchanged with ocean
hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0
hfx_bom_1d(ji) = hfx_bom_1d(ji) - zdE * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat used in this process [W.m-2], >0
wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux
sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux
! using s_i_1d and not sz_i_1d(jk) is ok
ENDIF
! update thickness
zh_i (jk) = MAX( 0._wp, zh_i (jk) + zdum )
h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum )
!
! update heat content (J.m-2) and layer thickness
ze_i_old(jk) = ze_i_old(jk) + zdum * e_i_1d(ji,jk)
zh_i_old(jk) = zh_i_old(jk) + zdum
! Remove snow if ice has melted entirely
! --------------------------------------
IF( h_i_1d(ji) == 0._wp ) THEN
DO jk = 0, nlay_s
hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(jk) * zh_s(jk) * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0
wfx_res_1d(ji) = wfx_res_1d(ji) + rhos * zh_s(jk) * a_i_1d(ji) * r1_Dt_ice ! mass flux
ze_s (jk) = 0._wp
zh_s (jk) = 0._wp
END DO
ENDIF
! Snow load on ice
! -----------------
! When snow load exceeds Archimede's limit and sst is positive,
! snow-ice formation (next bloc) can lead to negative ice enthalpy.
! Therefore we consider here that this excess of snow falls into the ocean
zdeltah = h_s_1d(ji) + h_i_1d(ji) * (rhoi-rho0) * r1_rhos
DO jk = 0, nlay_s
IF( zdeltah > 0._wp .AND. sst_1d(ji) > 0._wp ) THEN
hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0
wfx_res_1d(ji) = wfx_res_1d(ji) + rhos * zdum * a_i_1d(ji) * r1_Dt_ice ! mass flux
h_s_1d(ji) = MAX( 0._wp, h_s_1d(ji) - zdum )
zh_s (jk) = MAX( 0._wp, zh_s (jk) - zdum )
! Snow-Ice formation
! ------------------
! When snow load exceeds Archimede's limit, snow-ice interface goes down under sea-level,
! flooding of seawater transforms snow into ice. Thickness that is transformed is dh_snowice (positive for the ice)
!
dh_snowice(ji) = MAX( 0._wp , ( rhos * h_s_1d(ji) + (rhoi-rho0) * h_i_1d(ji) ) * z1_rho )
h_i_1d(ji) = h_i_1d(ji) + dh_snowice(ji)
h_s_1d(ji) = h_s_1d(ji) - dh_snowice(ji)
! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
zfmdt = ( rhos - rhoi ) * dh_snowice(ji) ! <0
zEw = rcp * sst_1d(ji)
zQm = zfmdt * zEw
hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux
sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Salt flux
! Case constant salinity in time: virtual salt flux to keep salinity constant
IF( nn_icesal /= 2 ) THEN
sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d(ji) * zfmdt * a_i_1d(ji) * r1_Dt_ice & ! put back sss_m into the ocean
& - s_i_1d(ji) * dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_Dt_ice ! and get rn_icesal from the ocean
ENDIF
! 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) - dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_Dt_ice
wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + dh_snowice(ji) * rhos * a_i_1d(ji) * r1_Dt_ice
! update thickness
zh_i(0) = zh_i(0) + dh_snowice(ji)
zdeltah = dh_snowice(ji)
zh_i_old(0) = zh_i_old(0) + dh_snowice(ji)
ze_i_old(0) = ze_i_old(0) + zfmdt * zEw ! 1st part (sea water enthalpy)
!
DO jk = nlay_s, 0, -1 ! flooding of snow starts from the base
zdum = MIN( zdeltah, zh_s(jk) ) ! amount of snw that floods, > 0
zh_s(jk) = MAX( 0._wp, zh_s(jk) - zdum ) ! remove some snow thickness
ze_i_old(0) = ze_i_old(0) + zdum * ze_s(jk) ! 2nd part (snow enthalpy)
!!$ ! --- Update snow diags --- !
!!$ !!clem: this is wrong. dh_s_tot is not used anyway
!!$ DO ji = 1, npti
!!$ dh_s_tot(ji) = dh_s_tot(ji) + dh_s_mlt(ji) + zdeltah + zdh_s_sub(ji) - dh_snowice(ji)
! Remapping of snw enthalpy on a regular grid
!--------------------------------------------
e_s_1d(ji,:) = snw_ent( zh_s(:), ze_s(:) )
! recalculate t_s_1d from e_s_1d
IF( h_s_1d(ji) > 0._wp ) THEN
DO jk = 1, nlay_s
t_s_1d(ji,jk) = rt0 + ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi )
END DO
ENDIF
! Remapping of ice enthalpy on a regular grid
!--------------------------------------------
! ! ================== !
! ! End main loop here !
! ! ================== !
! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 ---
WHERE( h_i_1d(1:npti) == 0._wp )
a_i_1d (1:npti) = 0._wp
h_s_1d (1:npti) = 0._wp
t_su_1d(1:npti) = rt0
END WHERE
END SUBROUTINE ice_thd_dh
!!-------------------------------------------------------------------
!! *** ROUTINE snw_ent ***
!!
!! ** Purpose :
!! This routine computes new vertical grids in the snow,
!! and consistently redistributes temperatures.
!! Redistribution is made so as to ensure to energy conservation
!!
!!
!! ** Method : linear conservative remapping
!!
!! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses
!! 2) linear remapping on the new layers
!!
!! ------------ cum0(0) ------------- cum1(0)
!! NEW -------------
!! ------------ cum0(1) ==> -------------
!! ... -------------
!! ------------ -------------
!! ------------ cum0(nlay_s+1) ------------- cum1(nlay_s)
!!
!!
!! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005
!!-------------------------------------------------------------------
REAL(wp), DIMENSION(0:nlay_s), INTENT(in) :: ph_old ! old thicknesses (m)
REAL(wp), DIMENSION(0:nlay_s), INTENT(in) :: pe_old ! old enthlapies (J.m-3)
REAL(wp), DIMENSION(1:nlay_s) :: snw_ent ! new enthlapies (J.m-3, remapped)
!
INTEGER :: ji ! dummy loop indices
INTEGER :: jk0, jk1 ! old/new layer indices
REAL(wp), DIMENSION(0:nlay_s+1) :: zeh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces
REAL(wp), DIMENSION(0:nlay_s) :: zeh_cum1, zh_cum1 ! new cumulative enthlapies and layers interfaces
REAL(wp) :: zhnew ! new layers thicknesses
!!-------------------------------------------------------------------
!--------------------------------------------------------------------------
! 1) Cumulative integral of old enthalpy * thickness and layers interfaces
!--------------------------------------------------------------------------
zeh_cum0(0) = 0._wp
zh_cum0 (0) = 0._wp
DO jk0 = 1, nlay_s+1
zeh_cum0(jk0) = zeh_cum0(jk0-1) + pe_old(jk0-1) * ph_old(jk0-1)
zh_cum0 (jk0) = zh_cum0 (jk0-1) + ph_old(jk0-1)
END DO
!------------------------------------
! 2) Interpolation on the new layers
!------------------------------------
! new layer thickesses
zhnew = SUM( ph_old(0:nlay_s) ) * r1_nlay_s
! new layers interfaces
zh_cum1(0) = 0._wp
DO jk1 = 1, nlay_s
zh_cum1(jk1) = zh_cum1(jk1-1) + zhnew
END DO
zeh_cum1(0:nlay_s) = 0._wp
! new cumulative q*h => linear interpolation
DO jk0 = 1, nlay_s+1
DO jk1 = 1, nlay_s-1
IF( zh_cum1(jk1) <= zh_cum0(jk0) .AND. zh_cum1(jk1) > zh_cum0(jk0-1) ) THEN
zeh_cum1(jk1) = ( zeh_cum0(jk0-1) * ( zh_cum0(jk0) - zh_cum1(jk1 ) ) + &
& zeh_cum0(jk0 ) * ( zh_cum1(jk1) - zh_cum0(jk0-1) ) ) &
& / ( zh_cum0(jk0) - zh_cum0(jk0-1) )
ENDIF
END DO
! to ensure that total heat content is strictly conserved, set:
zeh_cum1(nlay_s) = zeh_cum0(nlay_s+1)
! new enthalpies
DO jk1 = 1, nlay_s
zswitch = MAX( 0._wp , SIGN( 1._wp , zhnew - epsi20 ) )
snw_ent(jk1) = zswitch * ( zeh_cum1(jk1) - zeh_cum1(jk1-1) ) / MAX( zhnew, epsi20 )
END DO
!!-------------------------------------------------------------------
!!
!! ** Purpose :
!! This routine computes new vertical grids in the ice,
!! and consistently redistributes temperatures.
!! Redistribution is made so as to ensure to energy conservation
!!
!!
!! ** Method : linear conservative remapping
!!
!! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses
!! 2) linear remapping on the new layers
!!
!! ------------ cum0(0) ------------- cum1(0)
!! NEW -------------
!! ------------ cum0(1) ==> -------------
!! ... -------------
!! ------------ -------------
!! ------------ cum0(nlay_i+2) ------------- cum1(nlay_i)
!!
!!
!! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005
!!-------------------------------------------------------------------
REAL(wp), DIMENSION(0:nlay_i+1), INTENT(in) :: ph_old, pe_old ! old tickness and enthlapy
REAL(wp), DIMENSION(1:nlay_i) :: ice_ent1 ! new enthlapies (J.m-3, remapped)
!
INTEGER :: ji ! dummy loop indices
INTEGER :: jk0, jk1 ! old/new layer indices
REAL(wp) :: zswitch
!
REAL(wp), DIMENSION(0:nlay_i+2) :: zeh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces
REAL(wp), DIMENSION(0:nlay_i) :: zeh_cum1, zh_cum1 ! new cumulative enthlapies and layers interfaces
REAL(wp) :: zhnew ! new layers thicknesses
!!-------------------------------------------------------------------
!--------------------------------------------------------------------------
! 1) Cumulative integral of old enthalpy * thickness and layers interfaces
!--------------------------------------------------------------------------
zeh_cum0(0) = 0._wp
zh_cum0 (0) = 0._wp
DO jk0 = 1, nlay_i+2
zeh_cum0(jk0) = zeh_cum0(jk0-1) + pe_old(jk0-1)
zh_cum0 (jk0) = zh_cum0 (jk0-1) + ph_old(jk0-1)
!------------------------------------
! 2) Interpolation on the new layers
!------------------------------------
! new layer thickesses
zhnew = SUM( ph_old(0:nlay_i+1) ) * r1_nlay_i
! new layers interfaces
zh_cum1(0) = 0._wp
DO jk1 = 1, nlay_i
zh_cum1(jk1) = zh_cum1(jk1-1) + zhnew
END DO
zeh_cum1(0:nlay_i) = 0._wp
! new cumulative q*h => linear interpolation
DO jk0 = 1, nlay_i+2
DO jk1 = 1, nlay_i-1
IF( zh_cum1(jk1) <= zh_cum0(jk0) .AND. zh_cum1(jk1) > zh_cum0(jk0-1) ) THEN
zeh_cum1(jk1) = ( zeh_cum0(jk0-1) * ( zh_cum0(jk0) - zh_cum1(jk1 ) ) + &
& zeh_cum0(jk0 ) * ( zh_cum1(jk1) - zh_cum0(jk0-1) ) ) &
& / ( zh_cum0(jk0) - zh_cum0(jk0-1) )
ENDIF
END DO
END DO
! to ensure that total heat content is strictly conserved, set:
zeh_cum1(nlay_i) = zeh_cum0(nlay_i+2)
! new enthalpies
DO jk1 = 1, nlay_i
zswitch = MAX( 0._wp , SIGN( 1._wp , zhnew - epsi20 ) )
ice_ent1(jk1) = zswitch * MAX( 0._wp, zeh_cum1(jk1) - zeh_cum1(jk1-1) ) / MAX( zhnew, epsi20 ) ! max for roundoff error
END DO
! --- diag error on heat remapping --- !
! comment: if input h_old and eh_old are already multiplied by a_i (as in icethd_do),
! then we should not (* a_i) again but not important since this is just to check that remap error is ~0
! hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_Dt_ice * &
! & ( SUM( pe_new(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_old(ji,0:nlay_i+1) ) )
#else
!!----------------------------------------------------------------------
!! Default option NO SI3 sea-ice model
!!----------------------------------------------------------------------
#endif
!!======================================================================
END MODULE icethd_dh