diff --git a/src/ICE/icethd.F90 b/src/ICE/icethd.F90
index 218c628271c1368e199de9034744053ce311d222..448ae5aa9177444c6e73320eeae3a134211db853 100644
--- a/src/ICE/icethd.F90
+++ b/src/ICE/icethd.F90
@@ -71,7 +71,6 @@ CONTAINS
!! - call ice_thd_zdf for vertical heat diffusion
!! - call ice_thd_dh for vertical ice growth and melt
!! - call ice_thd_pnd for melt ponds
- !! - call ice_thd_ent for enthalpy remapping
!! - call ice_thd_sal for ice desalination
!! - call ice_thd_temp to retrieve temperature from ice enthalpy
!! - call ice_thd_mono for extra lateral ice melt if active virtual thickness distribution
diff --git a/src/ICE/icethd_dh.F90 b/src/ICE/icethd_dh.F90
index 6e51cdc54ac1e203553fa575317677827f0c1c1a..dc3c420181805bfd8e6d5c09eb8c4fe3b832ef4b 100644
--- a/src/ICE/icethd_dh.F90
+++ b/src/ICE/icethd_dh.F90
@@ -18,7 +18,6 @@ MODULE icethd_dh
USE ice ! sea-ice: variables
USE ice1D ! sea-ice: thermodynamics variables
USE icethd_sal ! sea-ice: salinity profiles
- USE icethd_ent ! sea-ice: enthalpy redistribution
USE icevar ! for CALL ice_var_snwblow
!
USE in_out_manager ! I/O manager
@@ -91,11 +90,11 @@ CONTAINS
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(jpij,0:nlay_i+1) :: zh_i ! ice layer thickness (m)
- REAL(wp), DIMENSION(jpij,0:nlay_s ) :: zh_s ! snw layer thickness (m)
- REAL(wp), DIMENSION(jpij,0:nlay_s ) :: ze_s ! snw layer enthalpy (J.m-3)
- REAL(wp), DIMENSION(jpij,0:nlay_i+1) :: zh_i_old ! old thickness
- REAL(wp), DIMENSION(jpij,0:nlay_i+1) :: ze_i_old ! old enthalpy
+ 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
REAL(wp) :: zswitch_sal
@@ -108,27 +107,6 @@ CONTAINS
CASE( 2 ) ; zswitch_sal = 1._wp ! varying salinity profile
END SELECT
- ! initialize ice layer thicknesses and enthalpies
- ze_i_old(1:npti,0:nlay_i+1) = 0._wp
- zh_i_old(1:npti,0:nlay_i+1) = 0._wp
- zh_i (1:npti,0:nlay_i+1) = 0._wp
- DO jk = 1, nlay_i
- DO ji = 1, npti
- ze_i_old(ji,jk) = h_i_1d(ji) * r1_nlay_i * e_i_1d(ji,jk)
- zh_i_old(ji,jk) = h_i_1d(ji) * r1_nlay_i
- zh_i (ji,jk) = h_i_1d(ji) * r1_nlay_i
- END DO
- END DO
- !
- ! initialize snw layer thicknesses and enthalpies
- zh_s(1:npti,0) = 0._wp
- ze_s(1:npti,0) = 0._wp
- DO jk = 1, nlay_s
- DO ji = 1, npti
- zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
- ze_s(ji,jk) = e_s_1d(ji,jk)
- END DO
- END DO
!
! ! ============================================== !
! ! Available heat for surface and bottom ablation !
@@ -166,6 +144,24 @@ CONTAINS
DO ji = 1, npti
+ ! initialize ice layer thicknesses and enthalpies
+ ze_i_old(0:nlay_i+1) = 0._wp
+ zh_i_old(0:nlay_i+1) = 0._wp
+ zh_i ( 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
+ zh_i ( 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 !
! ! ============ !
@@ -175,13 +171,13 @@ CONTAINS
! IF snow temperature is above freezing point, THEN snow melts (should not happen but sometimes it does)
DO jk = 1, nlay_s
IF( t_s_1d(ji,jk) > rt0 ) THEN
- hfx_res_1d (ji) = hfx_res_1d (ji) - ze_s(ji,jk) * zh_s(ji,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(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! mass flux
+ 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
! updates
- dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk)
- h_s_1d (ji) = MAX( 0._wp, h_s_1d (ji) - zh_s(ji,jk) )
- zh_s (ji,jk) = 0._wp
- ze_s (ji,jk) = 0._wp
+ 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 (jk) = 0._wp
+ ze_s (jk) = 0._wp
END IF
END DO
@@ -189,14 +185,14 @@ CONTAINS
!-------------------
IF( sprecip_1d(ji) > 0._wp ) THEN
- zh_s(ji,0) = zsnw(ji) * sprecip_1d(ji) * rDt_ice * r1_rhos / at_i_1d(ji) ! thickness of precip
- ze_s(ji,0) = MAX( 0._wp, - qprec_ice_1d(ji) ) ! enthalpy of the precip (>0, J.m-3)
+ 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(ji,0) * zh_s(ji,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(ji,0) * a_i_1d(ji) * r1_Dt_ice ! mass flux, <0
+ 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
!
! update thickness
- h_s_1d(ji) = h_s_1d(ji) + zh_s(ji,0)
+ h_s_1d(ji) = h_s_1d(ji) + zh_s(0)
ENDIF
! Snow melting
@@ -204,21 +200,21 @@ CONTAINS
! 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(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN
+ IF( zh_s(jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN
!
- rswitch = MAX( 0._wp , SIGN( 1._wp , ze_s(ji,jk) - epsi20 ) )
- zdum = - rswitch * zq_top(ji) / MAX( ze_s(ji,jk), epsi20 ) ! thickness change
- zdum = MAX( zdum , - zh_s(ji,jk) ) ! bound melting
+ 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(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! heat used to melt snow(W.m-2, >0)
+ 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
! updates available heat + thickness
dh_s_mlt(ji) = dh_s_mlt(ji) + zdum
- zq_top (ji) = MAX( 0._wp , zq_top (ji) + zdum * ze_s(ji,jk) )
+ 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 (ji,jk) = MAX( 0._wp , zh_s (ji,jk) + zdum )
-!!$ IF( zh_s(ji,jk) == 0._wp ) ze_s(ji,jk) = 0._wp
+ zh_s (jk) = MAX( 0._wp , zh_s (jk) + zdum )
+!!$ IF( zh_s(jk) == 0._wp ) ze_s(jk) = 0._wp
!
ENDIF
END DO
@@ -231,15 +227,15 @@ CONTAINS
zevap_rema = evap_ice_1d(ji) * rDt_ice + zdeltah * rhos ! remaining evap in kg.m-2 (used for ice sublimation later on)
DO jk = 0, nlay_s
- zdum = MAX( -zh_s(ji,jk), zdeltah ) ! snow layer thickness that sublimates, < 0
+ zdum = MAX( -zh_s(jk), zdeltah ) ! snow layer thickness that sublimates, < 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
+ 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
! update thickness
h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdum )
- zh_s (ji,jk) = MAX( 0._wp , zh_s (ji,jk) + zdum )
-!!$ IF( zh_s(ji,jk) == 0._wp ) ze_s(ji,jk) = 0._wp
+ zh_s (jk) = MAX( 0._wp , zh_s (jk) + zdum )
+!!$ IF( zh_s(jk) == 0._wp ) ze_s(jk) = 0._wp
! update sublimation left
zdeltah = MIN( zdeltah - zdum, 0._wp )
@@ -260,7 +256,7 @@ CONTAINS
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(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing
+ 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]
!
@@ -281,7 +277,7 @@ CONTAINS
zdum = - zfmdt * r1_rhoi ! Melt of layer jk [m, <0]
- zdum = MIN( 0._wp , MAX( zdum , - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [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
@@ -298,17 +294,17 @@ CONTAINS
! using s_i_1d and not sz_i_1d(jk) is ok)
END IF
! update thickness
- zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum )
+ 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(ji,jk) = ze_i_old(ji,jk) + zdum * e_i_1d(ji,jk)
- zh_i_old(ji,jk) = zh_i_old(ji,jk) + zdum
+ ze_i_old(jk) = ze_i_old(jk) + zdum * e_i_1d(ji,jk)
+ zh_i_old(jk) = zh_i_old(jk) + zdum
!
!
! Ice sublimation
! ---------------
- zdum = MAX( - zh_i(ji,jk) , - zevap_rema * r1_rhoi )
+ 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
@@ -318,18 +314,18 @@ CONTAINS
! if all ice is melted. => must be corrected
! update remaining mass flux and thickness
zevap_rema = zevap_rema + zdum * rhoi
- zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum )
+ 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
! update heat content (J.m-2) and layer thickness
- ze_i_old(ji,jk) = ze_i_old(ji,jk) + zdum * e_i_1d(ji,jk)
- zh_i_old(ji,jk) = zh_i_old(ji,jk) + 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(ji,jk) ) )
+ rswitch = MAX( 0._wp , SIGN( 1._wp , - zh_i(jk) ) )
icount(jk) = NINT( rswitch )
END DO
@@ -391,12 +387,12 @@ CONTAINS
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
- zh_i(ji,nlay_i+1) = zh_i(ji,nlay_i+1) + dh_i_bog(ji)
+ zh_i(nlay_i+1) = zh_i(nlay_i+1) + dh_i_bog(ji)
h_i_1d(ji) = h_i_1d(ji) + dh_i_bog(ji)
! update heat content (J.m-2) and layer thickness
- ze_i_old(ji,nlay_i+1) = ze_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi)
- zh_i_old(ji,nlay_i+1) = zh_i_old(ji,nlay_i+1) + dh_i_bog(ji)
+ 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)
ENDIF
@@ -413,7 +409,7 @@ CONTAINS
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(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing
+ 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
!
@@ -434,7 +430,7 @@ CONTAINS
zdum = - zfmdt * r1_rhoi ! Gross thickness change
- zdum = MIN( 0._wp , MAX( zdum, - zh_i(ji,jk) ) ) ! bound 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
@@ -451,29 +447,29 @@ CONTAINS
! using s_i_1d and not sz_i_1d(jk) is ok
ENDIF
! update thickness
- zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum )
+ 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(ji,jk) = ze_i_old(ji,jk) + zdum * e_i_1d(ji,jk)
- zh_i_old(ji,jk) = zh_i_old(ji,jk) + zdum
+ ze_i_old(jk) = ze_i_old(jk) + zdum * e_i_1d(ji,jk)
+ zh_i_old(jk) = zh_i_old(jk) + zdum
ENDIF
END DO
! Remove snow if ice has melted entirely
! --------------------------------------
- DO jk = 0, nlay_s
- IF( h_i_1d(ji) == 0._wp ) THEN
+ IF( h_i_1d(ji) == 0._wp ) THEN
+ DO jk = 0, nlay_s
! mass & energy loss to the ocean
- hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zh_s(ji,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(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! mass flux
+ 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
! update thickness and energy
h_s_1d(ji) = 0._wp
- ze_s (ji,jk) = 0._wp
- zh_s (ji,jk) = 0._wp
- ENDIF
- END DO
+ ze_s (jk) = 0._wp
+ zh_s (jk) = 0._wp
+ END DO
+ ENDIF
! Snow load on ice
! -----------------
@@ -484,13 +480,13 @@ CONTAINS
DO jk = 0, nlay_s
IF( zdeltah > 0._wp .AND. sst_1d(ji) > 0._wp ) THEN
! snow layer thickness that falls into the ocean
- zdum = MIN( zdeltah , zh_s(ji,jk) )
+ zdum = MIN( zdeltah , zh_s(jk) )
! mass & energy loss to the ocean
- hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0
+ 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
! update thickness and energy
h_s_1d(ji) = MAX( 0._wp, h_s_1d(ji) - zdum )
- zh_s (ji,jk) = MAX( 0._wp, zh_s(ji,jk) - zdum )
+ zh_s (jk) = MAX( 0._wp, zh_s(jk) - zdum )
! update snow thickness that still has to fall
zdeltah = MAX( 0._wp, zdeltah - zdum )
ENDIF
@@ -525,18 +521,18 @@ CONTAINS
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(ji,0) = zh_i(ji,0) + dh_snowice(ji)
+ zh_i(0) = zh_i(0) + dh_snowice(ji)
zdeltah = dh_snowice(ji)
! update heat content (J.m-2) and layer thickness
- zh_i_old(ji,0) = zh_i_old(ji,0) + dh_snowice(ji)
- ze_i_old(ji,0) = ze_i_old(ji,0) + zfmdt * zEw ! 1st part (sea water enthalpy)
+ 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(ji,jk) ) ! amount of snw that floods, > 0
- zh_s(ji,jk) = MAX( 0._wp, zh_s(ji,jk) - zdum ) ! remove some snow thickness
- ze_i_old(ji,0) = ze_i_old(ji,0) + zdum * ze_s(ji,jk) ! 2nd part (snow enthalpy)
+ 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 dh_snowice
zdeltah = MAX( 0._wp, zdeltah - zdum )
END DO
@@ -548,26 +544,28 @@ CONTAINS
!!$ dh_s_tot(ji) = dh_s_tot(ji) + dh_s_mlt(ji) + zdeltah + zdh_s_sub(ji) - dh_snowice(ji)
!!$ END DO
!
- END DO ! npti
- !
- ! Remapping of snw enthalpy on a regular grid
- !--------------------------------------------
- CALL snw_ent( zh_s, ze_s, e_s_1d )
-
- ! recalculate t_s_1d from e_s_1d
- DO jk = 1, nlay_s
- DO ji = 1,npti
- IF( h_s_1d(ji) > 0._wp ) THEN
+ ! 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 )
- ELSE
+ END DO
+ ELSE
+ DO jk = 1, nlay_s
t_s_1d(ji,jk) = rt0
- ENDIF
- END DO
- END DO
+ END DO
+ ENDIF
+
+ ! Remapping of ice enthalpy on a regular grid
+ !--------------------------------------------
+ e_i_1d(ji,:) = ice_ent( zh_i_old(:), ze_i_old(:) )
- ! Remapping of ice enthalpy on a regular grid
- !--------------------------------------------
- CALL ice_thd_ent( zh_i_old, ze_i_old, e_i_1d(1:npti,:) )
+ END DO ! npti
+
+
! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 ---
WHERE( h_i_1d(1:npti) == 0._wp )
@@ -578,7 +576,7 @@ CONTAINS
END SUBROUTINE ice_thd_dh
- SUBROUTINE snw_ent( ph_old, pe_old, pe_new )
+ PURE FUNCTION snw_ent( ph_old, pe_old )
!!-------------------------------------------------------------------
!! *** ROUTINE snw_ent ***
!!
@@ -603,66 +601,152 @@ CONTAINS
!!
!! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005
!!-------------------------------------------------------------------
- REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in ) :: ph_old ! old thicknesses (m)
- REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in ) :: pe_old ! old enthlapies (J.m-3)
- REAL(wp), DIMENSION(jpij,1:nlay_s), INTENT(inout) :: pe_new ! new enthlapies (J.m-3, remapped)
+ 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) :: zswitch
!
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
!!-------------------------------------------------------------------
- DO ji = 1, npti
- !--------------------------------------------------------------------------
- ! 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(ji,jk0-1) * ph_old(ji,jk0-1)
- zh_cum0 (jk0) = zh_cum0 (jk0-1) + ph_old(ji,jk0-1)
- END DO
+ !--------------------------------------------------------------------------
+ ! 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(ji,0:nlay_s) ) * r1_nlay_s
+ !------------------------------------
+ ! 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
+ ! 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
+ 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)
+ 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
- ! new enthalpies
- DO jk1 = 1, nlay_s
- rswitch = MAX( 0._wp , SIGN( 1._wp , zhnew - epsi20 ) )
- pe_new(ji,jk1) = rswitch * ( zeh_cum1(jk1) - zeh_cum1(jk1-1) ) / MAX( zhnew, epsi20 )
- END DO
+ END FUNCTION snw_ent
+
+ PURE FUNCTION ice_ent( ph_old, pe_old )
+ !!-------------------------------------------------------------------
+ !! *** ROUTINE ice_ent ***
+ !!
+ !! ** 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_ent ! 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)
END DO
- END SUBROUTINE snw_ent
+ !------------------------------------
+ ! 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_ent(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) ) )
+
+
+ END FUNCTION ice_ent
+
#else
!!----------------------------------------------------------------------
!! Default option NO SI3 sea-ice model
diff --git a/src/ICE/icethd_do.F90 b/src/ICE/icethd_do.F90
index 358e59ea2dfcbbbe28d2a8a51c23fcbaff75fbe3..6867b76a3aa22550c12f601cad9335117d8aa1ce 100644
--- a/src/ICE/icethd_do.F90
+++ b/src/ICE/icethd_do.F90
@@ -21,7 +21,6 @@ MODULE icethd_do
USE ice ! sea-ice: variables
USE icetab ! sea-ice: 2D <==> 1D
USE icectl ! sea-ice: conservation
- USE icethd_ent ! sea-ice: thermodynamics, enthalpy
USE icevar ! sea-ice: operations
USE icethd_sal ! sea-ice: salinity profiles
!
@@ -97,7 +96,7 @@ CONTAINS
REAL(wp), DIMENSION(jpij) :: zh_newice ! thickness of accreted ice
REAL(wp), DIMENSION(jpij) :: zfraz_frac_1d ! relative ice / frazil velocity (1D vector)
!
- REAL(wp), DIMENSION(jpij,0:nlay_i+1,jpl) :: zh_i_old, ze_i_old
+ REAL(wp), DIMENSION(0:nlay_i+1) :: zh_i_old, ze_i_old
!!-----------------------------------------------------------------------!
IF( ln_icediachk ) CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft )
@@ -254,11 +253,11 @@ CONTAINS
DO jl = 1, jpl
! for remapping
- zh_i_old(ji,0:nlay_i+1,jl) = 0._wp
- ze_i_old(ji,0:nlay_i+1,jl) = 0._wp
+ zh_i_old(0:nlay_i+1) = 0._wp
+ ze_i_old(0:nlay_i+1) = 0._wp
DO jk = 1, nlay_i
- zh_i_old(ji,jk,jl) = v_i_2d(ji,jl) * r1_nlay_i
- ze_i_old(ji,jk,jl) = e_i_2d(ji,jk,jl) * v_i_2d(ji,jl) * r1_nlay_i
+ zh_i_old(jk) = v_i_2d(ji,jl) * r1_nlay_i
+ ze_i_old(jk) = e_i_2d(ji,jk,jl) * v_i_2d(ji,jl) * r1_nlay_i
END DO
! new volumes including lateral/bottom accretion + residual
@@ -267,21 +266,18 @@ CONTAINS
a_i_2d(ji,jl) = rswitch * a_i_2d(ji,jl)
v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newfra
! for remapping
- zh_i_old(ji,nlay_i+1,jl) = zv_newfra
- ze_i_old(ji,nlay_i+1,jl) = ze_newice * zv_newfra
- END DO
+ zh_i_old(nlay_i+1) = zv_newfra
+ ze_i_old(nlay_i+1) = ze_newice * zv_newfra
- ! --- Update salinity --- !
- DO jl = 1, jpl
+ ! --- Update salinity --- !
sv_i_2d(ji,jl) = sv_i_2d(ji,jl) + zs_newice(ji) * ( v_i_2d(ji,jl) - zv_b(jl) )
+
+ ! --- Ice enthalpy remapping --- !
+ e_i_2d(ji,:,jl) = ice_ent( zh_i_old(:), ze_i_old(:) )
END DO
END DO ! npti
- ! --- Ice enthalpy remapping --- !
- DO jl = 1, jpl
- CALL ice_thd_ent( zh_i_old(:,:,jl), ze_i_old(:,:,jl), e_i_2d(1:npti,:,jl) )
- END DO
! Change units for e_i
DO jl = 1, jpl
@@ -397,6 +393,94 @@ CONTAINS
ENDIF
END SUBROUTINE ice_thd_frazil
+ PURE FUNCTION ice_ent( ph_old, pe_old )
+ !!-------------------------------------------------------------------
+ !! *** ROUTINE ice_ent ***
+ !!
+ !! ** 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_ent ! 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)
+ END DO
+
+ !------------------------------------
+ ! 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_ent(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) ) )
+
+
+ END FUNCTION ice_ent
+
SUBROUTINE ice_thd_do_init
!!-----------------------------------------------------------------------
diff --git a/src/ICE/iceupdate.F90 b/src/ICE/iceupdate.F90
index aac48c07cdb3eab6d8f77af8674db3215083f107..8efb5c2d0b2f2eb2c7c6d5e2f00afaa15264bcd5 100644
--- a/src/ICE/iceupdate.F90
+++ b/src/ICE/iceupdate.F90
@@ -55,7 +55,7 @@ CONTAINS
!!-------------------------------------------------------------------
!! *** ROUTINE ice_update_alloc ***
!!-------------------------------------------------------------------
- ALLOCATE( utau_oce(A2D(0)), vtau_oce(A2D(0)), tmod_io(A2D(1)), STAT=ice_update_alloc )
+ ALLOCATE( utau_oce(A2D(0)), vtau_oce(A2D(0)), tmod_io(jpi,jpj), STAT=ice_update_alloc )
!
CALL mpp_sum( 'iceupdate', ice_update_alloc )
IF( ice_update_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_update_alloc: failed to allocate arrays' )
diff --git a/src/ICE/icevar.F90 b/src/ICE/icevar.F90
index 935198289d2a3ca8728d494c873bf84a335a6938..89f5c1859a33566a59e8f4ebef4b4af49b858953 100644
--- a/src/ICE/icevar.F90
+++ b/src/ICE/icevar.F90
@@ -138,10 +138,8 @@ CONTAINS
st_i(ji,jj) = SUM( sv_i(ji,jj,:) )
et_s(ji,jj) = SUM( SUM( e_s (ji,jj,:,:), dim=2 ) )
et_i(ji,jj) = SUM( SUM( e_i (ji,jj,:,:), dim=2 ) )
- END_2D
- !
- !!GS: tm_su always needed by ABL over sea-ice
- DO_2D( 0, 0, 0, 0 )
+ !
+ !!GS: tm_su always needed by ABL over sea-ice
IF( at_i(ji,jj) <= epsi20 ) THEN
z1_at_i(ji,jj) = 0._wp
tm_su (ji,jj) = rt0
@@ -284,20 +282,24 @@ CONTAINS
!-------------------
zlay_i = REAL( nlay_i , wp ) ! number of layers
DO jl = 1, jpl
- DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_i )
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area
- !
- ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3]
- ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C]
- ! Conversion q(S,T) -> T (second order equation)
- zbbb = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus
- zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) )
- t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts
- !
+ DO jk = 1, nlay_i
+ !
+ ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3]
+ ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C]
+ ! Conversion q(S,T) -> T (second order equation)
+ zbbb = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus
+ zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) )
+ t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts
+ !
+ END DO
ELSE !--- no ice
- t_i(ji,jj,jk,jl) = rt0
+ DO jk = 1, nlay_i
+ t_i(ji,jj,jk,jl) = rt0
+ END DO
ENDIF
- END_3D
+ END_2D
END DO
!--------------------