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/2019/dev_r11842_SI3-10_EAP/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icethd_dh.F90 @ 13662

Last change on this file since 13662 was 13662, checked in by clem, 4 years ago
update to almost r4.0.4
Property svn:keywords set to `Id`
File size: 34.1 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
[13662]	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
[13662]	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	!!
	57	!! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
	58	!! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646
	59	!! Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let.
	60	!! Vancoppenolle et al.,2009, Ocean Modelling
	61	!!------------------------------------------------------------------
	62	INTEGER :: ji, jk ! dummy loop indices
	63	INTEGER :: iter ! local integer
	64
	65	REAL(wp) :: ztmelts ! local scalar
	66	REAL(wp) :: zdum
	67	REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment
	68	REAL(wp) :: zswi1 ! switch for computation of bottom salinity
	69	REAL(wp) :: zswi12 ! switch for computation of bottom salinity
	70	REAL(wp) :: zswi2 ! switch for computation of bottom salinity
	71	REAL(wp) :: zgrr ! bottom growth rate
	72	REAL(wp) :: zt_i_new ! bottom formation temperature
[9935]	73	REAL(wp) :: z1_rho ! 1/(rhos+rau0-rhoi)
[8586]	74
	75	REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
	76	REAL(wp) :: zEi ! specific enthalpy of sea ice (J/kg)
	77	REAL(wp) :: zEw ! specific enthalpy of exchanged water (J/kg)
	78	REAL(wp) :: zdE ! specific enthalpy difference (J/kg)
	79	REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean
	80
	81	REAL(wp), DIMENSION(jpij) :: zqprec ! energy of fallen snow (J.m-3)
[9922]	82	REAL(wp), DIMENSION(jpij) :: zq_top ! heat for surface ablation (J.m-2)
	83	REAL(wp), DIMENSION(jpij) :: zq_bot ! heat for bottom ablation (J.m-2)
[8586]	84	REAL(wp), DIMENSION(jpij) :: zq_rema ! remaining heat at the end of the routine (J.m-2)
	85	REAL(wp), DIMENSION(jpij) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2)
	86	REAL(wp), DIMENSION(jpij) :: zevap_rema ! remaining mass flux from sublimation (kg.m-2)
	87
	88	REAL(wp), DIMENSION(jpij) :: zdh_s_mel ! snow melt
	89	REAL(wp), DIMENSION(jpij) :: zdh_s_pre ! snow precipitation
	90	REAL(wp), DIMENSION(jpij) :: zdh_s_sub ! snow sublimation
	91
[9271]	92	REAL(wp), DIMENSION(jpij,nlay_s) :: zh_s ! snw layer thickness
	93	REAL(wp), DIMENSION(jpij,nlay_i) :: zh_i ! ice layer thickness
[8586]	94	REAL(wp), DIMENSION(jpij,nlay_i) :: zdeltah
	95	INTEGER , DIMENSION(jpij,nlay_i) :: icount ! number of layers vanished by melting
	96
	97	REAL(wp), DIMENSION(jpij) :: zsnw ! distribution of snow after wind blowing
	98
	99	REAL(wp) :: zswitch_sal
	100
	101	INTEGER :: num_iter_max ! Heat conservation
	102	!!------------------------------------------------------------------
	103
	104	! Discriminate between time varying salinity and constant
	105	SELECT CASE( nn_icesal ) ! varying salinity or not
	106	CASE( 1 , 3 ) ; zswitch_sal = 0._wp ! prescribed salinity profile
	107	CASE( 2 ) ; zswitch_sal = 1._wp ! varying salinity profile
	108	END SELECT
	109
	110	! initialize layer thicknesses and enthalpies
	111	h_i_old (1:npti,0:nlay_i+1) = 0._wp
	112	eh_i_old(1:npti,0:nlay_i+1) = 0._wp
	113	DO jk = 1, nlay_i
	114	DO ji = 1, npti
	115	h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i
	116	eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk)
	117	END DO
	118	END DO
	119	!
[9604]	120	! ! ============================================== !
	121	! ! Available heat for surface and bottom ablation !
	122	! ! ============================================== !
[8813]	123	!
[10534]	124	IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN
[8813]	125	!
	126	DO ji = 1, npti
[9922]	127	zq_top(ji) = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
[8813]	128	END DO
	129	!
[10534]	130	ELSE
[8813]	131	!
	132	DO ji = 1, npti
[9916]	133	zdum = qns_ice_1d(ji) + qsr_ice_1d(ji) - qtr_ice_top_1d(ji) - qcn_ice_top_1d(ji)
[9274]	134	qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
[9922]	135	zq_top(ji) = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
[8813]	136	END DO
	137	!
[10534]	138	ENDIF
[8813]	139	!
[8586]	140	DO ji = 1, npti
[13662]	141	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)
[9922]	142	zq_bot(ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
[8586]	143	END DO
	144
[9274]	145	! Ice and snow layer thicknesses
	146	!-------------------------------
[9271]	147	DO jk = 1, nlay_i
	148	DO ji = 1, npti
	149	zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i
	150	END DO
	151	END DO
	152
	153	DO jk = 1, nlay_s
	154	DO ji = 1, npti
	155	zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
	156	END DO
	157	END DO
	158
[9274]	159	! ! ============ !
	160	! ! Snow !
	161	! ! ============ !
	162	!
	163	! Internal melting
	164	! ----------------
	165	! IF snow temperature is above freezing point, THEN snow melts (should not happen but sometimes it does)
[9271]	166	DO jk = 1, nlay_s
[8586]	167	DO ji = 1, npti
[9274]	168	IF( t_s_1d(ji,jk) > rt0 ) THEN
	169	hfx_res_1d (ji) = hfx_res_1d (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice ! heat flux to the ocean [W.m-2], < 0
[9935]	170	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice ! mass flux
[9271]	171	! updates
[9274]	172	dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk)
	173	h_s_1d (ji) = h_s_1d(ji) - zh_s(ji,jk)
	174	zh_s (ji,jk) = 0._wp
	175	e_s_1d (ji,jk) = 0._wp
	176	t_s_1d (ji,jk) = rt0
[9271]	177	END IF
[8586]	178	END DO
[9274]	179	END DO
[8586]	180
[9274]	181	! Snow precipitation
	182	!-------------------
[13662]	183	CALL ice_var_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing
[8586]	184
	185	zdeltah(1:npti,:) = 0._wp
	186	DO ji = 1, npti
[9274]	187	IF( sprecip_1d(ji) > 0._wp ) THEN
	188	!
	189	! --- precipitation ---
[9935]	190	zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhos / at_i_1d(ji) ! thickness change
[9274]	191	zqprec (ji) = - qprec_ice_1d(ji) ! enthalpy of the precip (>0, J.m-3)
	192	!
	193	hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice ! heat flux from snow precip (>0, W.m-2)
[9935]	194	wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice ! mass flux, <0
[9274]	195
	196	! --- melt of falling snow ---
	197	rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
[9922]	198	zdeltah (ji,1) = - rswitch * zq_top(ji) / MAX( zqprec(ji) , epsi20 ) ! thickness change
	199	zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting
[9274]	200	hfx_snw_1d (ji) = hfx_snw_1d (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice ! heat used to melt snow (W.m-2, >0)
[9935]	201	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice ! snow melting only = water into the ocean (then without snow precip), >0
[9274]	202
	203	! updates available heat + precipitations after melting
	204	dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1)
[9922]	205	zq_top (ji) = MAX( 0._wp , zq_top (ji) + zdeltah(ji,1) * zqprec(ji) )
[9274]	206	zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
	207
	208	! update thickness
	209	h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) )
	210	!
	211	ELSE
	212	!
	213	zdh_s_pre(ji) = 0._wp
	214	zqprec (ji) = 0._wp
	215	!
	216	ENDIF
[8586]	217	END DO
	218
[9271]	219	! recalculate snow layers
	220	DO jk = 1, nlay_s
	221	DO ji = 1, npti
	222	zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
	223	END DO
	224	END DO
	225
[9274]	226	! Snow melting
	227	! ------------
[9922]	228	! If heat still available (zq_top > 0), then melt more snow
[8586]	229	zdeltah(1:npti,:) = 0._wp
	230	zdh_s_mel(1:npti) = 0._wp
	231	DO jk = 1, nlay_s
	232	DO ji = 1, npti
[9922]	233	IF( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN
[9274]	234	!
	235	rswitch = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) )
[9922]	236	zdeltah (ji,jk) = - rswitch * zq_top(ji) / MAX( e_s_1d(ji,jk), epsi20 ) ! thickness change
	237	zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) ) ! bound melting
[9274]	238	zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk)
	239
	240	hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_rdtice ! heat used to melt snow(W.m-2, >0)
[9935]	241	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! snow melting only = water into the ocean (then without snow precip)
[9274]	242
	243	! updates available heat + thickness
	244	dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,jk)
[9922]	245	zq_top (ji) = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
[9274]	246	h_s_1d (ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) )
	247	zh_s (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) )
	248	!
	249	ENDIF
[8586]	250	END DO
	251	END DO
	252
[9274]	253	! Snow sublimation
	254	!-----------------
[8586]	255	! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
[8885]	256	! comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
[8586]	257	zdeltah(1:npti,:) = 0._wp
	258	DO ji = 1, npti
[9274]	259	IF( evap_ice_1d(ji) > 0._wp ) THEN
	260	!
[9935]	261	zdh_s_sub (ji) = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos * rdt_ice )
	262	zevap_rema(ji) = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhos ! remaining evap in kg.m-2 (used for ice melting later on)
[9274]	263	zdeltah (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
	264
	265	hfx_sub_1d (ji) = hfx_sub_1d(ji) + & ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow)
	266	& ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) ) &
	267	& * a_i_1d(ji) * r1_rdtice
[9935]	268	wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice ! Mass flux by sublimation
[9274]	269
	270	! new snow thickness
	271	h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) )
	272	! update precipitations after sublimation and correct sublimation
	273	zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
	274	zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
	275	!
	276	ELSE
	277	!
	278	zdh_s_sub (ji) = 0._wp
	279	zevap_rema(ji) = 0._wp
	280	!
	281	ENDIF
[8586]	282	END DO
	283
	284	! --- Update snow diags --- !
	285	DO ji = 1, npti
	286	dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
	287	END DO
	288
[9274]	289	! Update temperature, energy
	290	!---------------------------
[8586]	291	! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
	292	DO jk = 1, nlay_s
	293	DO ji = 1,npti
	294	rswitch = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) )
[9274]	295	e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) * &
	296	& ( ( zdh_s_pre(ji) ) * zqprec(ji) + &
[9935]	297	& ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) )
[8586]	298	END DO
	299	END DO
[9274]	300
	301	! ! ============ !
	302	! ! Ice !
	303	! ! ============ !
[8586]	304
[9274]	305	! Surface ice melting
	306	!--------------------
[8586]	307	zdeltah(1:npti,:) = 0._wp ! important
	308	DO jk = 1, nlay_i
	309	DO ji = 1, npti
[9935]	310	ztmelts = - rTmlt * sz_i_1d(ji,jk) ! Melting point of layer k [C]
[8586]	311
[9274]	312	IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting
[8586]	313
[9935]	314	zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of layer k [J/kg, <0]
[9274]	315	zdE = 0._wp ! Specific enthalpy difference (J/kg, <0)
[8586]	316	! set up at 0 since no energy is needed to melt water...(it is already melted)
	317	zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing
	318	! this should normally not happen, but sometimes, heat diffusion leads to this
[9935]	319	zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0
[8586]	320
[9750]	321	dh_i_itm(ji) = dh_i_itm(ji) + zdeltah(ji,jk) ! Cumulate internal melting
[8586]	322
[9935]	323	zfmdt = - rhoi * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0]
[8586]	324
[9274]	325	hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice ! Heat flux to the ocean [W.m-2], <0
	326	! ice enthalpy zEi is "sent" to the ocean
[9935]	327	sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux
[9274]	328	! using s_i_1d and not sz_i_1d(jk) is ok
[9935]	329	wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux
[8586]	330
[9274]	331	ELSE !-- Surface melting
[8586]	332
[9935]	333	zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of layer k [J/kg, <0]
[8586]	334	zEw = rcp * ztmelts ! Specific enthalpy of resulting meltwater [J/kg, <0]
	335	zdE = zEi - zEw ! Specific enthalpy difference < 0
	336
[9922]	337	zfmdt = - zq_top(ji) / zdE ! Mass flux to the ocean [kg/m2, >0]
[8586]	338
[9935]	339	zdeltah(ji,jk) = - zfmdt * r1_rhoi ! Melt of layer jk [m, <0]
[8586]	340
	341	zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0]
	342
[9935]	343	zq_top(ji) = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat
[8586]	344
[9750]	345	dh_i_sum(ji) = dh_i_sum(ji) + zdeltah(ji,jk) ! Cumulate surface melt
[8586]	346
[9935]	347	zfmdt = - rhoi * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0]
[8586]	348
	349	zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0]
	350
[9935]	351	sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux >0
[9274]	352	! using s_i_1d and not sz_i_1d(jk) is ok)
	353	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux [W.m-2], < 0
	354	hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Heat flux used in this process [W.m-2], > 0
	355	!
[9935]	356	wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux
[8586]	357
	358	END IF
[9274]	359
	360	! Ice sublimation
	361	! ---------------
[9935]	362	zdum = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi )
[9274]	363	zdeltah (ji,jk) = zdeltah (ji,jk) + zdum
	364	dh_i_sub(ji) = dh_i_sub(ji) + zdum
	365
[9935]	366	sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice ! Salt flux >0
[9274]	367	! clem: flux is sent to the ocean for simplicity
	368	! but salt should remain in the ice except
	369	! if all ice is melted. => must be corrected
	370	hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice ! Heat flux [W.m-2], < 0
[8586]	371
[9935]	372	wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_rdtice ! Mass flux > 0
[9274]	373
[8586]	374	! update remaining mass flux
[9935]	375	zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi
[8586]	376
	377	! record which layers have disappeared (for bottom melting)
	378	! => icount=0 : no layer has vanished
	379	! => icount=5 : 5 layers have vanished
	380	rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )
	381	icount(ji,jk) = NINT( rswitch )
	382	zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
	383
	384	! update heat content (J.m-2) and layer thickness
	385	eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
	386	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
	387	END DO
	388	END DO
[9274]	389
[8586]	390	! update ice thickness
	391	DO ji = 1, npti
[9750]	392	h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_sum(ji) + dh_i_itm(ji) + dh_i_sub(ji) )
[8586]	393	END DO
	394
	395	! remaining "potential" evap is sent to ocean
	396	DO ji = 1, npti
	397	wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_rdtice ! <=0 (net evap for the ocean in kg.m-2.s-1)
	398	END DO
	399
	400
[9274]	401	! Ice Basal growth
[8586]	402	!------------------
	403	! Basal growth is driven by heat imbalance at the ice-ocean interface,
[9916]	404	! between the inner conductive flux (qcn_ice_bot), from the open water heat flux
[9913]	405	! (fhld) and the sensible ice-ocean flux (qsb_ice_bot).
[9916]	406	! qcn_ice_bot is positive downwards. qsb_ice_bot and fhld are positive to the ice
[8586]	407
	408	! If salinity varies in time, an iterative procedure is required, because
	409	! the involved quantities are inter-dependent.
[9750]	410	! Basal growth (dh_i_bog) depends upon new ice specific enthalpy (zEi),
	411	! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bog)
[8586]	412	! -> need for an iterative procedure, which converges quickly
	413
	414	num_iter_max = 1
	415	IF( nn_icesal == 2 ) num_iter_max = 5 ! salinity varying in time
	416
	417	DO ji = 1, npti
	418	IF( zf_tt(ji) < 0._wp ) THEN
[9274]	419	DO iter = 1, num_iter_max ! iterations
[8586]	420
	421	! New bottom ice salinity (Cox & Weeks, JGR88 )
	422	!--- zswi1 if dh/dt < 2.0e-8
	423	!--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
	424	!--- zswi2 if dh/dt > 3.6e-7
[9750]	425	zgrr = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_rdtice , epsi10 ) )
[9274]	426	zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
	427	zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
	428	zswi1 = 1. - zswi2 * zswi12
	429	zfracs = MIN( zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) &
	430	& + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 )
[8586]	431
[9274]	432	s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji) ! New ice salinity
[8586]	433
[9935]	434	ztmelts = - rTmlt * s_i_new(ji) ! New ice melting point (C)
[9274]	435
	436	zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
[8586]	437
[9935]	438	zEi = rcpi * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0)
	439	& - rLfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp * ztmelts
[8586]	440
[9274]	441	zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)
[8586]	442
[9274]	443	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
[8586]	444
[9935]	445	dh_i_bog(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) )
[8586]	446
	447	END DO
	448	! Contribution to Energy and Salt Fluxes
[9935]	449	zfmdt = - rhoi * dh_i_bog(ji) ! Mass flux x time step (kg/m2, < 0)
[8586]	450
[9274]	451	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux to the ocean [W.m-2], >0
	452	hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Heat flux used in this process [W.m-2], <0
[8586]	453
[9935]	454	sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_rdtice ! Salt flux, <0
[8586]	455
[9935]	456	wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_rdtice ! Mass flux, <0
[8586]	457
	458	! update heat content (J.m-2) and layer thickness
[9935]	459	eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi)
[9750]	460	h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bog(ji)
[8586]	461
	462	ENDIF
	463
	464	END DO
	465
[9274]	466	! Ice Basal melt
	467	!---------------
[8586]	468	zdeltah(1:npti,:) = 0._wp ! important
	469	DO jk = nlay_i, 1, -1
	470	DO ji = 1, npti
	471	IF( zf_tt(ji) > 0._wp .AND. jk > icount(ji,jk) ) THEN ! do not calculate where layer has already disappeared by surface melting
	472
[9935]	473	ztmelts = - rTmlt * sz_i_1d(ji,jk) ! Melting point of layer jk (C)
[8586]	474
[9274]	475	IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting
[8586]	476
[9935]	477	zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0)
[8586]	478	zdE = 0._wp ! Specific enthalpy difference (J/kg, <0)
[9274]	479	! set up at 0 since no energy is needed to melt water...(it is already melted)
[8586]	480	zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing
	481	! this should normally not happen, but sometimes, heat diffusion leads to this
	482
[9750]	483	dh_i_itm (ji) = dh_i_itm(ji) + zdeltah(ji,jk)
[8586]	484
[9935]	485	zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0
[8586]	486
[9274]	487	hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice ! Heat flux to the ocean [W.m-2], <0
	488	! ice enthalpy zEi is "sent" to the ocean
[9935]	489	sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux
[9274]	490	! using s_i_1d and not sz_i_1d(jk) is ok
[9935]	491	wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux
[8586]	492
	493	! update heat content (J.m-2) and layer thickness
	494	eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
	495	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
	496
[9274]	497	ELSE !-- Basal melting
[8586]	498
[9935]	499	zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0)
[9274]	500	zEw = rcp * ztmelts ! Specific enthalpy of meltwater (J/kg, <0)
	501	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
[8586]	502
[9922]	503	zfmdt = - zq_bot(ji) / zdE ! Mass flux x time step (kg/m2, >0)
[8586]	504
[9935]	505	zdeltah(ji,jk) = - zfmdt * r1_rhoi ! Gross thickness change
[8586]	506
[9274]	507	zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
[8586]	508
[9935]	509	zq_bot(ji) = MAX( 0._wp , zq_bot(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat. MAX is necessary for roundup errors
[8586]	510
[9922]	511	dh_i_bom(ji) = dh_i_bom(ji) + zdeltah(ji,jk) ! Update basal melt
[8586]	512
[9935]	513	zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0
[8586]	514
[9274]	515	zQm = zfmdt * zEw ! Heat exchanged with ocean
[8586]	516
[9274]	517	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux to the ocean [W.m-2], <0
	518	hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Heat used in this process [W.m-2], >0
[8586]	519
[9935]	520	sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux
[9274]	521	! using s_i_1d and not sz_i_1d(jk) is ok
[8586]	522
[9935]	523	wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux
[8586]	524
	525	! update heat content (J.m-2) and layer thickness
	526	eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
	527	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
	528	ENDIF
	529
	530	ENDIF
	531	END DO
	532	END DO
	533
	534	! Update temperature, energy
[9274]	535	! --------------------------
[8586]	536	DO ji = 1, npti
[9750]	537	h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) )
[8586]	538	END DO
	539
[9274]	540	! If heat still available then melt more snow
[8586]	541	!-------------------------------------------
	542	zdeltah(1:npti,:) = 0._wp ! important
	543	DO ji = 1, npti
[9922]	544	zq_rema (ji) = zq_top(ji) + zq_bot(ji)
[9274]	545	rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ! =1 if snow
	546	rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) )
	547	zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 )
	548	zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting
	549	dh_s_tot(ji) = dh_s_tot(ji) + zdeltah(ji,1)
	550	h_s_1d (ji) = h_s_1d (ji) + zdeltah(ji,1)
[8586]	551
[9274]	552	zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1) ! update available heat (J.m-2)
	553	hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! Heat used to melt snow, W.m-2 (>0)
[9935]	554	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice ! Mass flux
[8637]	555	dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,1)
[8586]	556	!
	557	! Remaining heat flux (W.m-2) is sent to the ocean heat budget
[13662]	558	!!!hfx_res_1d(ji) = hfx_res_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
[8586]	559
	560	IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
	561	END DO
	562
	563	!
[9274]	564	! Snow-Ice formation
	565	! ------------------
[8586]	566	! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
	567	! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
[9935]	568	z1_rho = 1._wp / ( rhos+rau0-rhoi )
[8586]	569	DO ji = 1, npti
	570	!
[9935]	571	dh_snowice(ji) = MAX( 0._wp , ( rhos * h_s_1d(ji) + (rhoi-rau0) * h_i_1d(ji) ) * z1_rho )
[8586]	572
	573	h_i_1d(ji) = h_i_1d(ji) + dh_snowice(ji)
	574	h_s_1d(ji) = h_s_1d(ji) - dh_snowice(ji)
	575
	576	! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
[9935]	577	zfmdt = ( rhos - rhoi ) * dh_snowice(ji) ! <0
[8586]	578	zEw = rcp * sst_1d(ji)
	579	zQm = zfmdt * zEw
	580
[9274]	581	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux
[8586]	582
[9274]	583	sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice ! Salt flux
[8586]	584
	585	! Case constant salinity in time: virtual salt flux to keep salinity constant
[8637]	586	IF( nn_icesal /= 2 ) THEN
[8586]	587	sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt * r1_rdtice & ! put back sss_m into the ocean
[9935]	588	& - s_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_rdtice ! and get rn_icesal from the ocean
[8586]	589	ENDIF
	590
[9274]	591	! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume
[9935]	592	wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_rdtice
	593	wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_rdtice
[8586]	594
	595	! update heat content (J.m-2) and layer thickness
	596	eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw
	597	h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
	598
	599	END DO
	600
	601	!
	602	! Update temperature, energy
[9274]	603	! --------------------------
[8586]	604	DO ji = 1, npti
[9274]	605	rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) )
	606	t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0
[8586]	607	END DO
	608
	609	DO jk = 1, nlay_s
	610	DO ji = 1,npti
[10786]	611	! where there is no ice or no snow
	612	rswitch = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ) * ( 1._wp - MAX( 0._wp, SIGN(1._wp, - h_i_1d(ji) ) ) )
	613	! mass & energy loss to the ocean
	614	hfx_res_1d(ji) = hfx_res_1d(ji) + ( 1._wp - rswitch ) * &
	615	& ( e_s_1d(ji,jk) * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_rdtice ) ! heat flux to the ocean [W.m-2], < 0
	616	wfx_res_1d(ji) = wfx_res_1d(ji) + ( 1._wp - rswitch ) * &
	617	& ( rhos * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_rdtice ) ! mass flux
	618	! update energy (mass is updated in the next loop)
[8586]	619	e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk)
	620	! recalculate t_s_1d from e_s_1d
[9935]	621	t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi )
[8586]	622	END DO
	623	END DO
	624
[10786]	625	! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 ---
	626	WHERE( h_i_1d(1:npti) == 0._wp )
	627	a_i_1d(1:npti) = 0._wp
	628	h_s_1d(1:npti) = 0._wp
	629	END WHERE
[8586]	630	!
	631	END SUBROUTINE ice_thd_dh
	632
	633	#else
	634	!!----------------------------------------------------------------------
[9570]	635	!! Default option NO SI3 sea-ice model
[8586]	636	!!----------------------------------------------------------------------
	637	#endif
	638
	639	!!======================================================================
	640	END MODULE icethd_dh

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

Download in other formats: