New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

icethd_dh.F90 in NEMO/branches/2020/tickets_icb_1900/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/tickets_icb_1900/src/ICE/icethd_dh.F90 @ 14016

Last change on this file since 14016 was 14016, checked in by mathiot, 3 years ago
ticket 1900: upgrade to trunk@r14015 (head trunk at 16h27)
Property svn:keywords set to `Id`
File size: 36.4 KB

Rev	Line
[8586]	1	MODULE icethd_dh
	2	!!======================================================================
	3	!! * MODULE icethd_dh *
	4	!! seaice : thermodynamic growth and melt
	5	!!======================================================================
[9604]	6	!! History : ! 2003-05 (M. Vancoppenolle) Original code in 1D
	7	!! ! 2005-06 (M. Vancoppenolle) 3D version
	8	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
[8586]	9	!!----------------------------------------------------------------------
[9570]	10	#if defined key_si3
[8586]	11	!!----------------------------------------------------------------------
[9570]	12	!! 'key_si3' SI3 sea-ice model
[8586]	13	!!----------------------------------------------------------------------
	14	!! ice_thd_dh : vertical sea-ice growth and melt
[13899]	15	!!----------------------------------------------------------------------
[8586]	16	USE dom_oce ! ocean space and time domain
	17	USE phycst ! physical constants
	18	USE ice ! sea-ice: variables
	19	USE ice1D ! sea-ice: thermodynamics variables
	20	USE icethd_sal ! sea-ice: salinity profiles
[13899]	21	USE icevar ! for CALL ice_var_snwblow
[8586]	22	!
	23	USE in_out_manager ! I/O manager
	24	USE lib_mpp ! MPP library
	25	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
	26
	27	IMPLICIT NONE
	28	PRIVATE
	29
	30	PUBLIC ice_thd_dh ! called by ice_thd
	31
	32	!!----------------------------------------------------------------------
[9598]	33	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	34	!! $Id$
[10068]	35	!! Software governed by the CeCILL license (see ./LICENSE)
[8586]	36	!!----------------------------------------------------------------------
	37	CONTAINS
	38
	39	SUBROUTINE ice_thd_dh
	40	!!------------------------------------------------------------------
	41	!! * ROUTINE ice_thd_dh *
	42	!!
[9274]	43	!! ** Purpose : compute ice and snow thickness changes due to growth/melting
[8586]	44	!!
	45	!! ** Method : Ice/Snow surface melting arises from imbalance in surface fluxes
[9274]	46	!! Bottom accretion/ablation arises from flux budget
	47	!! Snow thickness can increase by precipitation and decrease by sublimation
	48	!! If snow load excesses Archmiede limit, snow-ice is formed by
	49	!! the flooding of sea-water in the snow
[8586]	50	!!
[9274]	51	!! - Compute available flux of heat for surface ablation
	52	!! - Compute snow and sea ice enthalpies
	53	!! - Surface ablation and sublimation
	54	!! - Bottom accretion/ablation
	55	!! - Snow ice formation
[8586]	56	!!
[14016]	57	!! ** Note : h=max(0,h+dh) are often used to ensure positivity of h.
	58	!! very small negative values can occur otherwise (e.g. -1.e-20)
	59	!!
[8586]	60	!! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
	61	!! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646
	62	!! Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let.
	63	!! Vancoppenolle et al.,2009, Ocean Modelling
	64	!!------------------------------------------------------------------
	65	INTEGER :: ji, jk ! dummy loop indices
	66	INTEGER :: iter ! local integer
	67
	68	REAL(wp) :: ztmelts ! local scalar
	69	REAL(wp) :: zdum
	70	REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment
	71	REAL(wp) :: zswi1 ! switch for computation of bottom salinity
	72	REAL(wp) :: zswi12 ! switch for computation of bottom salinity
	73	REAL(wp) :: zswi2 ! switch for computation of bottom salinity
	74	REAL(wp) :: zgrr ! bottom growth rate
	75	REAL(wp) :: zt_i_new ! bottom formation temperature
[12489]	76	REAL(wp) :: z1_rho ! 1/(rhos+rho0-rhoi)
[8586]	77
	78	REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
	79	REAL(wp) :: zEi ! specific enthalpy of sea ice (J/kg)
	80	REAL(wp) :: zEw ! specific enthalpy of exchanged water (J/kg)
	81	REAL(wp) :: zdE ! specific enthalpy difference (J/kg)
	82	REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean
	83
[9922]	84	REAL(wp), DIMENSION(jpij) :: zq_top ! heat for surface ablation (J.m-2)
	85	REAL(wp), DIMENSION(jpij) :: zq_bot ! heat for bottom ablation (J.m-2)
[8586]	86	REAL(wp), DIMENSION(jpij) :: zq_rema ! remaining heat at the end of the routine (J.m-2)
	87	REAL(wp), DIMENSION(jpij) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2)
	88	REAL(wp), DIMENSION(jpij) :: zevap_rema ! remaining mass flux from sublimation (kg.m-2)
[14016]	89	REAL(wp), DIMENSION(jpij) :: zdeltah
	90	REAL(wp), DIMENSION(jpij) :: zsnw ! distribution of snow after wind blowing
[8586]	91
[14016]	92	INTEGER , DIMENSION(jpij,nlay_i) :: icount ! number of layers vanishing by melting
	93	REAL(wp), DIMENSION(jpij,0:nlay_i+1) :: zh_i ! ice layer thickness (m)
	94	REAL(wp), DIMENSION(jpij,0:nlay_s ) :: zh_s ! snw layer thickness (m)
	95	REAL(wp), DIMENSION(jpij,0:nlay_s ) :: ze_s ! snw layer enthalpy (J.m-3)
[8586]	96
	97	REAL(wp) :: zswitch_sal
	98
	99	INTEGER :: num_iter_max ! Heat conservation
	100	!!------------------------------------------------------------------
	101
	102	! Discriminate between time varying salinity and constant
	103	SELECT CASE( nn_icesal ) ! varying salinity or not
	104	CASE( 1 , 3 ) ; zswitch_sal = 0._wp ! prescribed salinity profile
	105	CASE( 2 ) ; zswitch_sal = 1._wp ! varying salinity profile
	106	END SELECT
	107
[14016]	108	! initialize ice layer thicknesses and enthalpies
	109	eh_i_old(1:npti,0:nlay_i+1) = 0._wp
[8586]	110	h_i_old (1:npti,0:nlay_i+1) = 0._wp
[14016]	111	zh_i (1:npti,0:nlay_i+1) = 0._wp
[8586]	112	DO jk = 1, nlay_i
	113	DO ji = 1, npti
[14016]	114	eh_i_old(ji,jk) = h_i_1d(ji) * r1_nlay_i * e_i_1d(ji,jk)
[8586]	115	h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i
[14016]	116	zh_i (ji,jk) = h_i_1d(ji) * r1_nlay_i
[8586]	117	END DO
	118	END DO
	119	!
[14016]	120	! initialize snw layer thicknesses and enthalpies
	121	zh_s(1:npti,0) = 0._wp
	122	ze_s(1:npti,0) = 0._wp
	123	DO jk = 1, nlay_s
	124	DO ji = 1, npti
	125	zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
	126	ze_s(ji,jk) = e_s_1d(ji,jk)
	127	END DO
	128	END DO
	129	!
[9604]	130	! ! ============================================== !
	131	! ! Available heat for surface and bottom ablation !
	132	! ! ============================================== !
[8813]	133	!
[10534]	134	IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN
[8813]	135	!
	136	DO ji = 1, npti
[12489]	137	zq_top(ji) = MAX( 0._wp, qml_ice_1d(ji) * rDt_ice )
[8813]	138	END DO
	139	!
[10534]	140	ELSE
[8813]	141	!
	142	DO ji = 1, npti
[9916]	143	zdum = qns_ice_1d(ji) + qsr_ice_1d(ji) - qtr_ice_top_1d(ji) - qcn_ice_top_1d(ji)
[9274]	144	qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
[12489]	145	zq_top(ji) = MAX( 0._wp, qml_ice_1d(ji) * rDt_ice )
[8813]	146	END DO
	147	!
[10534]	148	ENDIF
[8813]	149	!
[8586]	150	DO ji = 1, npti
[13899]	151	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)
[12489]	152	zq_bot(ji) = MAX( 0._wp, zf_tt(ji) * rDt_ice )
[8586]	153	END DO
	154
[9274]	155	! ! ============ !
	156	! ! Snow !
	157	! ! ============ !
	158	!
	159	! Internal melting
	160	! ----------------
	161	! IF snow temperature is above freezing point, THEN snow melts (should not happen but sometimes it does)
[9271]	162	DO jk = 1, nlay_s
[8586]	163	DO ji = 1, npti
[9274]	164	IF( t_s_1d(ji,jk) > rt0 ) THEN
[14016]	165	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
	166	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! mass flux
[9271]	167	! updates
[14016]	168	dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk)
	169	h_s_1d (ji) = MAX( 0._wp, h_s_1d (ji) - zh_s(ji,jk) )
[9274]	170	zh_s (ji,jk) = 0._wp
[14016]	171	ze_s (ji,jk) = 0._wp
[9271]	172	END IF
[8586]	173	END DO
[9274]	174	END DO
[8586]	175
[9274]	176	! Snow precipitation
	177	!-------------------
[14016]	178	CALL ice_var_snwblow( 1._wp - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing
[8586]	179
	180	DO ji = 1, npti
[9274]	181	IF( sprecip_1d(ji) > 0._wp ) THEN
[14016]	182	zh_s(ji,0) = zsnw(ji) * sprecip_1d(ji) * rDt_ice * r1_rhos / at_i_1d(ji) ! thickness of precip
	183	ze_s(ji,0) = MAX( 0._wp, - qprec_ice_1d(ji) ) ! enthalpy of the precip (>0, J.m-3)
[9274]	184	!
[14016]	185	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)
	186	wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos * zh_s(ji,0) * a_i_1d(ji) * r1_Dt_ice ! mass flux, <0
[9274]	187	!
	188	! update thickness
[14016]	189	h_s_1d(ji) = h_s_1d(ji) + zh_s(ji,0)
[9274]	190	ENDIF
[8586]	191	END DO
	192
[9274]	193	! Snow melting
	194	! ------------
[14016]	195	! If heat still available (zq_top > 0)
	196	! then all snw precip has been melted and we need to melt more snow
	197	DO jk = 0, nlay_s
[8586]	198	DO ji = 1, npti
[9922]	199	IF( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN
[9274]	200	!
[14016]	201	rswitch = MAX( 0._wp , SIGN( 1._wp , ze_s(ji,jk) - epsi20 ) )
	202	zdum = - rswitch * zq_top(ji) / MAX( ze_s(ji,jk), epsi20 ) ! thickness change
	203	zdum = MAX( zdum , - zh_s(ji,jk) ) ! bound melting
[9274]	204
[14016]	205	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)
	206	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
[9274]	207
	208	! updates available heat + thickness
[14016]	209	dh_s_mlt(ji) = dh_s_mlt(ji) + zdum
	210	zq_top (ji) = MAX( 0._wp , zq_top (ji) + zdum * ze_s(ji,jk) )
	211	h_s_1d (ji) = MAX( 0._wp , h_s_1d (ji) + zdum )
	212	zh_s (ji,jk) = MAX( 0._wp , zh_s (ji,jk) + zdum )
	213	!!$ IF( zh_s(ji,jk) == 0._wp ) ze_s(ji,jk) = 0._wp
[9274]	214	!
	215	ENDIF
[8586]	216	END DO
	217	END DO
	218
[9274]	219	! Snow sublimation
	220	!-----------------
[8586]	221	! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
[8885]	222	! comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
[14016]	223	zdeltah (1:npti) = 0._wp ! total snow thickness that sublimates, < 0
	224	zevap_rema(1:npti) = 0._wp
[8586]	225	DO ji = 1, npti
[9274]	226	IF( evap_ice_1d(ji) > 0._wp ) THEN
[14016]	227	zdeltah (ji) = MAX( - evap_ice_1d(ji) * r1_rhos * rDt_ice, - h_s_1d(ji) ) ! amount of snw that sublimates, < 0
	228	zevap_rema(ji) = MAX( 0._wp, evap_ice_1d(ji) * rDt_ice + zdeltah(ji) * rhos ) ! remaining evap in kg.m-2 (used for ice sublimation later on)
[9274]	229	ENDIF
[8586]	230	END DO
	231
[14016]	232	DO jk = 0, nlay_s
	233	DO ji = 1, npti
	234	zdum = MAX( -zh_s(ji,jk), zdeltah(ji) ) ! snow layer thickness that sublimates, < 0
	235	!
	236	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
	237	wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux by sublimation
[8586]	238
[14016]	239	! update thickness
	240	h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdum )
	241	zh_s (ji,jk) = MAX( 0._wp , zh_s (ji,jk) + zdum )
	242	!!$ IF( zh_s(ji,jk) == 0._wp ) ze_s(ji,jk) = 0._wp
	243
	244	! update sublimation left
	245	zdeltah(ji) = MIN( zdeltah(ji) - zdum, 0._wp )
[8586]	246	END DO
	247	END DO
[14016]	248
	249	!
[9274]	250	! ! ============ !
	251	! ! Ice !
	252	! ! ============ !
[8586]	253
[9274]	254	! Surface ice melting
	255	!--------------------
[8586]	256	DO jk = 1, nlay_i
	257	DO ji = 1, npti
[9935]	258	ztmelts = - rTmlt * sz_i_1d(ji,jk) ! Melting point of layer k [C]
[8586]	259
[9274]	260	IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting
[8586]	261
[9935]	262	zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of layer k [J/kg, <0]
[14016]	263	zdE = 0._wp ! Specific enthalpy difference (J/kg, <0)
	264	! set up at 0 since no energy is needed to melt water...(it is already melted)
	265	zdum = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing
	266	! this should normally not happen, but sometimes, heat diffusion leads to this
	267	zfmdt = - zdum * rhoi ! Recompute mass flux [kg/m2, >0]
	268	!
	269	dh_i_itm(ji) = dh_i_itm(ji) + zdum ! Cumulate internal melting
	270	!
	271	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
	272	! ice enthalpy zEi is "sent" to the ocean
	273	wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux
	274	sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux
	275	! using s_i_1d and not sz_i_1d(jk) is ok
[9274]	276	ELSE !-- Surface melting
[8586]	277
[9935]	278	zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of layer k [J/kg, <0]
[8586]	279	zEw = rcp * ztmelts ! Specific enthalpy of resulting meltwater [J/kg, <0]
	280	zdE = zEi - zEw ! Specific enthalpy difference < 0
	281
[9922]	282	zfmdt = - zq_top(ji) / zdE ! Mass flux to the ocean [kg/m2, >0]
[8586]	283
[14016]	284	zdum = - zfmdt * r1_rhoi ! Melt of layer jk [m, <0]
[8586]	285
[14016]	286	zdum = MIN( 0._wp , MAX( zdum , - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0]
	287
	288	zq_top(ji) = MAX( 0._wp , zq_top(ji) - zdum * rhoi * zdE ) ! update available heat
[8586]	289
[14016]	290	dh_i_sum(ji) = dh_i_sum(ji) + zdum ! Cumulate surface melt
[8586]	291
[14016]	292	zfmdt = - rhoi * zdum ! Recompute mass flux [kg/m2, >0]
[8586]	293
	294	zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0]
	295
[14016]	296	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux [W.m-2], < 0
	297	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
	298	wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux
	299	sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux >0
	300	! using s_i_1d and not sz_i_1d(jk) is ok)
[8586]	301	END IF
[14016]	302	! update thickness
	303	zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum )
	304	h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum )
	305	!
	306	! update heat content (J.m-2) and layer thickness
	307	eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk)
	308	h_i_old (ji,jk) = h_i_old (ji,jk) + zdum
	309	!
	310	!
[9274]	311	! Ice sublimation
	312	! ---------------
[14016]	313	zdum = MAX( - zh_i(ji,jk) , - zevap_rema(ji) * r1_rhoi )
	314	!
	315	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
	316	wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux > 0
	317	sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux >0
	318	! clem: flux is sent to the ocean for simplicity
	319	! but salt should remain in the ice except
	320	! if all ice is melted. => must be corrected
	321	! update remaining mass flux and thickness
	322	zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi
	323	zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum )
	324	h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum )
	325	dh_i_sub(ji) = dh_i_sub(ji) + zdum
[8586]	326
[14016]	327	! update heat content (J.m-2) and layer thickness
	328	eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk)
	329	h_i_old (ji,jk) = h_i_old (ji,jk) + zdum
[9274]	330
[8586]	331	! record which layers have disappeared (for bottom melting)
	332	! => icount=0 : no layer has vanished
	333	! => icount=5 : 5 layers have vanished
[14016]	334	rswitch = MAX( 0._wp , SIGN( 1._wp , - zh_i(ji,jk) ) )
[8586]	335	icount(ji,jk) = NINT( rswitch )
	336
	337	END DO
	338	END DO
[9274]	339
[8586]	340	! remaining "potential" evap is sent to ocean
	341	DO ji = 1, npti
[12489]	342	wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_Dt_ice ! <=0 (net evap for the ocean in kg.m-2.s-1)
[8586]	343	END DO
	344
	345
[9274]	346	! Ice Basal growth
[8586]	347	!------------------
	348	! Basal growth is driven by heat imbalance at the ice-ocean interface,
[9916]	349	! between the inner conductive flux (qcn_ice_bot), from the open water heat flux
[9913]	350	! (fhld) and the sensible ice-ocean flux (qsb_ice_bot).
[9916]	351	! qcn_ice_bot is positive downwards. qsb_ice_bot and fhld are positive to the ice
[8586]	352
	353	! If salinity varies in time, an iterative procedure is required, because
	354	! the involved quantities are inter-dependent.
[9750]	355	! Basal growth (dh_i_bog) depends upon new ice specific enthalpy (zEi),
	356	! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bog)
[8586]	357	! -> need for an iterative procedure, which converges quickly
	358
	359	num_iter_max = 1
	360	IF( nn_icesal == 2 ) num_iter_max = 5 ! salinity varying in time
	361
	362	DO ji = 1, npti
	363	IF( zf_tt(ji) < 0._wp ) THEN
[9274]	364	DO iter = 1, num_iter_max ! iterations
[8586]	365
	366	! New bottom ice salinity (Cox & Weeks, JGR88 )
	367	!--- zswi1 if dh/dt < 2.0e-8
	368	!--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
	369	!--- zswi2 if dh/dt > 3.6e-7
[12489]	370	zgrr = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_Dt_ice , epsi10 ) )
[9274]	371	zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
	372	zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
	373	zswi1 = 1. - zswi2 * zswi12
	374	zfracs = MIN( zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) &
	375	& + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 )
[8586]	376
[14016]	377	s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji) ! New ice salinity
[8586]	378
[14016]	379	ztmelts = - rTmlt * s_i_new(ji) ! New ice melting point (C)
[9274]	380
[14016]	381	zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
[8586]	382
[14016]	383	zEi = rcpi * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0)
	384	& - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp * ztmelts
[8586]	385
[14016]	386	zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)
[8586]	387
[14016]	388	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
[8586]	389
[14016]	390	dh_i_bog(ji) = rDt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) )
[8586]	391
	392	END DO
	393	! Contribution to Energy and Salt Fluxes
[14016]	394	zfmdt = - rhoi * dh_i_bog(ji) ! Mass flux x time step (kg/m2, < 0)
[8586]	395
[14016]	396	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
	397	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
	398	wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * dh_i_bog(ji) * a_i_1d(ji) * r1_Dt_ice ! Mass flux, <0
	399	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
[8586]	400
[14016]	401	! update thickness
	402	zh_i(ji,nlay_i+1) = zh_i(ji,nlay_i+1) + dh_i_bog(ji)
	403	h_i_1d(ji) = h_i_1d(ji) + dh_i_bog(ji)
[8586]	404
	405	! update heat content (J.m-2) and layer thickness
[9935]	406	eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi)
[9750]	407	h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bog(ji)
[8586]	408
	409	ENDIF
	410
	411	END DO
	412
[9274]	413	! Ice Basal melt
	414	!---------------
[8586]	415	DO jk = nlay_i, 1, -1
	416	DO ji = 1, npti
	417	IF( zf_tt(ji) > 0._wp .AND. jk > icount(ji,jk) ) THEN ! do not calculate where layer has already disappeared by surface melting
	418
[9935]	419	ztmelts = - rTmlt * sz_i_1d(ji,jk) ! Melting point of layer jk (C)
[8586]	420
[9274]	421	IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting
[8586]	422
[14016]	423	zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0)
	424	zdE = 0._wp ! Specific enthalpy difference (J/kg, <0)
	425	! set up at 0 since no energy is needed to melt water...(it is already melted)
	426	zdum = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing
	427	! this should normally not happen, but sometimes, heat diffusion leads to this
	428	dh_i_itm (ji) = dh_i_itm(ji) + zdum
	429	!
	430	zfmdt = - zdum * rhoi ! Mass flux x time step > 0
	431	!
	432	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
	433	! ice enthalpy zEi is "sent" to the ocean
	434	wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux
	435	sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux
	436	! using s_i_1d and not sz_i_1d(jk) is ok
[9274]	437	ELSE !-- Basal melting
[8586]	438
[14016]	439	zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0)
	440	zEw = rcp * ztmelts ! Specific enthalpy of meltwater (J/kg, <0)
	441	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
[8586]	442
[14016]	443	zfmdt = - zq_bot(ji) / zdE ! Mass flux x time step (kg/m2, >0)
[8586]	444
[14016]	445	zdum = - zfmdt * r1_rhoi ! Gross thickness change
[8586]	446
[14016]	447	zdum = MIN( 0._wp , MAX( zdum, - zh_i(ji,jk) ) ) ! bound thickness change
[8586]	448
[14016]	449	zq_bot(ji) = MAX( 0._wp , zq_bot(ji) - zdum * rhoi * zdE ) ! update available heat. MAX is necessary for roundup errors
[8586]	450
[14016]	451	dh_i_bom(ji) = dh_i_bom(ji) + zdum ! Update basal melt
[8586]	452
[14016]	453	zfmdt = - zdum * rhoi ! Mass flux x time step > 0
[8586]	454
[14016]	455	zQm = zfmdt * zEw ! Heat exchanged with ocean
[8586]	456
[14016]	457	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
	458	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
	459	wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux
	460	sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux
	461	! using s_i_1d and not sz_i_1d(jk) is ok
[8586]	462	ENDIF
[14016]	463	! update thickness
	464	zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum )
	465	h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum )
	466	!
	467	! update heat content (J.m-2) and layer thickness
	468	eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk)
	469	h_i_old (ji,jk) = h_i_old (ji,jk) + zdum
[8586]	470	ENDIF
	471	END DO
	472	END DO
	473
[14016]	474	! Remove snow if ice has melted entirely
	475	! --------------------------------------
	476	DO jk = 0, nlay_s
	477	DO ji = 1,npti
	478	IF( h_i_1d(ji) == 0._wp ) THEN
	479	! mass & energy loss to the ocean
	480	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
	481	wfx_res_1d(ji) = wfx_res_1d(ji) + rhos * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! mass flux
[8586]	482
[14016]	483	! update thickness and energy
	484	h_s_1d(ji) = 0._wp
	485	ze_s (ji,jk) = 0._wp
	486	zh_s (ji,jk) = 0._wp
	487	ENDIF
	488	END DO
[8586]	489	END DO
[14016]	490
	491	! Snow load on ice
	492	! -----------------
	493	! When snow load exceeds Archimede's limit and sst is positive,
	494	! snow-ice formation (next bloc) can lead to negative ice enthalpy.
	495	! Therefore we consider here that this excess of snow falls into the ocean
	496	zdeltah(1:npti) = h_s_1d(1:npti) + h_i_1d(1:npti) * (rhoi-rho0) * r1_rhos
	497	DO jk = 0, nlay_s
	498	DO ji = 1, npti
	499	IF( zdeltah(ji) > 0._wp .AND. sst_1d(ji) > 0._wp ) THEN
	500	! snow layer thickness that falls into the ocean
	501	zdum = MIN( zdeltah(ji) , zh_s(ji,jk) )
	502	! mass & energy loss to the ocean
	503	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
	504	wfx_res_1d(ji) = wfx_res_1d(ji) + rhos * zdum * a_i_1d(ji) * r1_Dt_ice ! mass flux
	505	! update thickness and energy
	506	h_s_1d(ji) = MAX( 0._wp, h_s_1d(ji) - zdum )
	507	zh_s (ji,jk) = MAX( 0._wp, zh_s(ji,jk) - zdum )
	508	! update snow thickness that still has to fall
	509	zdeltah(ji) = MAX( 0._wp, zdeltah(ji) - zdum )
	510	ENDIF
	511	END DO
	512	END DO
	513
[9274]	514	! Snow-Ice formation
	515	! ------------------
[14016]	516	! When snow load exceeds Archimede's limit, snow-ice interface goes down under sea-level,
	517	! flooding of seawater transforms snow into ice. Thickness that is transformed is dh_snowice (positive for the ice)
[12489]	518	z1_rho = 1._wp / ( rhos+rho0-rhoi )
[14016]	519	zdeltah(1:npti) = 0._wp
[8586]	520	DO ji = 1, npti
	521	!
[14016]	522	dh_snowice(ji) = MAX( 0._wp , ( rhos * h_s_1d(ji) + (rhoi-rho0) * h_i_1d(ji) ) * z1_rho )
[8586]	523
	524	h_i_1d(ji) = h_i_1d(ji) + dh_snowice(ji)
	525	h_s_1d(ji) = h_s_1d(ji) - dh_snowice(ji)
	526
	527	! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
[9935]	528	zfmdt = ( rhos - rhoi ) * dh_snowice(ji) ! <0
[8586]	529	zEw = rcp * sst_1d(ji)
	530	zQm = zfmdt * zEw
	531
[14016]	532	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux
	533	sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Salt flux
[8586]	534
	535	! Case constant salinity in time: virtual salt flux to keep salinity constant
[8637]	536	IF( nn_icesal /= 2 ) THEN
[14016]	537	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
	538	& - s_i_1d(ji) * dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_Dt_ice ! and get rn_icesal from the ocean
[8586]	539	ENDIF
	540
[9274]	541	! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume
[14016]	542	wfx_sni_1d (ji) = wfx_sni_1d (ji) - dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_Dt_ice
	543	wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + dh_snowice(ji) * rhos * a_i_1d(ji) * r1_Dt_ice
[8586]	544
[14016]	545	! update thickness
	546	zh_i(ji,0) = zh_i(ji,0) + dh_snowice(ji)
	547	zdeltah(ji) = dh_snowice(ji)
	548
[8586]	549	! update heat content (J.m-2) and layer thickness
	550	h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
[14016]	551	eh_i_old(ji,0) = eh_i_old(ji,0) + zfmdt * zEw ! 1st part (sea water enthalpy)
	552
[8586]	553	END DO
	554	!
[14016]	555	DO jk = nlay_s, 0, -1 ! flooding of snow starts from the base
	556	DO ji = 1, npti
	557	zdum = MIN( zdeltah(ji), zh_s(ji,jk) ) ! amount of snw that floods, > 0
	558	zh_s(ji,jk) = MAX( 0._wp, zh_s(ji,jk) - zdum ) ! remove some snow thickness
	559	eh_i_old(ji,0) = eh_i_old(ji,0) + zdum * ze_s(ji,jk) ! 2nd part (snow enthalpy)
	560	! update dh_snowice
	561	zdeltah(ji) = MAX( 0._wp, zdeltah(ji) - zdum )
	562	END DO
[8586]	563	END DO
[14016]	564	!
	565	!
	566	!!$ ! --- Update snow diags --- !
	567	!!$ !!clem: this is wrong. dh_s_tot is not used anyway
	568	!!$ DO ji = 1, npti
	569	!!$ dh_s_tot(ji) = dh_s_tot(ji) + dh_s_mlt(ji) + zdeltah(ji) + zdh_s_sub(ji) - dh_snowice(ji)
	570	!!$ END DO
	571	!
	572	!
	573	! Remapping of snw enthalpy on a regular grid
	574	!--------------------------------------------
	575	CALL snw_ent( zh_s, ze_s, e_s_1d )
	576
	577	! recalculate t_s_1d from e_s_1d
[8586]	578	DO jk = 1, nlay_s
	579	DO ji = 1,npti
[14016]	580	IF( h_s_1d(ji) > 0._wp ) THEN
	581	t_s_1d(ji,jk) = rt0 + ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi )
	582	ELSE
	583	t_s_1d(ji,jk) = rt0
	584	ENDIF
[8586]	585	END DO
	586	END DO
	587
[14016]	588	! Note: remapping of ice enthalpy is done in icethd.F90
	589
[10786]	590	! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 ---
	591	WHERE( h_i_1d(1:npti) == 0._wp )
[14016]	592	a_i_1d (1:npti) = 0._wp
	593	h_s_1d (1:npti) = 0._wp
	594	t_su_1d(1:npti) = rt0
[10786]	595	END WHERE
[14016]	596
[8586]	597	END SUBROUTINE ice_thd_dh
	598
[14016]	599	SUBROUTINE snw_ent( ph_old, pe_old, pe_new )
	600	!!-------------------------------------------------------------------
	601	!! * ROUTINE snw_ent *
	602	!!
	603	!! ** Purpose :
	604	!! This routine computes new vertical grids in the snow,
	605	!! and consistently redistributes temperatures.
	606	!! Redistribution is made so as to ensure to energy conservation
	607	!!
	608	!!
	609	!! ** Method : linear conservative remapping
	610	!!
	611	!! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses
	612	!! 2) linear remapping on the new layers
	613	!!
	614	!! ------------ cum0(0) ------------- cum1(0)
	615	!! NEW -------------
	616	!! ------------ cum0(1) ==> -------------
	617	!! ... -------------
	618	!! ------------ -------------
	619	!! ------------ cum0(nlay_s+1) ------------- cum1(nlay_s)
	620	!!
	621	!!
	622	!! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005
	623	!!-------------------------------------------------------------------
	624	REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in ) :: ph_old ! old thicknesses (m)
	625	REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in ) :: pe_old ! old enthlapies (J.m-3)
	626	REAL(wp), DIMENSION(jpij,1:nlay_s), INTENT(inout) :: pe_new ! new enthlapies (J.m-3, remapped)
	627	!
	628	INTEGER :: ji ! dummy loop indices
	629	INTEGER :: jk0, jk1 ! old/new layer indices
	630	!
	631	REAL(wp), DIMENSION(jpij,0:nlay_s+1) :: zeh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces
	632	REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeh_cum1, zh_cum1 ! new cumulative enthlapies and layers interfaces
	633	REAL(wp), DIMENSION(jpij) :: zhnew ! new layers thicknesses
	634	!!-------------------------------------------------------------------
	635
	636	!--------------------------------------------------------------------------
	637	! 1) Cumulative integral of old enthalpy * thickness and layers interfaces
	638	!--------------------------------------------------------------------------
	639	zeh_cum0(1:npti,0) = 0._wp
	640	zh_cum0 (1:npti,0) = 0._wp
	641	DO jk0 = 1, nlay_s+1
	642	DO ji = 1, npti
	643	zeh_cum0(ji,jk0) = zeh_cum0(ji,jk0-1) + pe_old(ji,jk0-1) * ph_old(ji,jk0-1)
	644	zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + ph_old(ji,jk0-1)
	645	END DO
	646	END DO
	647
	648	!------------------------------------
	649	! 2) Interpolation on the new layers
	650	!------------------------------------
	651	! new layer thickesses
	652	DO ji = 1, npti
	653	zhnew(ji) = SUM( ph_old(ji,0:nlay_s) ) * r1_nlay_s
	654	END DO
	655
	656	! new layers interfaces
	657	zh_cum1(1:npti,0) = 0._wp
	658	DO jk1 = 1, nlay_s
	659	DO ji = 1, npti
	660	zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji)
	661	END DO
	662	END DO
	663
	664	zeh_cum1(1:npti,0:nlay_s) = 0._wp
	665	! new cumulative q*h => linear interpolation
	666	DO jk0 = 1, nlay_s+1
	667	DO jk1 = 1, nlay_s-1
	668	DO ji = 1, npti
	669	IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN
	670	zeh_cum1(ji,jk1) = ( zeh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1 ) ) + &
	671	& zeh_cum0(ji,jk0 ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) ) &
	672	& / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) )
	673	ENDIF
	674	END DO
	675	END DO
	676	END DO
	677	! to ensure that total heat content is strictly conserved, set:
	678	zeh_cum1(1:npti,nlay_s) = zeh_cum0(1:npti,nlay_s+1)
	679
	680	! new enthalpies
	681	DO jk1 = 1, nlay_s
	682	DO ji = 1, npti
	683	rswitch = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) )
	684	pe_new(ji,jk1) = rswitch * ( zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 )
	685	END DO
	686	END DO
	687
	688	END SUBROUTINE snw_ent
	689
	690
[8586]	691	#else
	692	!!----------------------------------------------------------------------
[9570]	693	!! Default option NO SI3 sea-ice model
[8586]	694	!!----------------------------------------------------------------------
	695	#endif
	696
	697	!!======================================================================
	698	END MODULE icethd_dh

Note: See TracBrowser for help on using the repository browser.

Download in other formats: