diff --git a/src/ICE/icethd_dh.F90 b/src/ICE/icethd_dh.F90
index dc3c420181805bfd8e6d5c09eb8c4fe3b832ef4b..a0d40b7c8682ebb755dda8952ce74733f3d90bc8 100644
--- a/src/ICE/icethd_dh.F90
+++ b/src/ICE/icethd_dh.F90
@@ -106,52 +106,53 @@ CONTAINS
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
-
-
- CALL ice_var_snwblow( 1._wp - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing
-
+ 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
z1_rho = 1._wp / ( rhos+rho0-rhoi )
-
+ !
+ ! number of iterations for new sea ice
IF( nn_icesal == 2 ) THEN ; num_iter_max = 5 ! salinity varying in time
ELSE ; num_iter_max = 1
ENDIF
-
+ ! ! ==================== !
+ ! ! Start main loop here !
+ ! ! ==================== !
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
+ 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
+ zh_i (jk) = h_i_1d(ji) * r1_nlay_i
END DO
!
! initialize snw layer thicknesses and enthalpies
@@ -172,10 +173,10 @@ CONTAINS
DO jk = 1, nlay_s
IF( t_s_1d(ji,jk) > rt0 ) THEN
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
+ 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(jk)
- h_s_1d (ji) = MAX( 0._wp, h_s_1d (ji) - zh_s(jk) )
+ 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
@@ -183,13 +184,12 @@ CONTAINS
! Snow precipitation
!-------------------
-
IF( sprecip_1d(ji) > 0._wp ) THEN
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
+ 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(0)
@@ -207,12 +207,12 @@ CONTAINS
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
+ 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(jk) )
- h_s_1d (ji) = MAX( 0._wp , h_s_1d (ji) + zdum )
+ 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
!
@@ -224,16 +224,15 @@ CONTAINS
! 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)
-
+ 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(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
+ 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 )
+ 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
@@ -256,7 +255,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(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]
!
@@ -294,8 +293,8 @@ CONTAINS
! 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 )
+ 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)
@@ -313,10 +312,10 @@ CONTAINS
! 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
+ 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
! update heat content (J.m-2) and layer thickness
ze_i_old(jk) = ze_i_old(jk) + zdum * e_i_1d(ji,jk)
@@ -347,7 +346,6 @@ CONTAINS
! 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
-
IF( zf_tt(ji) < 0._wp ) THEN
DO iter = 1, num_iter_max ! iterations
@@ -388,7 +386,7 @@ CONTAINS
! update thickness
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)
+ h_i_1d(ji) = h_i_1d(ji) + dh_i_bog(ji)
! update heat content (J.m-2) and layer thickness
ze_i_old(nlay_i+1) = ze_i_old(nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi)
@@ -396,7 +394,6 @@ CONTAINS
ENDIF
-
! Ice Basal melt
!---------------
DO jk = nlay_i, 1, -1
@@ -447,8 +444,8 @@ CONTAINS
! 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 )
+ 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)
@@ -462,10 +459,10 @@ CONTAINS
DO jk = 0, nlay_s
! mass & energy loss to the ocean
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
+ 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
+ h_s_1d(ji) = 0._wp
ze_s (jk) = 0._wp
zh_s (jk) = 0._wp
END DO
@@ -483,10 +480,10 @@ CONTAINS
zdum = MIN( zdeltah , zh_s(jk) )
! mass & energy loss to the ocean
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
+ 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 (jk) = MAX( 0._wp, zh_s(jk) - zdum )
+ h_s_1d(ji) = MAX( 0._wp, h_s_1d(ji) - 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
@@ -521,8 +518,8 @@ 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(0) = zh_i(0) + dh_snowice(ji)
- zdeltah = 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(0) = zh_i_old(0) + dh_snowice(ji)
@@ -530,11 +527,11 @@ CONTAINS
!
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
+ 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)
+ 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 )
+ zdeltah = MAX( 0._wp, zdeltah - zdum )
END DO
!
!
@@ -561,12 +558,13 @@ CONTAINS
! Remapping of ice enthalpy on a regular grid
!--------------------------------------------
- e_i_1d(ji,:) = ice_ent( zh_i_old(:), ze_i_old(:) )
+ e_i_1d(ji,:) = ice_ent1( zh_i_old(:), ze_i_old(:) )
END DO ! npti
-
-
-
+ ! ! ================== !
+ ! ! 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
@@ -576,7 +574,7 @@ CONTAINS
END SUBROUTINE ice_thd_dh
- PURE FUNCTION snw_ent( ph_old, pe_old )
+ FUNCTION snw_ent( ph_old, pe_old )
!!-------------------------------------------------------------------
!! *** ROUTINE snw_ent ***
!!
@@ -659,9 +657,9 @@ CONTAINS
END FUNCTION snw_ent
- PURE FUNCTION ice_ent( ph_old, pe_old )
+ FUNCTION ice_ent1( ph_old, pe_old )
!!-------------------------------------------------------------------
- !! *** ROUTINE ice_ent ***
+ !! *** ROUTINE ice_ent1 ***
!!
!! ** Purpose :
!! This routine computes new vertical grids in the ice,
@@ -685,7 +683,7 @@ CONTAINS
!! 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)
+ 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
@@ -735,7 +733,7 @@ CONTAINS
! 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
+ 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 --- !
@@ -745,7 +743,7 @@ CONTAINS
! & ( SUM( pe_new(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_old(ji,0:nlay_i+1) ) )
- END FUNCTION ice_ent
+ END FUNCTION ice_ent1
#else
!!----------------------------------------------------------------------
diff --git a/src/ICE/icethd_do.F90 b/src/ICE/icethd_do.F90
index 6867b76a3aa22550c12f601cad9335117d8aa1ce..063c3e903702c89f77810d94beae028aab95620a 100644
--- a/src/ICE/icethd_do.F90
+++ b/src/ICE/icethd_do.F90
@@ -102,11 +102,11 @@ CONTAINS
IF( ln_icediachk ) CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft )
IF( ln_icediachk ) CALL ice_cons2D ( 0, 'icethd_do', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft )
- at_i(A2D(0)) = SUM( a_i(A2D(0),:), dim=3 )
!------------------------------------------------------------------------------!
! 1) Compute thickness, salinity, enthalpy, age, area and volume of new ice
!------------------------------------------------------------------------------!
! it occurs if cooling
+ at_i(A2D(0)) = SUM( a_i(A2D(0),:), dim=3 )
! Identify grid points where new ice forms
npti = 0 ; nptidx(:) = 0
@@ -160,7 +160,9 @@ CONTAINS
zs_newice(1:npti) = 2.3
END SELECT
!
- !
+ ! ! ==================== !
+ ! ! Start main loop here !
+ ! ! ==================== !
DO ji = 1, npti
! Keep old ice areas and volume in memory
@@ -273,12 +275,14 @@ CONTAINS
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(:) )
+ e_i_2d(ji,:,jl) = ice_ent2( zh_i_old(:), ze_i_old(:) )
END DO
END DO ! npti
-
-
+ ! ! ================== !
+ ! ! End main loop here !
+ ! ! ================== !
+ !
! Change units for e_i
DO jl = 1, jpl
DO jk = 1, nlay_i
@@ -393,9 +397,9 @@ CONTAINS
ENDIF
END SUBROUTINE ice_thd_frazil
- PURE FUNCTION ice_ent( ph_old, pe_old )
+ FUNCTION ice_ent2( ph_old, pe_old )
!!-------------------------------------------------------------------
- !! *** ROUTINE ice_ent ***
+ !! *** ROUTINE ice_ent2 ***
!!
!! ** Purpose :
!! This routine computes new vertical grids in the ice,
@@ -419,7 +423,7 @@ CONTAINS
!! 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)
+ REAL(wp), DIMENSION(1:nlay_i) :: ice_ent2 ! new enthlapies (J.m-3, remapped)
!
INTEGER :: ji ! dummy loop indices
INTEGER :: jk0, jk1 ! old/new layer indices
@@ -468,8 +472,8 @@ CONTAINS
! 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
+ zswitch = MAX( 0._wp , SIGN( 1._wp , zhnew - epsi20 ) )
+ ice_ent2(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 --- !
@@ -479,7 +483,7 @@ CONTAINS
! & ( SUM( pe_new(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_old(ji,0:nlay_i+1) ) )
- END FUNCTION ice_ent
+ END FUNCTION ice_ent2
SUBROUTINE ice_thd_do_init