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 @ 4780

Last change on this file since 4780 was 4765, checked in by rblod, 10 years ago
Compilation issue, see ticket #1379
Property svn:keywords set to `Id`
File size: 40.2 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
[4765]	335	ztcond_i(ji,layer) = rcdic + &
	336	& 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) &
	337	& / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) &
	338	& - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )
[2715]	339	ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin )
	340	END DO
	341	END DO
	342	DO ji = kideb , kiut
[4333]	343	ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt) &
[2715]	344	& - 0.011_wp * ( t_bo_b(ji) - rtt )
	345	ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
	346	END DO
[921]	347	ENDIF
	348	!
	349	!------------------------------------------------------------------------------\|
	350	! 5) kappa factors \|
	351	!------------------------------------------------------------------------------\|
	352	!
	353	DO ji = kideb, kiut
[825]	354
[921]	355	!-- Snow kappa factors
[4333]	356	zkappa_s(ji,0) = rcdsn / MAX(epsi10,zh_s(ji))
	357	zkappa_s(ji,nlay_s) = rcdsn / MAX(epsi10,zh_s(ji))
[921]	358	END DO
[825]	359
[921]	360	DO layer = 1, nlay_s-1
	361	DO ji = kideb , kiut
	362	zkappa_s(ji,layer) = 2.0 * rcdsn / &
[4333]	363	MAX(epsi10,2.0*zh_s(ji))
[921]	364	END DO
	365	END DO
[825]	366
[921]	367	DO layer = 1, nlay_i-1
	368	DO ji = kideb , kiut
	369	!-- Ice kappa factors
	370	zkappa_i(ji,layer) = 2.0*ztcond_i(ji,layer)/ &
[4333]	371	MAX(epsi10,2.0*zh_i(ji))
[921]	372	END DO
	373	END DO
[825]	374
[921]	375	DO ji = kideb , kiut
[4333]	376	zkappa_i(ji,0) = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji))
	377	zkappa_i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji))
[921]	378	!-- Interface
[4333]	379	zkappa_s(ji,nlay_s) = 2.0rcdsnztcond_i(ji,0)/MAX(epsi10, &
[921]	380	(ztcond_i(ji,0)zh_s(ji) + rcdsnzh_i(ji)))
[4161]	381	zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) &
	382	+ zkappa_i(ji,0)*REAL( 1 - isnow(ji) )
[921]	383	END DO
	384	!
	385	!------------------------------------------------------------------------------\|
	386	! 6) Sea ice specific heat, eta factors \|
	387	!------------------------------------------------------------------------------\|
	388	!
	389	DO layer = 1, nlay_i
	390	DO ji = kideb , kiut
	391	ztitemp(ji,layer) = t_i_b(ji,layer)
	392	zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ &
[4333]	393	MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10)
[921]	394	zeta_i(ji,layer) = rdt_ice / MAX(rhoiczspeche_i(ji,layer)zh_i(ji), &
[4333]	395	epsi10)
[921]	396	END DO
	397	END DO
[825]	398
[921]	399	DO layer = 1, nlay_s
	400	DO ji = kideb , kiut
	401	ztstemp(ji,layer) = t_s_b(ji,layer)
[4333]	402	zeta_s(ji,layer) = rdt_ice / MAX(rhosncpiczh_s(ji),epsi10)
[921]	403	END DO
	404	END DO
	405	!
	406	!------------------------------------------------------------------------------\|
	407	! 7) surface flux computation \|
	408	!------------------------------------------------------------------------------\|
	409	!
	410	DO ji = kideb , kiut
	411	! update of the non solar flux according to the update in T_su
[4688]	412	qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) )
[825]	413
[921]	414	! update incoming flux
	415	zf(ji) = zfsw(ji) & ! net absorbed solar radiation
[4688]	416	+ qns_ice_1d(ji) ! non solar total flux
[921]	417	! (LWup, LWdw, SH, LH)
	418	END DO
[825]	419
[921]	420	!
	421	!------------------------------------------------------------------------------\|
	422	! 8) tridiagonal system terms \|
	423	!------------------------------------------------------------------------------\|
	424	!
	425	!!layer denotes the number of the layer in the snow or in the ice
	426	!!numeq denotes the reference number of the equation in the tridiagonal
	427	!!system, terms of tridiagonal system are indexed as following :
	428	!!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
[825]	429
[921]	430	!!ice interior terms (top equation has the same form as the others)
	431
	432	DO numeq=1,jkmax+2
	433	DO ji = kideb , kiut
	434	ztrid(ji,numeq,1) = 0.
	435	ztrid(ji,numeq,2) = 0.
	436	ztrid(ji,numeq,3) = 0.
	437	zindterm(ji,numeq)= 0.
	438	zindtbis(ji,numeq)= 0.
	439	zdiagbis(ji,numeq)= 0.
	440	ENDDO
	441	ENDDO
	442
	443	DO numeq = nlay_s + 2, nlay_s + nlay_i
	444	DO ji = kideb , kiut
	445	layer = numeq - nlay_s - 1
	446	ztrid(ji,numeq,1) = - zeta_i(ji,layer)*zkappa_i(ji,layer-1)
	447	ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,layer)*(zkappa_i(ji,layer-1) + &
	448	zkappa_i(ji,layer))
	449	ztrid(ji,numeq,3) = - zeta_i(ji,layer)*zkappa_i(ji,layer)
	450	zindterm(ji,numeq) = ztiold(ji,layer) + zeta_i(ji,layer)* &
	451	zradab_i(ji,layer)
	452	END DO
	453	ENDDO
	454
	455	numeq = nlay_s + nlay_i + 1
	456	DO ji = kideb , kiut
[825]	457	!!ice bottom term
	458	ztrid(ji,numeq,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)
	459	ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,nlay_i)( zkappa_i(ji,nlay_i)zg1 &
[921]	460	+ zkappa_i(ji,nlay_i-1) )
[825]	461	ztrid(ji,numeq,3) = 0.0
	462	zindterm(ji,numeq) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i)* &
[921]	463	( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 &
	464	* t_bo_b(ji) )
	465	ENDDO
[825]	466
	467
[921]	468	DO ji = kideb , kiut
[825]	469	IF ( ht_s_b(ji).gt.0.0 ) THEN
[921]	470	!
	471	!------------------------------------------------------------------------------\|
	472	! snow-covered cells \|
	473	!------------------------------------------------------------------------------\|
	474	!
	475	!!snow interior terms (bottom equation has the same form as the others)
	476	DO numeq = 3, nlay_s + 1
	477	layer = numeq - 1
	478	ztrid(ji,numeq,1) = - zeta_s(ji,layer)*zkappa_s(ji,layer-1)
	479	ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,layer)*( zkappa_s(ji,layer-1) + &
	480	zkappa_s(ji,layer) )
	481	ztrid(ji,numeq,3) = - zeta_s(ji,layer)*zkappa_s(ji,layer)
	482	zindterm(ji,numeq) = ztsold(ji,layer) + zeta_s(ji,layer)* &
	483	zradab_s(ji,layer)
	484	END DO
[825]	485
[921]	486	!!case of only one layer in the ice (ice equation is altered)
	487	IF ( nlay_i.eq.1 ) THEN
	488	ztrid(ji,nlay_s+2,3) = 0.0
	489	zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* &
	490	t_bo_b(ji)
	491	ENDIF
[834]	492
[921]	493	IF ( t_su_b(ji) .LT. rtt ) THEN
[825]	494
[921]	495	!------------------------------------------------------------------------------\|
	496	! case 1 : no surface melting - snow present \|
	497	!------------------------------------------------------------------------------\|
	498	zdifcase(ji) = 1.0
	499	numeqmin(ji) = 1
	500	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	501
[921]	502	!!surface equation
	503	ztrid(ji,1,1) = 0.0
	504	ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0)
	505	ztrid(ji,1,3) = zg1s*zkappa_s(ji,0)
	506	zindterm(ji,1) = dzf(ji)*t_su_b(ji) - zf(ji)
[825]	507
[921]	508	!!first layer of snow equation
	509	ztrid(ji,2,1) = - zkappa_s(ji,0)zg1szeta_s(ji,1)
	510	ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)(zkappa_s(ji,1) + zkappa_s(ji,0)zg1s)
	511	ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1)
	512	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)
[825]	513
[921]	514	ELSE
	515	!
	516	!------------------------------------------------------------------------------\|
	517	! case 2 : surface is melting - snow present \|
	518	!------------------------------------------------------------------------------\|
	519	!
	520	zdifcase(ji) = 2.0
	521	numeqmin(ji) = 2
	522	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	523
[921]	524	!!first layer of snow equation
	525	ztrid(ji,2,1) = 0.0
	526	ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + &
	527	zkappa_s(ji,0) * zg1s )
	528	ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1)
	529	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * &
	530	( zradab_s(ji,1) + &
	531	zkappa_s(ji,0) * zg1s * t_su_b(ji) )
	532	ENDIF
	533	ELSE
	534	!
	535	!------------------------------------------------------------------------------\|
	536	! cells without snow \|
	537	!------------------------------------------------------------------------------\|
	538	!
	539	IF (t_su_b(ji) .LT. rtt) THEN
	540	!
	541	!------------------------------------------------------------------------------\|
	542	! case 3 : no surface melting - no snow \|
	543	!------------------------------------------------------------------------------\|
	544	!
	545	zdifcase(ji) = 3.0
	546	numeqmin(ji) = nlay_s + 1
	547	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	548
[921]	549	!!surface equation
	550	ztrid(ji,numeqmin(ji),1) = 0.0
	551	ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1
	552	ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1
	553	zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_b(ji) - zf(ji)
[825]	554
[921]	555	!!first layer of ice equation
	556	ztrid(ji,numeqmin(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
	557	ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) &
	558	+ zkappa_i(ji,0) * zg1 )
	559	ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1)
	560	zindterm(ji,numeqmin(ji)+1)= ztiold(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)
[825]	561
[921]	562	!!case of only one layer in the ice (surface & ice equations are altered)
[825]	563
[921]	564	IF (nlay_i.eq.1) THEN
	565	ztrid(ji,numeqmin(ji),1) = 0.0
	566	ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*2.0
	567	ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*2.0
	568	ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0)2.0zeta_i(ji,1)
	569	ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,0)2.0 + &
	570	zkappa_i(ji,1))
	571	ztrid(ji,numeqmin(ji)+1,3) = 0.0
[825]	572
[921]	573	zindterm(ji,numeqmin(ji)+1) = ztiold(ji,1) + zeta_i(ji,1)* &
	574	( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji) )
	575	ENDIF
[825]	576
[921]	577	ELSE
[825]	578
[921]	579	!
	580	!------------------------------------------------------------------------------\|
	581	! case 4 : surface is melting - no snow \|
	582	!------------------------------------------------------------------------------\|
	583	!
	584	zdifcase(ji) = 4.0
	585	numeqmin(ji) = nlay_s + 2
	586	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	587
[921]	588	!!first layer of ice equation
	589	ztrid(ji,numeqmin(ji),1) = 0.0
	590	ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,1) + zkappa_i(ji,0) &
	591	zg1)
	592	ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)
	593	zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &
	594	zkappa_i(ji,0) * zg1 * t_su_b(ji) )
[825]	595
[921]	596	!!case of only one layer in the ice (surface & ice equations are altered)
	597	IF (nlay_i.eq.1) THEN
	598	ztrid(ji,numeqmin(ji),1) = 0.0
	599	ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,0)2.0 + &
	600	zkappa_i(ji,1))
	601	ztrid(ji,numeqmin(ji),3) = 0.0
	602	zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)* &
	603	(zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji)) &
	604	+ t_su_b(ji)zeta_i(ji,1)zkappa_i(ji,0)*2.0
	605	ENDIF
[825]	606
[921]	607	ENDIF
	608	ENDIF
[825]	609
[921]	610	END DO
[825]	611
[921]	612	!
	613	!------------------------------------------------------------------------------\|
	614	! 9) tridiagonal system solving \|
	615	!------------------------------------------------------------------------------\|
	616	!
[825]	617
[921]	618	! Solve the tridiagonal system with Gauss elimination method.
	619	! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
	620	! McGraw-Hill 1984.
[825]	621
[921]	622	maxnumeqmax = 0
	623	minnumeqmin = jkmax+4
[825]	624
[921]	625	DO ji = kideb , kiut
	626	zindtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji))
	627	zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2)
	628	minnumeqmin = MIN(numeqmin(ji),minnumeqmin)
	629	maxnumeqmax = MAX(numeqmax(ji),maxnumeqmax)
	630	END DO
	631
	632	DO layer = minnumeqmin+1, maxnumeqmax
	633	DO ji = kideb , kiut
	634	numeq = min(max(numeqmin(ji)+1,layer),numeqmax(ji))
	635	zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* &
	636	ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1)
	637	zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1)* &
	638	zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1)
	639	END DO
	640	END DO
	641
	642	DO ji = kideb , kiut
	643	! ice temperatures
	644	t_i_b(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))
	645	END DO
	646
	647	DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1
	648	DO ji = kideb , kiut
	649	layer = numeq - nlay_s - 1
	650	t_i_b(ji,layer) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &
	651	t_i_b(ji,layer+1))/zdiagbis(ji,numeq)
	652	END DO
	653	END DO
	654
	655	DO ji = kideb , kiut
[825]	656	! snow temperatures
[4205]	657	IF (ht_s_b(ji).GT.0._wp) &
[921]	658	t_s_b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
	659	* t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) &
[4161]	660	* MAX(0.0,SIGN(1.0,ht_s_b(ji)))
[825]	661
	662	! surface temperature
[4161]	663	isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) ) ) )
[2715]	664	ztsuoldit(ji) = t_su_b(ji)
[4161]	665	IF( t_su_b(ji) < ztfs(ji) ) &
	666	t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1) &
	667	& + REAL( 1 - isnow(ji) )*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))
[921]	668	END DO
	669	!
	670	!--------------------------------------------------------------------------
	671	! 10) Has the scheme converged ?, end of the iterative procedure \|
	672	!--------------------------------------------------------------------------
	673	!
	674	! check that nowhere it has started to melt
	675	! zerrit(ji) is a measure of error, it has to be under maxer_i_thd
	676	DO ji = kideb , kiut
[2715]	677	t_su_b(ji) = MAX( MIN( t_su_b(ji) , ztfs(ji) ) , 190._wp )
	678	zerrit(ji) = ABS( t_su_b(ji) - ztsuoldit(ji) )
[921]	679	END DO
[825]	680
[921]	681	DO layer = 1, nlay_s
	682	DO ji = kideb , kiut
[2715]	683	t_s_b(ji,layer) = MAX( MIN( t_s_b(ji,layer), rtt ), 190._wp )
	684	zerrit(ji) = MAX(zerrit(ji),ABS(t_s_b(ji,layer) - ztstemp(ji,layer)))
[921]	685	END DO
	686	END DO
[825]	687
[921]	688	DO layer = 1, nlay_i
	689	DO ji = kideb , kiut
[2715]	690	ztmelt_i = -tmut * s_i_b(ji,layer) + rtt
[4205]	691	t_i_b(ji,layer) = MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp)
[2715]	692	zerrit(ji) = MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer)))
[921]	693	END DO
	694	END DO
[825]	695
[921]	696	! Compute spatial maximum over all errors
[2715]	697	! note that this could be optimized substantially by iterating only the non-converging points
	698	zerritmax = 0._wp
	699	DO ji = kideb, kiut
	700	zerritmax = MAX( zerritmax, zerrit(ji) )
[921]	701	END DO
[2715]	702	IF( lk_mpp ) CALL mpp_max( zerritmax, kcom=ncomm_ice )
[825]	703
	704	END DO ! End of the do while iterative procedure
	705
[4333]	706	IF( ln_nicep .AND. lwp ) THEN
[1055]	707	WRITE(numout,*) ' zerritmax : ', zerritmax
	708	WRITE(numout,*) ' nconv : ', nconv
	709	ENDIF
[825]	710
[921]	711	!
[2715]	712	!-------------------------------------------------------------------------!
	713	! 11) Fluxes at the interfaces !
	714	!-------------------------------------------------------------------------!
[921]	715	DO ji = kideb, kiut
[3808]	716	! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)
[4688]	717	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]	718	! ! surface ice conduction flux
[4161]	719	isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) ) )
	720	fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji)) &
	721	& - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_b(ji,1) - t_su_b(ji))
[2715]	722	! ! bottom ice conduction flux
	723	fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
[921]	724	END DO
[825]	725
[4688]	726	!-----------------------------------------
	727	! Heat flux used to warm/cool ice in W.m-2
	728	!-----------------------------------------
	729	DO ji = kideb, kiut
	730	IF( t_su_b(ji) < rtt ) THEN ! case T_su < 0degC
[4765]	731	hfx_dif_1d(ji) = hfx_dif_1d(ji) + &
	732	& ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji)
[4688]	733	ELSE ! case T_su = 0degC
[4765]	734	hfx_dif_1d(ji) = hfx_dif_1d(ji) + &
	735	& ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji)
[4688]	736	ENDIF
	737	END DO
	738
	739	! --- computes sea ice energy of melting compulsory for limthd_dh --- !
	740	CALL lim_thd_enmelt( kideb, kiut )
	741
	742	! --- diag error on heat diffusion - PART 2 --- !
	743	DO ji = kideb, kiut
	744	zdq(ji) = - zq_ini(ji) + ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + &
	745	& SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )
	746	zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice )
	747	hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err * a_i_b(ji)
	748	! --- correction of qns_ice and surface conduction flux --- !
	749	qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err
	750	fc_su (ji) = fc_su (ji) - zhfx_err
	751	! --- Heat flux at the ice surface in W.m-2 --- !
	752	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
	753	hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) )
	754	END DO
	755
	756	!
	757	CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow )
	758	CALL wrk_dealloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw )
	759	CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu )
[4765]	760	CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, &
	761	& ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 )
[4688]	762	CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart = 0 )
	763	CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis )
	764	CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid )
	765	CALL wrk_dealloc( jpij, zdq, zq_ini )
	766
	767	END SUBROUTINE lim_thd_dif
	768
	769	SUBROUTINE lim_thd_enmelt( kideb, kiut )
	770	!!-----------------------------------------------------------------------
	771	!! * ROUTINE lim_thd_enmelt *
	772	!!
	773	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature
	774	!!
	775	!! ** Method : Formula (Bitz and Lipscomb, 1999)
	776	!!-------------------------------------------------------------------
	777	INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop
	778	!
	779	INTEGER :: ji, jk ! dummy loop indices
	780	REAL(wp) :: ztmelts, zindb ! local scalar
	781	!!-------------------------------------------------------------------
	782	!
	783	DO jk = 1, nlay_i ! Sea ice energy of melting
[921]	784	DO ji = kideb, kiut
[4688]	785	ztmelts = - tmut * s_i_b(ji,jk) + rtt
	786	zindb = MAX( 0._wp , SIGN( 1._wp , -(t_i_b(ji,jk) - rtt) - epsi10 ) )
	787	q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) ) &
	788	& + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) ) &
	789	& - rcp * ( ztmelts-rtt ) )
[921]	790	END DO
[4688]	791	END DO
	792	DO jk = 1, nlay_s ! Snow energy of melting
	793	DO ji = kideb, kiut
	794	q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
[921]	795	END DO
[4688]	796	END DO
[2715]	797	!
[4688]	798	END SUBROUTINE lim_thd_enmelt
[825]	799
	800	#else
[2715]	801	!!----------------------------------------------------------------------
	802	!! Dummy Module No LIM-3 sea-ice model
	803	!!----------------------------------------------------------------------
[825]	804	CONTAINS
	805	SUBROUTINE lim_thd_dif ! Empty routine
	806	END SUBROUTINE lim_thd_dif
	807	#endif
[2528]	808	!!======================================================================
[921]	809	END MODULE limthd_dif

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

Download in other formats: