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 branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 5047

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

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

Download in other formats: