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

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

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

Download in other formats: