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_dif.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 4741

Last change on this file since 4741 was 4688, checked in by clem, 10 years ago
new version of LIM3 with perfect conservation of heat, see ticket #1352
Property svn:keywords set to `Id`
File size: 40.0 KB

Rev	Line
[825]	1	MODULE limthd_dif
	2	!!======================================================================
	3	!! * MODULE limthd_dif *
	4	!! heat diffusion in sea ice
	5	!! computation of surface and inner T
	6	!!======================================================================
[2715]	7	!! History : LIM ! 02-2003 (M. Vancoppenolle) original 1D code
	8	!! ! 06-2005 (M. Vancoppenolle) 3d version
	9	!! ! 11-2006 (X Fettweis) Vectorization by Xavier
	10	!! ! 04-2007 (M. Vancoppenolle) Energy conservation
	11	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
[4161]	12	!! - ! 2012-05 (C. Rousset) add penetration solar flux
[825]	13	!!----------------------------------------------------------------------
[2528]	14	#if defined key_lim3
	15	!!----------------------------------------------------------------------
	16	!! 'key_lim3' LIM3 sea-ice model
	17	!!----------------------------------------------------------------------
[3625]	18	USE par_oce ! ocean parameters
	19	USE phycst ! physical constants (ocean directory)
	20	USE ice ! LIM-3 variables
	21	USE par_ice ! LIM-3 parameters
	22	USE thd_ice ! LIM-3: thermodynamics
	23	USE in_out_manager ! I/O manager
	24	USE lib_mpp ! MPP library
	25	USE wrk_nemo ! work arrays
	26	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[4688]	27	USE cpl_oasis3, ONLY : lk_cpl
[921]	28
[825]	29	IMPLICIT NONE
	30	PRIVATE
	31
[2528]	32	PUBLIC lim_thd_dif ! called by lim_thd
[825]	33
[4688]	34	REAL(wp) :: epsi10 = 1.e-10_wp !
[825]	35	!!----------------------------------------------------------------------
[4161]	36	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	37	!! $Id$
[2591]	38	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	39	!!----------------------------------------------------------------------
	40	CONTAINS
	41
[4688]	42	SUBROUTINE lim_thd_dif( kideb , kiut )
[921]	43	!!------------------------------------------------------------------
	44	!! * ROUTINE lim_thd_dif *
	45	!! ** Purpose :
	46	!! This routine determines the time evolution of snow and sea-ice
	47	!! temperature profiles.
	48	!! ** Method :
	49	!! This is done by solving the heat equation diffusion with
	50	!! a Neumann boundary condition at the surface and a Dirichlet one
	51	!! at the bottom. Solar radiation is partially absorbed into the ice.
	52	!! The specific heat and thermal conductivities depend on ice salinity
	53	!! and temperature to take into account brine pocket melting. The
	54	!! numerical
	55	!! scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid
	56	!! in the ice and snow system.
	57	!!
	58	!! The successive steps of this routine are
	59	!! 1. Thermal conductivity at the interfaces of the ice layers
	60	!! 2. Internal absorbed radiation
	61	!! 3. Scale factors due to non-uniform grid
	62	!! 4. Kappa factors
	63	!! Then iterative procedure begins
	64	!! 5. specific heat in the ice
	65	!! 6. eta factors
	66	!! 7. surface flux computation
	67	!! 8. tridiagonal system terms
	68	!! 9. solving the tridiagonal system with Gauss elimination
	69	!! Iterative procedure ends according to a criterion on evolution
	70	!! of temperature
	71	!!
	72	!! ** Arguments :
	73	!! kideb , kiut : Starting and ending points on which the
	74	!! the computation is applied
	75	!!
	76	!! ** Inputs / Ouputs : (global commons)
	77	!! surface temperature : t_su_b
	78	!! ice/snow temperatures : t_i_b, t_s_b
	79	!! ice salinities : s_i_b
	80	!! number of layers in the ice/snow: nlay_i, nlay_s
	81	!! profile of the ice/snow layers : z_i, z_s
	82	!! total ice/snow thickness : ht_i_b, ht_s_b
	83	!!
	84	!! ** External :
	85	!!
	86	!! ** References :
	87	!!
	88	!! ** History :
	89	!! (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium
	90	!! (06-2005) Martin Vancoppenolle, 3d version
	91	!! (11-2006) Vectorized by Xavier Fettweis (UCL-ASTR)
	92	!! (04-2007) Energy conservation tested by M. Vancoppenolle
	93	!!------------------------------------------------------------------
[4688]	94	INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied
[825]	95
[921]	96	!! * Local variables
[3294]	97	INTEGER :: ji ! spatial loop index
	98	INTEGER :: ii, ij ! temporary dummy loop index
	99	INTEGER :: numeq ! current reference number of equation
	100	INTEGER :: layer ! vertical dummy loop index
	101	INTEGER :: nconv ! number of iterations in iterative procedure
	102	INTEGER :: minnumeqmin, maxnumeqmax
[4688]	103	INTEGER, POINTER, DIMENSION(:) :: numeqmin ! reference number of top equation
	104	INTEGER, POINTER, DIMENSION(:) :: numeqmax ! reference number of bottom equation
	105	INTEGER, POINTER, DIMENSION(:) :: isnow ! switch for presence (1) or absence (0) of snow
[3294]	106	REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system
	107	REAL(wp) :: zg1 = 2._wp !
	108	REAL(wp) :: zgamma = 18009._wp ! for specific heat
	109	REAL(wp) :: zbeta = 0.117_wp ! for thermal conductivity (could be 0.13)
	110	REAL(wp) :: zraext_s = 1.e+8_wp ! extinction coefficient of radiation in the snow
	111	REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity
[4688]	112	REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered as 0°C
[2715]	113	REAL(wp) :: ztmelt_i ! ice melting temperature
	114	REAL(wp) :: zerritmax ! current maximal error on temperature
[4688]	115	REAL(wp), POINTER, DIMENSION(:) :: ztfs ! ice melting point
	116	REAL(wp), POINTER, DIMENSION(:) :: ztsuold ! old surface temperature (before the iterative procedure )
	117	REAL(wp), POINTER, DIMENSION(:) :: ztsuoldit ! surface temperature at previous iteration
	118	REAL(wp), POINTER, DIMENSION(:) :: zh_i ! ice layer thickness
	119	REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness
	120	REAL(wp), POINTER, DIMENSION(:) :: zfsw ! solar radiation absorbed at the surface
	121	REAL(wp), POINTER, DIMENSION(:) :: zf ! surface flux function
	122	REAL(wp), POINTER, DIMENSION(:) :: dzf ! derivative of the surface flux function
	123	REAL(wp), POINTER, DIMENSION(:) :: zerrit ! current error on temperature
	124	REAL(wp), POINTER, DIMENSION(:) :: zdifcase ! case of the equation resolution (1->4)
	125	REAL(wp), POINTER, DIMENSION(:) :: zftrice ! solar radiation transmitted through the ice
	126	REAL(wp), POINTER, DIMENSION(:) :: zihic, zhsu
	127	REAL(wp), POINTER, DIMENSION(:,:) :: ztcond_i ! Ice thermal conductivity
	128	REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_i ! Radiation transmitted through the ice
	129	REAL(wp), POINTER, DIMENSION(:,:) :: zradab_i ! Radiation absorbed in the ice
	130	REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_i ! Kappa factor in the ice
	131	REAL(wp), POINTER, DIMENSION(:,:) :: ztiold ! Old temperature in the ice
	132	REAL(wp), POINTER, DIMENSION(:,:) :: zeta_i ! Eta factor in the ice
	133	REAL(wp), POINTER, DIMENSION(:,:) :: ztitemp ! Temporary temperature in the ice to check the convergence
	134	REAL(wp), POINTER, DIMENSION(:,:) :: zspeche_i ! Ice specific heat
	135	REAL(wp), POINTER, DIMENSION(:,:) :: z_i ! Vertical cotes of the layers in the ice
	136	REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_s ! Radiation transmited through the snow
	137	REAL(wp), POINTER, DIMENSION(:,:) :: zradab_s ! Radiation absorbed in the snow
	138	REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_s ! Kappa factor in the snow
	139	REAL(wp), POINTER, DIMENSION(:,:) :: zeta_s ! Eta factor in the snow
	140	REAL(wp), POINTER, DIMENSION(:,:) :: ztstemp ! Temporary temperature in the snow to check the convergence
	141	REAL(wp), POINTER, DIMENSION(:,:) :: ztsold ! Temporary temperature in the snow
	142	REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! Vertical cotes of the layers in the snow
	143	REAL(wp), POINTER, DIMENSION(:,:) :: zindterm ! Independent term
	144	REAL(wp), POINTER, DIMENSION(:,:) :: zindtbis ! temporary independent term
	145	REAL(wp), POINTER, DIMENSION(:,:) :: zdiagbis
	146	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! tridiagonal system terms
	147	! diag errors on heat
	148	REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini
	149	REAL(wp) :: zhfx_err
[3625]	150	!!------------------------------------------------------------------
[3610]	151	!
[4688]	152	CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow )
	153	CALL wrk_alloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw )
	154	CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu )
	155	CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0)
	156	CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart=0)
	157	CALL wrk_alloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis )
	158	CALL wrk_alloc( jpij, jkmax+2, 3, ztrid )
	159
	160	CALL wrk_alloc( jpij, zdq, zq_ini )
	161
	162	! --- diag error on heat diffusion - PART 1 --- !
	163	zdq(:) = 0._wp ; zq_ini(:) = 0._wp
	164	DO ji = kideb, kiut
	165	zq_ini(ji) = ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + &
	166	& SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )
	167	END DO
	168
[921]	169	!------------------------------------------------------------------------------!
	170	! 1) Initialization !
	171	!------------------------------------------------------------------------------!
[4688]	172	! clem clean: replace just ztfs by rtt
[921]	173	DO ji = kideb , kiut
	174	! is there snow or not
[4161]	175	isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) ) )
[921]	176	! surface temperature of fusion
[4161]	177	ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt
[921]	178	! layer thickness
[4161]	179	zh_i(ji) = ht_i_b(ji) / REAL( nlay_i )
	180	zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
[921]	181	END DO
[825]	182
[921]	183	!--------------------
	184	! Ice / snow layers
	185	!--------------------
[825]	186
[2715]	187	z_s(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st snow layer
	188	z_i(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st ice layer
[825]	189
[2715]	190	DO layer = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer
[921]	191	DO ji = kideb , kiut
[4161]	192	z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s )
[921]	193	END DO
	194	END DO
[825]	195
[2715]	196	DO layer = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer
[921]	197	DO ji = kideb , kiut
[4161]	198	z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i )
[921]	199	END DO
	200	END DO
	201	!
	202	!------------------------------------------------------------------------------\|
	203	! 2) Radiations \|
	204	!------------------------------------------------------------------------------\|
	205	!
	206	!-------------------
	207	! Computation of i0
	208	!-------------------
	209	! i0 describes the fraction of solar radiation which does not contribute
	210	! to the surface energy budget but rather penetrates inside the ice.
	211	! We assume that no radiation is transmitted through the snow
	212	! If there is no no snow
	213	! zfsw = (1-i0).qsr_ice is absorbed at the surface
	214	! zftrice = io.qsr_ice is below the surface
[4688]	215	! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
[825]	216
[921]	217	DO ji = kideb , kiut
	218	! switches
[4161]	219	isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) ) )
[921]	220	! hs > 0, isnow = 1
[2715]	221	zhsu (ji) = hnzst ! threshold for the computation of i0
	222	zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) )
[825]	223
[4161]	224	i0(ji) = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
[921]	225	!fr1_i0_1d = i0 for a thin ice surface
	226	!fr1_i0_2d = i0 for a thick ice surface
	227	! a function of the cloud cover
	228	!
	229	!i0(ji) = (1.0-FLOAT(isnow(ji)))3.0/(100ht_s_b(ji)+10.0)
	230	!formula used in Cice
	231	END DO
[825]	232
[921]	233	!-------------------------------------------------------
	234	! Solar radiation absorbed / transmitted at the surface
	235	! Derivative of the non solar flux
	236	!-------------------------------------------------------
	237	DO ji = kideb , kiut
[2715]	238	zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface
	239	zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer
	240	dzf (ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux
[921]	241	END DO
[825]	242
[921]	243	!---------------------------------------------------------
	244	! Transmission - absorption of solar radiation in the ice
	245	!---------------------------------------------------------
[825]	246
[2715]	247	DO ji = kideb, kiut ! snow initialization
	248	zradtr_s(ji,0) = zftrice(ji) ! radiation penetrating through snow
[921]	249	END DO
[825]	250
[2715]	251	DO layer = 1, nlay_s ! Radiation through snow
	252	DO ji = kideb, kiut
	253	! ! radiation transmitted below the layer-th snow layer
	254	zradtr_s(ji,layer) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,layer) ) ) )
	255	! ! radiation absorbed by the layer-th snow layer
[921]	256	zradab_s(ji,layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer)
	257	END DO
	258	END DO
[825]	259
[2715]	260	DO ji = kideb, kiut ! ice initialization
[4161]	261	zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * REAL( isnow(ji) ) + zftrice(ji) * REAL( 1 - isnow(ji) )
[921]	262	END DO
[825]	263
[2715]	264	DO layer = 1, nlay_i ! Radiation through ice
	265	DO ji = kideb, kiut
	266	! ! radiation transmitted below the layer-th ice layer
	267	zradtr_i(ji,layer) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,layer) ) ) )
	268	! ! radiation absorbed by the layer-th ice layer
[921]	269	zradab_i(ji,layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer)
	270	END DO
	271	END DO
[825]	272
[2715]	273	DO ji = kideb, kiut ! Radiation transmitted below the ice
[4688]	274	!!!ftr_ice_1d(ji) = ftr_ice_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) ! clem modif
	275	ftr_ice_1d(ji) = zradtr_i(ji,nlay_i)
[921]	276	END DO
[834]	277
[921]	278	!
	279	!------------------------------------------------------------------------------\|
	280	! 3) Iterative procedure begins \|
	281	!------------------------------------------------------------------------------\|
	282	!
[2715]	283	DO ji = kideb, kiut ! Old surface temperature
	284	ztsuold (ji) = t_su_b(ji) ! temperature at the beg of iter pr.
	285	ztsuoldit(ji) = t_su_b(ji) ! temperature at the previous iter
[4688]	286	t_su_b (ji) = MIN( t_su_b(ji), ztfs(ji) - ztsu_err ) ! necessary
[2715]	287	zerrit (ji) = 1000._wp ! initial value of error
[921]	288	END DO
[825]	289
[2715]	290	DO layer = 1, nlay_s ! Old snow temperature
[921]	291	DO ji = kideb , kiut
[2715]	292	ztsold(ji,layer) = t_s_b(ji,layer)
[921]	293	END DO
	294	END DO
[825]	295
[2715]	296	DO layer = 1, nlay_i ! Old ice temperature
[921]	297	DO ji = kideb , kiut
[2715]	298	ztiold(ji,layer) = t_i_b(ji,layer)
[921]	299	END DO
	300	END DO
[825]	301
[2715]	302	nconv = 0 ! number of iterations
	303	zerritmax = 1000._wp ! maximal value of error on all points
[825]	304
[2715]	305	DO WHILE ( zerritmax > maxer_i_thd .AND. nconv < nconv_i_thd )
[921]	306	!
[2715]	307	nconv = nconv + 1
	308	!
[921]	309	!------------------------------------------------------------------------------\|
	310	! 4) Sea ice thermal conductivity \|
	311	!------------------------------------------------------------------------------\|
	312	!
[2715]	313	IF( thcon_i_swi == 0 ) THEN ! Untersteiner (1964) formula
[921]	314	DO ji = kideb , kiut
[4333]	315	ztcond_i(ji,0) = rcdic + zbeta*s_i_b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt)
[921]	316	ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin)
	317	END DO
	318	DO layer = 1, nlay_i-1
	319	DO ji = kideb , kiut
[2777]	320	ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) / &
[4333]	321	MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)
[921]	322	ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin)
	323	END DO
	324	END DO
	325	ENDIF
[825]	326
[2715]	327	IF( thcon_i_swi == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
[921]	328	DO ji = kideb , kiut
[4333]	329	ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt ) &
[2715]	330	& - 0.011_wp * ( t_i_b(ji,1) - rtt )
	331	ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
[921]	332	END DO
[2715]	333	DO layer = 1, nlay_i-1
	334	DO ji = kideb , kiut
	335	ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) &
[4333]	336	& / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) &
[2715]	337	& - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )
	338	ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin )
	339	END DO
	340	END DO
	341	DO ji = kideb , kiut
[4333]	342	ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt) &
[2715]	343	& - 0.011_wp * ( t_bo_b(ji) - rtt )
	344	ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
	345	END DO
[921]	346	ENDIF
	347	!
	348	!------------------------------------------------------------------------------\|
	349	! 5) kappa factors \|
	350	!------------------------------------------------------------------------------\|
	351	!
	352	DO ji = kideb, kiut
[825]	353
[921]	354	!-- Snow kappa factors
[4333]	355	zkappa_s(ji,0) = rcdsn / MAX(epsi10,zh_s(ji))
	356	zkappa_s(ji,nlay_s) = rcdsn / MAX(epsi10,zh_s(ji))
[921]	357	END DO
[825]	358
[921]	359	DO layer = 1, nlay_s-1
	360	DO ji = kideb , kiut
	361	zkappa_s(ji,layer) = 2.0 * rcdsn / &
[4333]	362	MAX(epsi10,2.0*zh_s(ji))
[921]	363	END DO
	364	END DO
[825]	365
[921]	366	DO layer = 1, nlay_i-1
	367	DO ji = kideb , kiut
	368	!-- Ice kappa factors
	369	zkappa_i(ji,layer) = 2.0*ztcond_i(ji,layer)/ &
[4333]	370	MAX(epsi10,2.0*zh_i(ji))
[921]	371	END DO
	372	END DO
[825]	373
[921]	374	DO ji = kideb , kiut
[4333]	375	zkappa_i(ji,0) = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji))
	376	zkappa_i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji))
[921]	377	!-- Interface
[4333]	378	zkappa_s(ji,nlay_s) = 2.0rcdsnztcond_i(ji,0)/MAX(epsi10, &
[921]	379	(ztcond_i(ji,0)zh_s(ji) + rcdsnzh_i(ji)))
[4161]	380	zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) &
	381	+ zkappa_i(ji,0)*REAL( 1 - isnow(ji) )
[921]	382	END DO
	383	!
	384	!------------------------------------------------------------------------------\|
	385	! 6) Sea ice specific heat, eta factors \|
	386	!------------------------------------------------------------------------------\|
	387	!
	388	DO layer = 1, nlay_i
	389	DO ji = kideb , kiut
	390	ztitemp(ji,layer) = t_i_b(ji,layer)
	391	zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ &
[4333]	392	MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10)
[921]	393	zeta_i(ji,layer) = rdt_ice / MAX(rhoiczspeche_i(ji,layer)zh_i(ji), &
[4333]	394	epsi10)
[921]	395	END DO
	396	END DO
[825]	397
[921]	398	DO layer = 1, nlay_s
	399	DO ji = kideb , kiut
	400	ztstemp(ji,layer) = t_s_b(ji,layer)
[4333]	401	zeta_s(ji,layer) = rdt_ice / MAX(rhosncpiczh_s(ji),epsi10)
[921]	402	END DO
	403	END DO
	404	!
	405	!------------------------------------------------------------------------------\|
	406	! 7) surface flux computation \|
	407	!------------------------------------------------------------------------------\|
	408	!
	409	DO ji = kideb , kiut
	410	! update of the non solar flux according to the update in T_su
[4688]	411	qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) )
[825]	412
[921]	413	! update incoming flux
	414	zf(ji) = zfsw(ji) & ! net absorbed solar radiation
[4688]	415	+ qns_ice_1d(ji) ! non solar total flux
[921]	416	! (LWup, LWdw, SH, LH)
	417	END DO
[825]	418
[921]	419	!
	420	!------------------------------------------------------------------------------\|
	421	! 8) tridiagonal system terms \|
	422	!------------------------------------------------------------------------------\|
	423	!
	424	!!layer denotes the number of the layer in the snow or in the ice
	425	!!numeq denotes the reference number of the equation in the tridiagonal
	426	!!system, terms of tridiagonal system are indexed as following :
	427	!!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
[825]	428
[921]	429	!!ice interior terms (top equation has the same form as the others)
	430
	431	DO numeq=1,jkmax+2
	432	DO ji = kideb , kiut
	433	ztrid(ji,numeq,1) = 0.
	434	ztrid(ji,numeq,2) = 0.
	435	ztrid(ji,numeq,3) = 0.
	436	zindterm(ji,numeq)= 0.
	437	zindtbis(ji,numeq)= 0.
	438	zdiagbis(ji,numeq)= 0.
	439	ENDDO
	440	ENDDO
	441
	442	DO numeq = nlay_s + 2, nlay_s + nlay_i
	443	DO ji = kideb , kiut
	444	layer = numeq - nlay_s - 1
	445	ztrid(ji,numeq,1) = - zeta_i(ji,layer)*zkappa_i(ji,layer-1)
	446	ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,layer)*(zkappa_i(ji,layer-1) + &
	447	zkappa_i(ji,layer))
	448	ztrid(ji,numeq,3) = - zeta_i(ji,layer)*zkappa_i(ji,layer)
	449	zindterm(ji,numeq) = ztiold(ji,layer) + zeta_i(ji,layer)* &
	450	zradab_i(ji,layer)
	451	END DO
	452	ENDDO
	453
	454	numeq = nlay_s + nlay_i + 1
	455	DO ji = kideb , kiut
[825]	456	!!ice bottom term
	457	ztrid(ji,numeq,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)
	458	ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,nlay_i)( zkappa_i(ji,nlay_i)zg1 &
[921]	459	+ zkappa_i(ji,nlay_i-1) )
[825]	460	ztrid(ji,numeq,3) = 0.0
	461	zindterm(ji,numeq) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i)* &
[921]	462	( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 &
	463	* t_bo_b(ji) )
	464	ENDDO
[825]	465
	466
[921]	467	DO ji = kideb , kiut
[825]	468	IF ( ht_s_b(ji).gt.0.0 ) THEN
[921]	469	!
	470	!------------------------------------------------------------------------------\|
	471	! snow-covered cells \|
	472	!------------------------------------------------------------------------------\|
	473	!
	474	!!snow interior terms (bottom equation has the same form as the others)
	475	DO numeq = 3, nlay_s + 1
	476	layer = numeq - 1
	477	ztrid(ji,numeq,1) = - zeta_s(ji,layer)*zkappa_s(ji,layer-1)
	478	ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,layer)*( zkappa_s(ji,layer-1) + &
	479	zkappa_s(ji,layer) )
	480	ztrid(ji,numeq,3) = - zeta_s(ji,layer)*zkappa_s(ji,layer)
	481	zindterm(ji,numeq) = ztsold(ji,layer) + zeta_s(ji,layer)* &
	482	zradab_s(ji,layer)
	483	END DO
[825]	484
[921]	485	!!case of only one layer in the ice (ice equation is altered)
	486	IF ( nlay_i.eq.1 ) THEN
	487	ztrid(ji,nlay_s+2,3) = 0.0
	488	zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* &
	489	t_bo_b(ji)
	490	ENDIF
[834]	491
[921]	492	IF ( t_su_b(ji) .LT. rtt ) THEN
[825]	493
[921]	494	!------------------------------------------------------------------------------\|
	495	! case 1 : no surface melting - snow present \|
	496	!------------------------------------------------------------------------------\|
	497	zdifcase(ji) = 1.0
	498	numeqmin(ji) = 1
	499	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	500
[921]	501	!!surface equation
	502	ztrid(ji,1,1) = 0.0
	503	ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0)
	504	ztrid(ji,1,3) = zg1s*zkappa_s(ji,0)
	505	zindterm(ji,1) = dzf(ji)*t_su_b(ji) - zf(ji)
[825]	506
[921]	507	!!first layer of snow equation
	508	ztrid(ji,2,1) = - zkappa_s(ji,0)zg1szeta_s(ji,1)
	509	ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)(zkappa_s(ji,1) + zkappa_s(ji,0)zg1s)
	510	ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1)
	511	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)
[825]	512
[921]	513	ELSE
	514	!
	515	!------------------------------------------------------------------------------\|
	516	! case 2 : surface is melting - snow present \|
	517	!------------------------------------------------------------------------------\|
	518	!
	519	zdifcase(ji) = 2.0
	520	numeqmin(ji) = 2
	521	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	522
[921]	523	!!first layer of snow equation
	524	ztrid(ji,2,1) = 0.0
	525	ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + &
	526	zkappa_s(ji,0) * zg1s )
	527	ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1)
	528	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * &
	529	( zradab_s(ji,1) + &
	530	zkappa_s(ji,0) * zg1s * t_su_b(ji) )
	531	ENDIF
	532	ELSE
	533	!
	534	!------------------------------------------------------------------------------\|
	535	! cells without snow \|
	536	!------------------------------------------------------------------------------\|
	537	!
	538	IF (t_su_b(ji) .LT. rtt) THEN
	539	!
	540	!------------------------------------------------------------------------------\|
	541	! case 3 : no surface melting - no snow \|
	542	!------------------------------------------------------------------------------\|
	543	!
	544	zdifcase(ji) = 3.0
	545	numeqmin(ji) = nlay_s + 1
	546	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	547
[921]	548	!!surface equation
	549	ztrid(ji,numeqmin(ji),1) = 0.0
	550	ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1
	551	ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1
	552	zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_b(ji) - zf(ji)
[825]	553
[921]	554	!!first layer of ice equation
	555	ztrid(ji,numeqmin(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
	556	ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) &
	557	+ zkappa_i(ji,0) * zg1 )
	558	ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1)
	559	zindterm(ji,numeqmin(ji)+1)= ztiold(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)
[825]	560
[921]	561	!!case of only one layer in the ice (surface & ice equations are altered)
[825]	562
[921]	563	IF (nlay_i.eq.1) THEN
	564	ztrid(ji,numeqmin(ji),1) = 0.0
	565	ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*2.0
	566	ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*2.0
	567	ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0)2.0zeta_i(ji,1)
	568	ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,0)2.0 + &
	569	zkappa_i(ji,1))
	570	ztrid(ji,numeqmin(ji)+1,3) = 0.0
[825]	571
[921]	572	zindterm(ji,numeqmin(ji)+1) = ztiold(ji,1) + zeta_i(ji,1)* &
	573	( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji) )
	574	ENDIF
[825]	575
[921]	576	ELSE
[825]	577
[921]	578	!
	579	!------------------------------------------------------------------------------\|
	580	! case 4 : surface is melting - no snow \|
	581	!------------------------------------------------------------------------------\|
	582	!
	583	zdifcase(ji) = 4.0
	584	numeqmin(ji) = nlay_s + 2
	585	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	586
[921]	587	!!first layer of ice equation
	588	ztrid(ji,numeqmin(ji),1) = 0.0
	589	ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,1) + zkappa_i(ji,0) &
	590	zg1)
	591	ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)
	592	zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &
	593	zkappa_i(ji,0) * zg1 * t_su_b(ji) )
[825]	594
[921]	595	!!case of only one layer in the ice (surface & ice equations are altered)
	596	IF (nlay_i.eq.1) THEN
	597	ztrid(ji,numeqmin(ji),1) = 0.0
	598	ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,0)2.0 + &
	599	zkappa_i(ji,1))
	600	ztrid(ji,numeqmin(ji),3) = 0.0
	601	zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)* &
	602	(zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji)) &
	603	+ t_su_b(ji)zeta_i(ji,1)zkappa_i(ji,0)*2.0
	604	ENDIF
[825]	605
[921]	606	ENDIF
	607	ENDIF
[825]	608
[921]	609	END DO
[825]	610
[921]	611	!
	612	!------------------------------------------------------------------------------\|
	613	! 9) tridiagonal system solving \|
	614	!------------------------------------------------------------------------------\|
	615	!
[825]	616
[921]	617	! Solve the tridiagonal system with Gauss elimination method.
	618	! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
	619	! McGraw-Hill 1984.
[825]	620
[921]	621	maxnumeqmax = 0
	622	minnumeqmin = jkmax+4
[825]	623
[921]	624	DO ji = kideb , kiut
	625	zindtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji))
	626	zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2)
	627	minnumeqmin = MIN(numeqmin(ji),minnumeqmin)
	628	maxnumeqmax = MAX(numeqmax(ji),maxnumeqmax)
	629	END DO
	630
	631	DO layer = minnumeqmin+1, maxnumeqmax
	632	DO ji = kideb , kiut
	633	numeq = min(max(numeqmin(ji)+1,layer),numeqmax(ji))
	634	zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* &
	635	ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1)
	636	zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1)* &
	637	zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1)
	638	END DO
	639	END DO
	640
	641	DO ji = kideb , kiut
	642	! ice temperatures
	643	t_i_b(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))
	644	END DO
	645
	646	DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1
	647	DO ji = kideb , kiut
	648	layer = numeq - nlay_s - 1
	649	t_i_b(ji,layer) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &
	650	t_i_b(ji,layer+1))/zdiagbis(ji,numeq)
	651	END DO
	652	END DO
	653
	654	DO ji = kideb , kiut
[825]	655	! snow temperatures
[4205]	656	IF (ht_s_b(ji).GT.0._wp) &
[921]	657	t_s_b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
	658	* t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) &
[4161]	659	* MAX(0.0,SIGN(1.0,ht_s_b(ji)))
[825]	660
	661	! surface temperature
[4161]	662	isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) ) ) )
[2715]	663	ztsuoldit(ji) = t_su_b(ji)
[4161]	664	IF( t_su_b(ji) < ztfs(ji) ) &
	665	t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1) &
	666	& + REAL( 1 - isnow(ji) )*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))
[921]	667	END DO
	668	!
	669	!--------------------------------------------------------------------------
	670	! 10) Has the scheme converged ?, end of the iterative procedure \|
	671	!--------------------------------------------------------------------------
	672	!
	673	! check that nowhere it has started to melt
	674	! zerrit(ji) is a measure of error, it has to be under maxer_i_thd
	675	DO ji = kideb , kiut
[2715]	676	t_su_b(ji) = MAX( MIN( t_su_b(ji) , ztfs(ji) ) , 190._wp )
	677	zerrit(ji) = ABS( t_su_b(ji) - ztsuoldit(ji) )
[921]	678	END DO
[825]	679
[921]	680	DO layer = 1, nlay_s
	681	DO ji = kideb , kiut
[2715]	682	t_s_b(ji,layer) = MAX( MIN( t_s_b(ji,layer), rtt ), 190._wp )
	683	zerrit(ji) = MAX(zerrit(ji),ABS(t_s_b(ji,layer) - ztstemp(ji,layer)))
[921]	684	END DO
	685	END DO
[825]	686
[921]	687	DO layer = 1, nlay_i
	688	DO ji = kideb , kiut
[2715]	689	ztmelt_i = -tmut * s_i_b(ji,layer) + rtt
[4205]	690	t_i_b(ji,layer) = MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp)
[2715]	691	zerrit(ji) = MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer)))
[921]	692	END DO
	693	END DO
[825]	694
[921]	695	! Compute spatial maximum over all errors
[2715]	696	! note that this could be optimized substantially by iterating only the non-converging points
	697	zerritmax = 0._wp
	698	DO ji = kideb, kiut
	699	zerritmax = MAX( zerritmax, zerrit(ji) )
[921]	700	END DO
[2715]	701	IF( lk_mpp ) CALL mpp_max( zerritmax, kcom=ncomm_ice )
[825]	702
	703	END DO ! End of the do while iterative procedure
	704
[4333]	705	IF( ln_nicep .AND. lwp ) THEN
[1055]	706	WRITE(numout,*) ' zerritmax : ', zerritmax
	707	WRITE(numout,*) ' nconv : ', nconv
	708	ENDIF
[825]	709
[921]	710	!
[2715]	711	!-------------------------------------------------------------------------!
	712	! 11) Fluxes at the interfaces !
	713	!-------------------------------------------------------------------------!
[921]	714	DO ji = kideb, kiut
[3808]	715	! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)
[4688]	716	IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) )
[2715]	717	! ! surface ice conduction flux
[4161]	718	isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) ) )
	719	fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji)) &
	720	& - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_b(ji,1) - t_su_b(ji))
[2715]	721	! ! bottom ice conduction flux
	722	fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
[921]	723	END DO
[825]	724
[4688]	725	!-----------------------------------------
	726	! Heat flux used to warm/cool ice in W.m-2
	727	!-----------------------------------------
	728	DO ji = kideb, kiut
	729	IF( t_su_b(ji) < rtt ) THEN ! case T_su < 0degC
	730	hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji)
	731	ELSE ! case T_su = 0degC
	732	hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji)
	733	ENDIF
	734	END DO
	735
	736	! --- computes sea ice energy of melting compulsory for limthd_dh --- !
	737	CALL lim_thd_enmelt( kideb, kiut )
	738
	739	! --- diag error on heat diffusion - PART 2 --- !
	740	DO ji = kideb, kiut
	741	zdq(ji) = - zq_ini(ji) + ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + &
	742	& SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )
	743	zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice )
	744	hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err * a_i_b(ji)
	745	! --- correction of qns_ice and surface conduction flux --- !
	746	qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err
	747	fc_su (ji) = fc_su (ji) - zhfx_err
	748	! --- Heat flux at the ice surface in W.m-2 --- !
	749	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
	750	hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) )
	751	END DO
	752
	753	!
	754	CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow )
	755	CALL wrk_dealloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw )
	756	CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu )
	757	CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 )
	758	CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart = 0 )
	759	CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis )
	760	CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid )
	761	CALL wrk_dealloc( jpij, zdq, zq_ini )
	762
	763	END SUBROUTINE lim_thd_dif
	764
	765	SUBROUTINE lim_thd_enmelt( kideb, kiut )
	766	!!-----------------------------------------------------------------------
	767	!! * ROUTINE lim_thd_enmelt *
	768	!!
	769	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature
	770	!!
	771	!! ** Method : Formula (Bitz and Lipscomb, 1999)
	772	!!-------------------------------------------------------------------
	773	INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop
	774	!
	775	INTEGER :: ji, jk ! dummy loop indices
	776	REAL(wp) :: ztmelts, zindb ! local scalar
	777	!!-------------------------------------------------------------------
	778	!
	779	DO jk = 1, nlay_i ! Sea ice energy of melting
[921]	780	DO ji = kideb, kiut
[4688]	781	ztmelts = - tmut * s_i_b(ji,jk) + rtt
	782	zindb = MAX( 0._wp , SIGN( 1._wp , -(t_i_b(ji,jk) - rtt) - epsi10 ) )
	783	q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) ) &
	784	& + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) ) &
	785	& - rcp * ( ztmelts-rtt ) )
[921]	786	END DO
[4688]	787	END DO
	788	DO jk = 1, nlay_s ! Snow energy of melting
	789	DO ji = kideb, kiut
	790	q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
[921]	791	END DO
[4688]	792	END DO
[2715]	793	!
[4688]	794	END SUBROUTINE lim_thd_enmelt
[825]	795
	796	#else
[2715]	797	!!----------------------------------------------------------------------
	798	!! Dummy Module No LIM-3 sea-ice model
	799	!!----------------------------------------------------------------------
[825]	800	CONTAINS
	801	SUBROUTINE lim_thd_dif ! Empty routine
	802	END SUBROUTINE lim_thd_dif
	803	#endif
[2528]	804	!!======================================================================
[921]	805	END MODULE limthd_dif

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

Download in other formats: