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.

limthd_dh.F90 in branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 4635

Last change on this file since 4635 was 4634, checked in by clem, 10 years ago
major changes in heat budget
Property svn:keywords set to `Id`
File size: 37.0 KB

Rev	Line
[825]	1	MODULE limthd_dh
[1572]	2	!!======================================================================
	3	!! * MODULE limthd_dh *
	4	!! LIM-3 : thermodynamic growth and decay of the ice
	5	!!======================================================================
	6	!! History : LIM ! 2003-05 (M. Vancoppenolle) Original code in 1D
	7	!! ! 2005-06 (M. Vancoppenolle) 3D version
[4634]	8	!! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw & wfx_ice
[3625]	9	!! 3.4 ! 2011-02 (G. Madec) dynamical allocation
	10	!! 3.5 ! 2012-10 (G. Madec & co) salt flux + bug fixes
[1572]	11	!!----------------------------------------------------------------------
[825]	12	#if defined key_lim3
[834]	13	!!----------------------------------------------------------------------
	14	!! 'key_lim3' LIM3 sea-ice model
	15	!!----------------------------------------------------------------------
[3625]	16	!! lim_thd_dh : vertical accr./abl. and lateral ablation of sea ice
[825]	17	!!----------------------------------------------------------------------
[3625]	18	USE par_oce ! ocean parameters
	19	USE phycst ! physical constants (OCE directory)
	20	USE sbc_oce ! Surface boundary condition: ocean fields
	21	USE ice ! LIM variables
	22	USE par_ice ! LIM parameters
	23	USE thd_ice ! LIM thermodynamics
	24	USE in_out_manager ! I/O manager
	25	USE lib_mpp ! MPP library
	26	USE wrk_nemo ! work arrays
	27	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[4634]	28	USE cpl_oasis3, ONLY : lk_cpl
	29
[825]	30	IMPLICIT NONE
	31	PRIVATE
	32
[1572]	33	PUBLIC lim_thd_dh ! called by lim_thd
[825]	34
[4332]	35	REAL(wp) :: epsi20 = 1.e-20 ! constant values
	36	REAL(wp) :: epsi10 = 1.e-10 !
[825]	37
	38	!!----------------------------------------------------------------------
[4045]	39	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
[1156]	40	!! $Id$
[2715]	41	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	42	!!----------------------------------------------------------------------
	43	CONTAINS
	44
[2715]	45	SUBROUTINE lim_thd_dh( kideb, kiut, jl )
[921]	46	!!------------------------------------------------------------------
	47	!! * ROUTINE lim_thd_dh *
	48	!!
[1572]	49	!! ** Purpose : determines variations of ice and snow thicknesses.
[921]	50	!!
[1572]	51	!! ** Method : Ice/Snow surface melting arises from imbalance in surface fluxes
	52	!! Bottom accretion/ablation arises from flux budget
	53	!! Snow thickness can increase by precipitation and decrease by sublimation
	54	!! If snow load excesses Archmiede limit, snow-ice is formed by
	55	!! the flooding of sea-water in the snow
[921]	56	!!
[1572]	57	!! 1) Compute available flux of heat for surface ablation
	58	!! 2) Compute snow and sea ice enthalpies
	59	!! 3) Surface ablation and sublimation
	60	!! 4) Bottom accretion/ablation
	61	!! 5) Case of Total ablation
	62	!! 6) Snow ice formation
[921]	63	!!
[1572]	64	!! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
	65	!! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646
	66	!! Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let.
	67	!! Vancoppenolle et al.,2009, Ocean Modelling
[921]	68	!!------------------------------------------------------------------
[1572]	69	INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied
	70	INTEGER , INTENT(in) :: jl ! Thickness cateogry number
	71	!!
	72	INTEGER :: ji , jk ! dummy loop indices
[4045]	73	INTEGER :: ii, ij ! 2D corresponding indices to ji
[1572]	74	INTEGER :: i_ice_switch ! ice thickness above a certain treshold or not
	75	INTEGER :: iter
[825]	76
[4634]	77	REAL(wp) :: ztmelts ! local scalar
	78	REAL(wp) :: zdh, zfdum !
[1572]	79	REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment
	80	REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads
[4634]	81	REAL(wp) :: zs_snic ! snow-ice salinity
[1572]	82	REAL(wp) :: zswi1 ! switch for computation of bottom salinity
	83	REAL(wp) :: zswi12 ! switch for computation of bottom salinity
	84	REAL(wp) :: zswi2 ! switch for computation of bottom salinity
	85	REAL(wp) :: zgrr ! bottom growth rate
[4634]	86	REAL(wp) :: zt_i_new ! bottom formation temperature
	87
	88	REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
	89	REAL(wp) :: zEi ! specific enthalpy of sea ice (J/kg)
	90	REAL(wp) :: zEw ! specific enthalpy of exchanged water (J/kg)
	91	REAL(wp) :: zdE ! specific enthalpy difference (J/kg)
	92	REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean
	93	REAL(wp) :: zsstK ! SST in Kelvin
	94
[3294]	95	REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness
[4634]	96	REAL(wp), POINTER, DIMENSION(:) :: zqprec ! energy of fallen snow (J.m-3)
	97	REAL(wp), POINTER, DIMENSION(:) :: zq_su ! heat for surface ablation (J.m-2)
	98	REAL(wp), POINTER, DIMENSION(:) :: zq_bo ! heat for bottom ablation (J.m-2)
	99	REAL(wp), POINTER, DIMENSION(:) :: zq_1cat ! corrected heat in case 1-cat and hmelt>15cm (J.m-2)
	100	REAL(wp), POINTER, DIMENSION(:) :: zq_rema ! remaining heat at the end of the routine (J.m-2)
	101	REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2)
	102	INTEGER , POINTER, DIMENSION(:) :: icount ! number of layers vanished by melting
[3294]	103
[3625]	104	REAL(wp), POINTER, DIMENSION(:) :: zdh_s_mel ! snow melt
	105	REAL(wp), POINTER, DIMENSION(:) :: zdh_s_pre ! snow precipitation
	106	REAL(wp), POINTER, DIMENSION(:) :: zdh_s_sub ! snow sublimation
[3294]	107
	108	REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah
[4634]	109	REAL(wp), POINTER, DIMENSION(:,:) :: zh_i ! ice layer thickness
[3294]	110
[4634]	111	REAL(wp), POINTER, DIMENSION(:) :: zqh_i ! total ice heat content (J.m-2)
	112	REAL(wp), POINTER, DIMENSION(:) :: zqh_s ! total snow heat content (J.m-2)
	113	REAL(wp), POINTER, DIMENSION(:) :: zq_s ! total snow enthalpy (J.m-3)
[3294]	114
[4045]	115	! mass and salt flux (clem)
[4634]	116	REAL(wp) :: zdvres, zswitch_sal
[4045]	117
[3294]	118	! Heat conservation
[4634]	119	INTEGER :: num_iter_max
	120	REAL(wp) :: zinda, zindq, zindh
	121	REAL(wp), POINTER, DIMENSION(:) :: zintermelt ! debug
	122
[1572]	123	!!------------------------------------------------------------------
[825]	124
[4634]	125	! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values)
	126	SELECT CASE( num_sal ) ! varying salinity or not
	127	CASE( 1, 3, 4 ) ; zswitch_sal = 0 ! prescribed salinity profile
	128	CASE( 2 ) ; zswitch_sal = 1 ! varying salinity profile
	129	END SELECT
[825]	130
[4634]	131	CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema )
	132	CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
	133	CALL wrk_alloc( jpij, zintermelt )
	134	CALL wrk_alloc( jpij, jkmax, zdeltah, zh_i )
	135	CALL wrk_alloc( jpij, icount )
[4045]	136
[4634]	137	dh_i_surf (:) = 0._wp ; dh_i_bott (:) = 0._wp ; dh_snowice(:) = 0._wp
	138	dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp
	139
	140	zqprec (:) = 0._wp ; zq_su (:) = 0._wp ; zq_bo (:) = 0._wp ; zf_tt (:) = 0._wp
	141	zq_1cat(:) = 0._wp ; zq_rema(:) = 0._wp
[2715]	142
[4634]	143	zh_s (:) = 0._wp
	144	zdh_s_pre(:) = 0._wp
	145	zdh_s_mel(:) = 0._wp
	146	zdh_s_sub(:) = 0._wp
	147	zqh_s (:) = 0._wp
	148	zqh_i (:) = 0._wp
[4045]	149
[4634]	150	zh_i (:,:) = 0._wp
	151	zdeltah (:,:) = 0._wp
	152	zintermelt(:) = 0._wp
	153	icount (:) = 0
	154
	155	! debug
	156	dq_i(:) = 0._wp
	157	dq_s(:) = 0._wp
	158
	159	! initialize layer thicknesses and enthalpies
	160	h_i_old (:,0:nlay_i+1) = 0._wp
	161	qh_i_old(:,0:nlay_i+1) = 0._wp
	162	DO jk = 1, nlay_i
	163	DO ji = kideb, kiut
	164	h_i_old (ji,jk) = ht_i_b(ji) / REAL( nlay_i )
	165	qh_i_old(ji,jk) = q_i_b(ji,jk) * h_i_old(ji,jk)
	166	ENDDO
	167	ENDDO
[921]	168	!
	169	!------------------------------------------------------------------------------!
[4634]	170	! 1) Calculate available heat for surface and bottom ablation !
[921]	171	!------------------------------------------------------------------------------!
	172	!
[2715]	173	DO ji = kideb, kiut
[4634]	174	zinda = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )
	175	ztmelts = zinda * rtt + ( 1._wp - zinda ) * rtt
	176
	177	zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)
	178	zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)
	179
	180	zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_b(ji) - ztmelts ) )
	181	zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
[825]	182	END DO ! ji
	183
[921]	184	!
	185	!------------------------------------------------------------------------------!
[4634]	186	! If snow temperature is above freezing point, then snow melts
	187	! (should not happen but sometimes it does)
[921]	188	!------------------------------------------------------------------------------!
[4634]	189	DO ji = kideb, kiut
	190	IF( t_s_b(ji,1) > rtt ) THEN !!! Internal melting
	191	! Contribution to heat flux to the ocean [W.m-2], < 0
	192	hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_b(ji,1) * ht_s_b(ji) * a_i_b(ji) * r1_rdtice
	193	! Contribution to mass flux
	194	wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * ht_s_b(ji) * a_i_b(ji) * r1_rdtice
	195	! updates
	196	ht_s_b(ji) = 0._wp
	197	q_s_b (ji,1) = 0._wp
	198	t_s_b (ji,1) = rtt
	199	END IF
	200	END DO
	201
	202	!------------------------------------------------------------!
	203	! 2) Computing layer thicknesses and enthalpies. !
	204	!------------------------------------------------------------!
[921]	205	!
[4634]	206	DO ji = kideb, kiut
[4045]	207	zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
[825]	208	END DO
[2715]	209	!
[825]	210	DO jk = 1, nlay_s
[2715]	211	DO ji = kideb, kiut
[4634]	212	zqh_s(ji) = zqh_s(ji) + q_s_b(ji,jk) * zh_s(ji)
[825]	213	END DO
	214	END DO
[2715]	215	!
[825]	216	DO jk = 1, nlay_i
[2715]	217	DO ji = kideb, kiut
[4634]	218	zh_i(ji,jk) = ht_i_b(ji) / REAL( nlay_i )
	219	zqh_i(ji) = zqh_i(ji) + q_i_b(ji,jk) * zh_i(ji,jk)
[825]	220	END DO
	221	END DO
[921]	222	!
	223	!------------------------------------------------------------------------------\|
	224	! 3) Surface ablation and sublimation \|
	225	!------------------------------------------------------------------------------\|
	226	!
[834]	227	!-------------------------
	228	! 3.1 Snow precips / melt
	229	!-------------------------
[825]	230	! Snow accumulation in one thermodynamic time step
	231	! snowfall is partitionned between leads and ice
	232	! if snow fall was uniform, a fraction (1-at_i) would fall into leads
	233	! but because of the winds, more snow falls on leads than on sea ice
	234	! and a greater fraction (1-at_i)^beta of the total mass of snow
[834]	235	! (beta < 1) falls in leads.
[825]	236	! In reality, beta depends on wind speed,
	237	! and should decrease with increasing wind speed but here, it is
[834]	238	! considered as a constant. an average value is 0.66
[825]	239	! Martin Vancoppenolle, December 2006
	240
	241	DO ji = kideb, kiut
[4634]	242	!-----------
	243	! Snow fall
	244	!-----------
	245	! thickness change
	246	zcoeff = ( 1._wp - ( 1._wp - at_i_b(ji) )**betas ) / at_i_b(ji)
[825]	247	zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
[4634]	248	! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K)
	249	zqprec (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )
	250	IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
	251	! heat flux from snow precip (>0, W.m-2)
	252	hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice
	253	! update thickness
	254	ht_s_b (ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_pre(ji) )
[825]	255
[4634]	256	!---------------------
	257	! Melt of falling snow
	258	!---------------------
	259	! thickness change
	260	zindq = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) )
	261	zdh_s_mel (ji) = - zindq * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
	262	zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting
	263	! Heat flux associated with snow melt of falling snow (W.m-2, <0)
	264	hfx_snw_1d(ji) = hfx_snw_1d(ji) + zdh_s_mel(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice
	265	! heat used to melt snow (W.m-2, >0)
	266	hfx_tot_1d(ji) = hfx_tot_1d(ji) - zdh_s_mel(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice
	267	! snow melting only = water into the ocean (then without snow precip)
	268	wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdh_s_mel(ji) * r1_rdtice
	269
	270	! updates available heat + thickness
	271	zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) )
	272	ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_mel(ji) )
	273	zh_s (ji) = ht_s_b(ji) / REAL( nlay_s )
	274
	275	! clem debug: variation of enthalpy (J.m-2)
	276	dq_s(ji) = dq_s(ji) + ( zdh_s_pre(ji) + zdh_s_mel(ji) ) * zqprec(ji)
	277
[825]	278	END DO
	279
[4634]	280	! If heat still available, then melt more snow
	281	zdeltah(:,:) = 0._wp ! important
[825]	282	DO jk = 1, nlay_s
	283	DO ji = kideb, kiut
[4634]	284	! thickness change
	285	zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_b(ji) ) )
	286	zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_b(ji,jk) + epsi20 ) )
	287	zdeltah (ji,jk) = - zindh * zindq * zq_su(ji) / MAX( q_s_b(ji,jk), epsi20 )
	288	zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting
	289	zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk)
	290	! heat flux associated with snow melt(W.m-2, <0)
	291	hfx_snw_1d(ji) = hfx_snw_1d(ji) + zdeltah(ji,jk) * a_i_b(ji) * q_s_b(ji,jk) * r1_rdtice
	292	! heat used to melt snow(W.m-2, >0)
	293	hfx_tot_1d(ji) = hfx_tot_1d(ji) - zdeltah(ji,jk) * a_i_b(ji) * q_s_b(ji,jk) * r1_rdtice
	294	! snow melting only = water into the ocean (then without snow precip)
	295	wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
[825]	296
[4634]	297	! updates available heat + thickness
	298	zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_b(ji,jk) )
	299	ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdeltah(ji,jk) )
[825]	300
[4634]	301	! clem debug: variation of enthalpy (J.m-2)
	302	dq_s(ji) = dq_s(ji) + zdeltah(ji,jk) * q_s_b(ji,jk)
[1572]	303	END DO
	304	END DO
[825]	305
[4634]	306	!----------------------
	307	! 3.2 Snow sublimation
	308	!----------------------
	309	! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
	310	IF( lk_cpl ) THEN
	311	! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice)
	312	zdh_s_sub(:) = 0._wp
	313	ELSE
	314	! forced mode: snow thickness change due to sublimation
[921]	315	DO ji = kideb, kiut
[4634]	316	zdh_s_sub(ji) = MAX( - ht_s_b(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice )
	317	! Heat flux by sublimation [W.m-2], < 0
	318	! sublimate first snow that had fallen, then pre-existing snow
	319	zcoeff = ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) * zqprec(ji) + &
	320	& ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_b(ji,1) ) &
	321	& * a_i_b(ji) * r1_rdtice
	322	hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff ! diag only (to close heat budget)
	323	! heat used for sublimation (>0, W.m-2)
	324	!!? hfx_tot_1d(ji) = hfx_tot_1d(ji) - zcoeff
	325	! Mass flux by sublimation
	326	wfx_sub_1d(ji) = wfx_sub_1d(ji) + rhosn * a_i_b(ji) * zdh_s_sub(ji) * r1_rdtice ! diag only
	327	wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdh_s_sub(ji) * r1_rdtice
	328	! new snow thickness
	329	ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_sub(ji) )
	330	! clem debug: variation of enthalpy (J.m-2)
	331	dq_s(ji) = dq_s(ji) + zdh_s_sub(ji) * q_s_b(ji,1)
[1572]	332	END DO
	333	ENDIF
[825]	334
[4634]	335	! --- Update snow diags --- !
[825]	336	DO ji = kideb, kiut
[4634]	337	dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
	338	zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
	339	END DO ! ji
[825]	340
[4634]	341	!-------------------------------------------
	342	! 3.3 Update temperature, energy
	343	!-------------------------------------------
	344	! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
	345	zq_s(:) = 0._wp
[825]	346	DO jk = 1, nlay_s
	347	DO ji = kideb,kiut
[4634]	348	zindh = MAX( 0._wp , SIGN( 1._wp, - ht_s_b(ji) + epsi20 ) )
	349	q_s_b(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_b(ji), epsi20 ) * &
	350	& ( ( MAX( 0._wp, dh_s_tot(ji) ) ) * zqprec(ji) + &
	351	& ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_b(ji) ) * rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) )
	352	zq_s(ji) = zq_s(ji) + q_s_b(ji,jk)
[825]	353	END DO
	354	END DO
	355
[4634]	356	!--------------------------
	357	! 3.4 Surface ice ablation
	358	!--------------------------
	359	zdeltah(:,:) = 0._wp ! important
	360	DO jk = 1, nlay_i
	361	DO ji = kideb, kiut
	362	zEi = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0]
	363
	364	ztmelts = - tmut * s_i_b(ji,jk) + rtt ! Melting point of layer k [K]
	365
	366	zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0]
	367
	368	zdE = zEi - zEw ! Specific enthalpy difference < 0
	369
	370	zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0]
	371
	372	zdeltah(ji,jk) = - zfmdt / rhoic ! Melt of layer jk [m, <0]
	373
	374	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]
	375
	376	zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat
	377
	378	dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt
	379
	380	zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0]
	381
	382	zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0]
	383
	384	! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
	385	sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
	386
	387	! Contribution to heat flux [W.m-2], < 0
	388	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
	389
	390	! Total heat flux used in this process [W.m-2], < 0
	391	hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
	392
	393	! Contribution to mass flux
	394	wfx_sum_1d(ji) = wfx_sum_1d(ji) + rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
	395
	396	! record which layers have disappeared (for bottom melting)
	397	! => icount=0 : no layer has vanished
	398	! => icount=5 : 5 layers have vanished
	399	zindh = NINT( MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) )
	400	icount(ji) = icount(ji) + zindh
	401	zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
	402
	403	! clem debug: variation of enthalpy (J.m-2)
	404	dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)
	405
	406	! update heat content (J.m-2) and layer thickness
	407	qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
	408	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
[825]	409	END DO
[921]	410	END DO
[4634]	411	! update ice thickness
	412	DO ji = kideb, kiut
	413	ht_i_b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_surf(ji) )
	414	END DO
[825]	415
[921]	416	!
	417	!------------------------------------------------------------------------------!
	418	! 4) Basal growth / melt !
	419	!------------------------------------------------------------------------------!
	420	!
[4634]	421	!------------------
	422	! 4.1 Basal growth
	423	!------------------
	424	! Basal growth is driven by heat imbalance at the ice-ocean interface,
	425	! between the inner conductive flux (fc_bo_i), from the open water heat flux
	426	! (fhldb) and the turbulent ocean flux (fhtur).
	427	! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice
[825]	428
[4634]	429	! If salinity varies in time, an iterative procedure is required, because
	430	! the involved quantities are inter-dependent.
	431	! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi),
	432	! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott)
	433	! -> need for an iterative procedure, which converges quickly
	434
	435	IF ( num_sal == 2 ) THEN
	436	num_iter_max = 5
	437	ELSE
	438	num_iter_max = 1
[1572]	439	ENDIF
[825]	440
[4634]	441	!clem debug. Just to be sure that enthalpy at nlay_i+1 is null
	442	DO ji = kideb, kiut
	443	q_i_b(ji,nlay_i+1) = 0._wp
	444	END DO
[825]	445
[4634]	446	! Iterative procedure
	447	DO iter = 1, num_iter_max
	448	DO ji = kideb, kiut
	449	IF( zf_tt(ji) < 0._wp ) THEN
[825]	450
[4634]	451	! New bottom ice salinity (Cox & Weeks, JGR88 )
	452	!--- zswi1 if dh/dt < 2.0e-8
	453	!--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
	454	!--- zswi2 if dh/dt > 3.6e-7
	455	zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) )
	456	zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
	457	zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
	458	zswi1 = 1. - zswi2 * zswi12
	459	zfracs = MIN ( zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) &
	460	& + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 )
[825]	461
[4634]	462	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
[825]	463
[4634]	464	s_i_new(ji) = zswitch_sal * zfracs * sss_m(ii,ij) & ! New ice salinity
	465	+ ( 1. - zswitch_sal ) * sm_i_b(ji)
	466	! New ice growth
	467	ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K)
[825]	468
[4634]	469	zt_i_new = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)
	470
	471	zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0)
	472	& - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) ) &
	473	& + rcp * ( ztmelts-rtt )
	474
	475	zEw = rcp * ( t_bo_b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)
	476
	477	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
	478
	479	dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
	480
	481	q_i_b(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0)
	482
	483	ENDIF ! fc_bo_i
	484	END DO ! ji
	485	END DO ! iter
	486
	487	! Contribution to Energy and Salt Fluxes
[825]	488	DO ji = kideb, kiut
[4634]	489	IF( zf_tt(ji) < 0._wp ) THEN
	490	! New ice growth
	491
	492	zfmdt = - rhoic * dh_i_bott(ji) ! Mass flux x time step (kg/m2, < 0)
	493
	494	! Contribution to heat flux to the ocean [W.m-2], >0
	495	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
	496	! Total heat flux used in this process [W.m-2]
	497	hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
	498
	499	! Contribution to salt flux ()
	500	sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_b(ji) * zfmdt * r1_rdtice
	501
	502	! Contribution to mass flux
	503	wfx_bog_1d(ji) = wfx_bog_1d(ji) + rhoic * a_i_b(ji) * dh_i_bott(ji) * r1_rdtice
	504
	505	! clem debug: variation of enthalpy (J.m-2)
	506	dq_i(ji) = dq_i(ji) + dh_i_bott(ji) * q_i_b(ji,nlay_i+1)
	507
	508	! update heat content (J.m-2) and layer thickness
	509	qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_b(ji,nlay_i+1)
	510	h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
[825]	511	ENDIF
	512	END DO
	513
[4634]	514	!----------------
	515	! 4.2 Basal melt
	516	!----------------
	517	zdeltah(:,:) = 0._wp ! important
[825]	518	DO jk = nlay_i, 1, -1
	519	DO ji = kideb, kiut
[4634]	520	IF( zf_tt(ji) >= 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared from surface melting
[825]	521
[4634]	522	ztmelts = - tmut * s_i_b(ji,jk) + rtt ! Melting point of layer jk (K)
[825]	523
[4634]	524	IF( t_i_b(ji,jk) >= ztmelts ) THEN !!! Internal melting
	525	zintermelt(ji) = 1._wp
[825]	526
[4634]	527	zEi = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)
[825]	528
[4634]	529	!!zEw = rcp * ( t_i_b(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_b (J/kg, <0)
[4045]	530
[4634]	531	zdE = 0._wp ! Specific enthalpy difference (J/kg, <0)
	532	! set up at 0 since no energy is needed to melt water...(it is already melted)
[4045]	533
[4634]	534	zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing
	535	! this should normally not happen, but sometimes, heat diffusion leads to this
[825]	536
[4634]	537	dh_i_bott (ji) = dh_i_bott(ji) + zdeltah(ji,jk)
[825]	538
[4634]	539	zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0
[825]	540
[4634]	541	! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
	542	hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_b(ji) * zEi * r1_rdtice
[825]	543
[4634]	544	! clem debug: variation of enthalpy (J.m-2)
	545	dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)
[825]	546
[4634]	547	! update heat content (J.m-2) and layer thickness
	548	qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
	549	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
[825]	550
[4634]	551	ELSE !!! Basal melting
[825]	552
[4634]	553	zEi = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)
[825]	554
[4634]	555	zEw = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0)
[825]	556
[4634]	557	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
	558
	559	zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0)
	560
	561	zdeltah(ji,jk) = - zfmdt / rhoic ! Gross thickness change
	562
	563	zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
	564
	565	zq_bo(ji) = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
	566
	567	dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) ! Update basal melt
	568
	569	zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0
	570
	571	zQm = zfmdt * zEw ! Heat exchanged with ocean
	572
	573	! Contribution to heat flux to the ocean [W.m-2], <0
	574	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
	575
	576	! clem debug: variation of enthalpy (J.m-2)
	577	dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)
	578
	579	! update heat content (J.m-2) and layer thickness
	580	qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
	581	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
	582	ENDIF
	583
	584	! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
	585	sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
	586
	587	! Total heat flux used in this process [W.m-2]
	588	hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
	589
	590	! Contribution to mass flux
	591	wfx_bom_1d(ji) = wfx_bom_1d(ji) + rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
	592
	593	ENDIF
	594	END DO ! ji
	595	END DO ! jk
	596
	597	!------------------------------------------------------------------------------!
	598	! Excessive ablation in a 1-category model
	599	! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15)
	600	!------------------------------------------------------------------------------!
	601	! ??? keep ???
	602	! clem bug: I think this should be included above, so we would not have to
	603	! track heat/salt/mass fluxes backwards
	604	IF( jpl == 1 ) THEN
	605	DO ji = kideb, kiut
	606	IF( zf_tt(ji) >= 0._wp ) THEN
	607	zdh = MAX( hmelt , dh_i_bott(ji) )
	608	zdvres = zdh - dh_i_bott(ji) ! >=0
	609	dh_i_bott(ji) = zdh
	610
	611	! excessive energy is sent to lateral ablation
	612	zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_b(ji) - epsi20 ) )
	613	zq_1cat(ji) = zinda * rhoic * lfus * at_i_b(ji) / MAX( 1._wp - at_i_b(ji) , epsi20 ) * zdvres ! J.m-2 >=0
	614
	615	! correct salt and mass fluxes
	616	sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation
	617	wfx_bom_1d(ji) = wfx_bom_1d(ji) + rhoic * a_i_b(ji) * zdvres * r1_rdtice
	618	ENDIF
	619	END DO
	620	ENDIF
	621
[825]	622	!-------------------------------------------
[4634]	623	! Update temperature, energy
[825]	624	!-------------------------------------------
	625	DO ji = kideb, kiut
[4634]	626	ht_i_b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_bott(ji) )
	627	END DO
[825]	628
[4634]	629	!-------------------------------------------
	630	! 5. What to do with remaining energy
	631	!-------------------------------------------
	632	! If heat still available for melting and snow remains, then melt more snow
	633	!-------------------------------------------
	634	zdeltah(:,:) = 0._wp ! important
	635	DO ji = kideb, kiut
	636	zq_rema(ji) = zq_su(ji) + zq_bo(ji)
	637	! zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_b(ji) ) ) ! =1 if snow
	638	! zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zq_s(ji) + epsi20 ) )
	639	! zdeltah (ji,1) = - zindh * zindq * zq_rema(ji) / MAX( zq_s(ji), epsi20 )
	640	! zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_b(ji) ) ) ! bound melting
	641	! zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1)
	642	! dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1)
	643	! ht_s_b (ji) = ht_s_b(ji) + zdeltah(ji,1)
	644	!
	645	! zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji) ! update available heat (J.m-2)
	646	! ! Heat flux associated with snow melt
	647	! hfx_snw_1d(ji) = hfx_snw_1d(ji) + zdeltah(ji,1) * a_i_b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (<0)
	648	! ! heat used to melt snow
	649	! hfx_tot_1d(ji) = hfx_tot_1d(ji) - zdeltah(ji,1) * a_i_b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0)
	650	! ! Contribution to mass flux
	651	! wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdeltah(ji,1) * r1_rdtice
	652	! ! clem debug: variation of enthalpy (J.m-2)
	653	! dq_s(ji) = dq_s(ji) + zdeltah(ji,1) * q_s_b(ji,1)
	654	!
	655	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
	656	! Remaining heat flux (W.m-2) is sent to the ocean heat budget
	657	hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_b(ji) ) * r1_rdtice
[825]	658
[4634]	659	IF( ln_nicep .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
	660	END DO
	661
[921]	662	!
	663	!------------------------------------------------------------------------------\|
	664	! 6) Snow-Ice formation \|
	665	!------------------------------------------------------------------------------\|
[1572]	666	! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
	667	! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
[825]	668	DO ji = kideb, kiut
[1572]	669	!
[4634]	670	dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic ) )
[825]	671
[4634]	672	ht_i_b(ji) = ht_i_b(ji) + dh_snowice(ji)
	673	ht_s_b(ji) = ht_s_b(ji) - dh_snowice(ji)
[825]	674
[4634]	675	! Salinity of snow ice
	676	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
	677	zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_b(ji)
[825]	678
	679	! entrapment during snow ice formation
[4634]	680	! new salinity difference stored (to be used in limthd_ent.F90)
[4045]	681	IF ( num_sal == 2 ) THEN
[4634]	682	i_ice_switch = MAX( 0._wp , SIGN( 1._wp , ht_i_b(ji) - epsi10 ) )
[4045]	683	! salinity dif due to snow-ice formation
[4634]	684	dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( ht_i_b(ji), epsi10 ) * i_ice_switch
[4045]	685	! salinity dif due to bottom growth
[4634]	686	IF ( zf_tt(ji) < 0._wp ) THEN
	687	dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( ht_i_b(ji), epsi10 ) * i_ice_switch
[4045]	688	ENDIF
	689	ENDIF
[825]	690
[4634]	691	! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
	692	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
	693	zfmdt = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0._wp ) ! <0
	694	zsstK = sst_m(ii,ij) + rt0
	695	zEw = rcp * ( zsstK - rt0 )
	696	zQm = zfmdt * zEw
	697
	698	! Contribution to heat flux
	699	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
[825]	700
[4634]	701	! Contribution to salt flux
	702	sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_b(ji) * zfmdt * r1_rdtice
	703
	704	! Contribution to mass flux
	705	! All snow is thrown in the ocean, and seawater is taken to replace the volume
	706	wfx_sni_1d(ji) = wfx_sni_1d(ji) + a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice
	707	wfx_snw_1d(ji) = wfx_snw_1d(ji) - a_i_b(ji) * dh_snowice(ji) * rhosn * r1_rdtice
	708
	709	! clem debug: variation of enthalpy (J.m-2)
	710	dq_s(ji) = dq_s(ji) - dh_snowice(ji) * q_s_b(ji,1)
	711	dq_i(ji) = dq_i(ji) + dh_snowice(ji) * q_s_b(ji,1) + zfmdt * zEw
	712
	713	! update heat content (J.m-2) and layer thickness
	714	qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_b(ji,1) + zfmdt * zEw
	715	h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
	716
	717	! Total ablation (to debug)
[2715]	718	IF( ht_i_b(ji) <= 0._wp ) a_i_b(ji) = 0._wp
[825]	719
[4634]	720	END DO !ji
[4045]	721
[2715]	722	!
[4634]	723	!-------------------------------------------
	724	! Update temperature, energy
	725	!-------------------------------------------
	726	!clem bug: we should take snow into account here
	727	DO ji = kideb, kiut
	728	zindh = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) )
	729	t_su_b(ji) = zindh * t_su_b(ji) + ( 1.0 - zindh ) * rtt
	730	END DO ! ji
	731
	732	DO jk = 1, nlay_s
	733	DO ji = kideb,kiut
	734	! mask enthalpy
	735	zinda = MAX( 0._wp , SIGN( 1._wp, - ht_s_b(ji) ) )
	736	q_s_b(ji,jk) = ( 1.0 - zinda ) * q_s_b(ji,jk)
	737	! recalculate t_s_b from q_s_b
	738	t_s_b(ji,jk) = rtt + ( 1._wp - zinda ) * ( - q_s_b(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
	739	END DO
	740	END DO
	741
	742	CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema )
	743	CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
	744	CALL wrk_dealloc( jpij, zintermelt )
	745	CALL wrk_dealloc( jpij, jkmax, zdeltah, zh_i )
	746	CALL wrk_dealloc( jpij, icount )
[2715]	747	!
[4045]	748	!
[921]	749	END SUBROUTINE lim_thd_dh
[1572]	750
[825]	751	#else
[1572]	752	!!----------------------------------------------------------------------
	753	!! Default option NO LIM3 sea-ice model
	754	!!----------------------------------------------------------------------
[825]	755	CONTAINS
	756	SUBROUTINE lim_thd_dh ! Empty routine
	757	END SUBROUTINE lim_thd_dh
	758	#endif
[1572]	759
	760	!!======================================================================
[921]	761	END MODULE limthd_dh

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

Download in other formats: