Skip to content
Snippets Groups Projects
Commit c06e429d authored by Clement Rousset's avatar Clement Rousset
Browse files

forgotten commits

parent 7437c60b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/ICE/icethd_dh.F90
+ 58
− 60
View file @ c06e429d
......@@ -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
!!----------------------------------------------------------------------
......
src/ICE/icethd_do.F90
+ 15
− 11
View file @ c06e429d
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment