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_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 5357

Last change on this file since 5357 was 5357, checked in by clem, 9 years ago
LIM3: change the interface between the ice and atm for both coupled and forced modes. Some work still needs to be done to deal with sublimation in coupled mode.
Property svn:keywords set to `Id`
File size: 42.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 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)
[921]	26
[825]	27	IMPLICIT NONE
	28	PRIVATE
	29
[2528]	30	PUBLIC lim_thd_dif ! called by lim_thd
[825]	31
	32	!!----------------------------------------------------------------------
[4161]	33	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	34	!! $Id$
[2591]	35	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	36	!!----------------------------------------------------------------------
	37	CONTAINS
	38
[4688]	39	SUBROUTINE lim_thd_dif( kideb , kiut )
[921]	40	!!------------------------------------------------------------------
	41	!! * ROUTINE lim_thd_dif *
	42	!! ** Purpose :
	43	!! This routine determines the time evolution of snow and sea-ice
	44	!! temperature profiles.
	45	!! ** Method :
	46	!! This is done by solving the heat equation diffusion with
	47	!! a Neumann boundary condition at the surface and a Dirichlet one
	48	!! at the bottom. Solar radiation is partially absorbed into the ice.
	49	!! The specific heat and thermal conductivities depend on ice salinity
	50	!! and temperature to take into account brine pocket melting. The
	51	!! numerical
	52	!! scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid
	53	!! in the ice and snow system.
	54	!!
	55	!! The successive steps of this routine are
	56	!! 1. Thermal conductivity at the interfaces of the ice layers
	57	!! 2. Internal absorbed radiation
	58	!! 3. Scale factors due to non-uniform grid
	59	!! 4. Kappa factors
	60	!! Then iterative procedure begins
	61	!! 5. specific heat in the ice
	62	!! 6. eta factors
	63	!! 7. surface flux computation
	64	!! 8. tridiagonal system terms
	65	!! 9. solving the tridiagonal system with Gauss elimination
	66	!! Iterative procedure ends according to a criterion on evolution
	67	!! of temperature
	68	!!
	69	!! ** Arguments :
	70	!! kideb , kiut : Starting and ending points on which the
	71	!! the computation is applied
	72	!!
	73	!! ** Inputs / Ouputs : (global commons)
[4872]	74	!! surface temperature : t_su_1d
	75	!! ice/snow temperatures : t_i_1d, t_s_1d
	76	!! ice salinities : s_i_1d
[921]	77	!! number of layers in the ice/snow: nlay_i, nlay_s
	78	!! profile of the ice/snow layers : z_i, z_s
[4872]	79	!! total ice/snow thickness : ht_i_1d, ht_s_1d
[921]	80	!!
	81	!! ** External :
	82	!!
	83	!! ** References :
	84	!!
	85	!! ** History :
	86	!! (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium
	87	!! (06-2005) Martin Vancoppenolle, 3d version
	88	!! (11-2006) Vectorized by Xavier Fettweis (UCL-ASTR)
	89	!! (04-2007) Energy conservation tested by M. Vancoppenolle
	90	!!------------------------------------------------------------------
[4688]	91	INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied
[825]	92
[921]	93	!! * Local variables
[3294]	94	INTEGER :: ji ! spatial loop index
	95	INTEGER :: ii, ij ! temporary dummy loop index
	96	INTEGER :: numeq ! current reference number of equation
[4870]	97	INTEGER :: jk ! vertical dummy loop index
[3294]	98	INTEGER :: nconv ! number of iterations in iterative procedure
	99	INTEGER :: minnumeqmin, maxnumeqmax
[5123]	100
[4688]	101	INTEGER, POINTER, DIMENSION(:) :: numeqmin ! reference number of top equation
	102	INTEGER, POINTER, DIMENSION(:) :: numeqmax ! reference number of bottom equation
[5123]	103
[3294]	104	REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system
	105	REAL(wp) :: zg1 = 2._wp !
	106	REAL(wp) :: zgamma = 18009._wp ! for specific heat
	107	REAL(wp) :: zbeta = 0.117_wp ! for thermal conductivity (could be 0.13)
[4990]	108	REAL(wp) :: zraext_s = 10._wp ! extinction coefficient of radiation in the snow
[3294]	109	REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity
[4688]	110	REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered as 0°C
[2715]	111	REAL(wp) :: ztmelt_i ! ice melting temperature
	112	REAL(wp) :: zerritmax ! current maximal error on temperature
[5123]	113	REAL(wp) :: zhsu
	114
	115	REAL(wp), POINTER, DIMENSION(:) :: isnow ! switch for presence (1) or absence (0) of snow
	116	REAL(wp), POINTER, DIMENSION(:) :: ztsub ! old surface temperature (before the iterative procedure )
	117	REAL(wp), POINTER, DIMENSION(:) :: ztsubit ! 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
[5146]	121	REAL(wp), POINTER, DIMENSION(:) :: zqns_ice_b ! solar radiation absorbed at the surface
[5123]	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
[5123]	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	!
[5123]	168	CALL wrk_alloc( jpij, numeqmin, numeqmax )
	169	CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw )
[5146]	170	CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe )
[5167]	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 )
	173	CALL wrk_alloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis )
	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
[5123]	181	zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + &
	182	& SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
[4688]	183	END DO
	184
[921]	185	!------------------------------------------------------------------------------!
	186	! 1) Initialization !
	187	!------------------------------------------------------------------------------!
	188	DO ji = kideb , kiut
[5123]	189	isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) ) ! is there snow or not
[921]	190	! layer thickness
[5123]	191	zh_i(ji) = ht_i_1d(ji) * r1_nlay_i
	192	zh_s(ji) = ht_s_1d(ji) * r1_nlay_s
[921]	193	END DO
[825]	194
[921]	195	!--------------------
	196	! Ice / snow layers
	197	!--------------------
[825]	198
[2715]	199	z_s(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st snow layer
	200	z_i(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st ice layer
[825]	201
[4870]	202	DO jk = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer
[921]	203	DO ji = kideb , kiut
[5123]	204	z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s
[921]	205	END DO
	206	END DO
[825]	207
[4870]	208	DO jk = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer
[921]	209	DO ji = kideb , kiut
[5123]	210	z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i
[921]	211	END DO
	212	END DO
	213	!
	214	!------------------------------------------------------------------------------\|
[5123]	215	! 2) Radiation \|
[921]	216	!------------------------------------------------------------------------------\|
	217	!
	218	!-------------------
	219	! Computation of i0
	220	!-------------------
	221	! i0 describes the fraction of solar radiation which does not contribute
	222	! to the surface energy budget but rather penetrates inside the ice.
	223	! We assume that no radiation is transmitted through the snow
	224	! If there is no no snow
	225	! zfsw = (1-i0).qsr_ice is absorbed at the surface
	226	! zftrice = io.qsr_ice is below the surface
[4688]	227	! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
[5123]	228	! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
	229	zhsu = 0.1_wp ! threshold for the computation of i0
[921]	230	DO ji = kideb , kiut
	231	! switches
[5123]	232	isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )
[921]	233	! hs > 0, isnow = 1
[5123]	234	zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) )
[825]	235
[5123]	236	i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
[921]	237	END DO
[825]	238
[921]	239	!-------------------------------------------------------
	240	! Solar radiation absorbed / transmitted at the surface
	241	! Derivative of the non solar flux
	242	!-------------------------------------------------------
	243	DO ji = kideb , kiut
[5146]	244	zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface
	245	zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer
	246	dzf (ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux
	247	zqns_ice_b(ji) = qns_ice_1d(ji) ! store previous qns_ice_1d value
[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
[5123]	268	zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - 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
[5123]	274	zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_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
[5123]	292	t_su_1d(ji) = MIN( t_su_1d(ji), rt0 - 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
[5123]	311	DO WHILE ( zerritmax > rn_terr_dif .AND. nconv < nn_conv_dif )
[921]	312	!
[2715]	313	nconv = nconv + 1
	314	!
[921]	315	!------------------------------------------------------------------------------\|
	316	! 4) Sea ice thermal conductivity \|
	317	!------------------------------------------------------------------------------\|
	318	!
[5123]	319	IF( nn_ice_thcon == 0 ) THEN ! Untersteiner (1964) formula
[921]	320	DO ji = kideb , kiut
[5123]	321	ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
	322	ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
[921]	323	END DO
[4870]	324	DO jk = 1, nlay_i-1
[921]	325	DO ji = kideb , kiut
[5123]	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 * rt0)
	328	ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
[921]	329	END DO
	330	END DO
	331	ENDIF
[825]	332
[5123]	333	IF( nn_ice_thcon == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
[921]	334	DO ji = kideb , kiut
[5123]	335	ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) &
	336	& - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
[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
[5123]	341	ztcond_i(ji,jk) = rcdic + &
	342	& 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) &
	343	& / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 ) &
	344	& - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 )
[4870]	345	ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
[2715]	346	END DO
	347	END DO
	348	DO ji = kideb , kiut
[5123]	349	ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) &
	350	& - 0.011_wp * ( t_bo_1d(ji) - rt0 )
[2715]	351	ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
	352	END DO
[921]	353	ENDIF
[5123]	354
[921]	355	!
	356	!------------------------------------------------------------------------------\|
[5123]	357	! 5) G(he) - enhancement of thermal conductivity in mono-category case \|
[921]	358	!------------------------------------------------------------------------------\|
	359	!
[5123]	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
	364	zghe(:) = 1._wp
	365
	366	SELECT CASE ( nn_monocat )
	367
	368	CASE (1,3) ! LIM3
	369
	370	zepsilon = 0.1_wp
	371	zh_thres = EXP( 1._wp ) * zepsilon * 0.5_wp
	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)
	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._wp * zhe / zepsilon ) )
	387
	388	! Impose G(he) < 2.
	389	zghe(ji) = MIN( zghe(ji), 2._wp )
	390
	391	END DO
	392
	393	END SELECT
	394
	395	!
	396	!------------------------------------------------------------------------------\|
	397	! 6) kappa factors \|
	398	!------------------------------------------------------------------------------\|
	399	!
	400	!--- Snow
[921]	401	DO ji = kideb, kiut
[5123]	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
[5123]	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
[5123]	413	!--- Ice
[4870]	414	DO jk = 1, nlay_i-1
[921]	415	DO ji = kideb , kiut
[5123]	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
[5123]	420	!--- Snow-ice interface
[921]	421	DO ji = kideb , kiut
[5123]	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) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
[921]	428	END DO
[5123]	429
[921]	430	!
	431	!------------------------------------------------------------------------------\|
[5123]	432	! 7) 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)
[5123]	438	zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 )
	439	zeta_i(ji,jk) = rdt_ice / MAX( rhoic * zspeche_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)
[5123]	446	zeta_s(ji,jk) = rdt_ice / MAX( rhosn * cpic * zh_s(ji), epsi10 )
[921]	447	END DO
	448	END DO
[5123]	449
[921]	450	!
	451	!------------------------------------------------------------------------------\|
[5123]	452	! 8) surface flux computation \|
[921]	453	!------------------------------------------------------------------------------\|
	454	!
[5146]	455	IF ( ln_it_qnsice ) THEN
[4990]	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
[5123]	465	zf(ji) = zfsw(ji) & ! net absorbed solar radiation
	466	& + qns_ice_1d(ji) ! non solar total flux (LWup, LWdw, SH, LH)
[921]	467	END DO
[825]	468
[921]	469	!
	470	!------------------------------------------------------------------------------\|
[5123]	471	! 9) tridiagonal system terms \|
[921]	472	!------------------------------------------------------------------------------\|
	473	!
	474	!!layer denotes the number of the layer in the snow or in the ice
	475	!!numeq denotes the reference number of the equation in the tridiagonal
	476	!!system, terms of tridiagonal system are indexed as following :
	477	!!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
[825]	478
[921]	479	!!ice interior terms (top equation has the same form as the others)
	480
[4873]	481	DO numeq=1,nlay_i+3
[921]	482	DO ji = kideb , kiut
	483	ztrid(ji,numeq,1) = 0.
	484	ztrid(ji,numeq,2) = 0.
	485	ztrid(ji,numeq,3) = 0.
[5123]	486	zindterm(ji,numeq)= 0.
	487	zindtbis(ji,numeq)= 0.
[921]	488	zdiagbis(ji,numeq)= 0.
	489	ENDDO
	490	ENDDO
	491
	492	DO numeq = nlay_s + 2, nlay_s + nlay_i
	493	DO ji = kideb , kiut
[5123]	494	jk = numeq - nlay_s - 1
	495	ztrid(ji,numeq,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
	496	ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
	497	ztrid(ji,numeq,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk)
	498	zindterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
[921]	499	END DO
	500	ENDDO
	501
	502	numeq = nlay_s + nlay_i + 1
	503	DO ji = kideb , kiut
[825]	504	!!ice bottom term
	505	ztrid(ji,numeq,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)
[5123]	506	ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
[825]	507	ztrid(ji,numeq,3) = 0.0
[5123]	508	zindterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i) * &
	509	& ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) )
[921]	510	ENDDO
[825]	511
	512
[921]	513	DO ji = kideb , kiut
[5123]	514	IF ( ht_s_1d(ji) > 0.0 ) THEN
[921]	515	!
	516	!------------------------------------------------------------------------------\|
	517	! snow-covered cells \|
	518	!------------------------------------------------------------------------------\|
	519	!
	520	!!snow interior terms (bottom equation has the same form as the others)
	521	DO numeq = 3, nlay_s + 1
[5123]	522	jk = numeq - 1
	523	ztrid(ji,numeq,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
	524	ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
[4870]	525	ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk)
[5123]	526	zindterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
[921]	527	END DO
[825]	528
[921]	529	!!case of only one layer in the ice (ice equation is altered)
	530	IF ( nlay_i.eq.1 ) THEN
	531	ztrid(ji,nlay_s+2,3) = 0.0
[5123]	532	zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)
[921]	533	ENDIF
[834]	534
[5123]	535	IF ( t_su_1d(ji) < rt0 ) THEN
[825]	536
[921]	537	!------------------------------------------------------------------------------\|
	538	! case 1 : no surface melting - snow present \|
	539	!------------------------------------------------------------------------------\|
	540	zdifcase(ji) = 1.0
	541	numeqmin(ji) = 1
	542	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	543
[921]	544	!!surface equation
[5123]	545	ztrid(ji,1,1) = 0.0
	546	ztrid(ji,1,2) = dzf(ji) - zg1s * zkappa_s(ji,0)
	547	ztrid(ji,1,3) = zg1s * zkappa_s(ji,0)
	548	zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji)
[825]	549
[921]	550	!!first layer of snow equation
[5123]	551	ztrid(ji,2,1) = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
	552	ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
[921]	553	ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1)
[5123]	554	zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
[825]	555
[921]	556	ELSE
	557	!
	558	!------------------------------------------------------------------------------\|
	559	! case 2 : surface is melting - snow present \|
	560	!------------------------------------------------------------------------------\|
	561	!
	562	zdifcase(ji) = 2.0
	563	numeqmin(ji) = 2
	564	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	565
[921]	566	!!first layer of snow equation
	567	ztrid(ji,2,1) = 0.0
[5123]	568	ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
[921]	569	ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1)
[5123]	570	zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * &
	571	& ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )
[921]	572	ENDIF
	573	ELSE
	574	!
	575	!------------------------------------------------------------------------------\|
	576	! cells without snow \|
	577	!------------------------------------------------------------------------------\|
	578	!
[5123]	579	IF ( t_su_1d(ji) < rt0 ) THEN
[921]	580	!
	581	!------------------------------------------------------------------------------\|
	582	! case 3 : no surface melting - no snow \|
	583	!------------------------------------------------------------------------------\|
	584	!
	585	zdifcase(ji) = 3.0
	586	numeqmin(ji) = nlay_s + 1
	587	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	588
[921]	589	!!surface equation
	590	ztrid(ji,numeqmin(ji),1) = 0.0
	591	ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1
	592	ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1
[5123]	593	zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji)
[825]	594
[921]	595	!!first layer of ice equation
	596	ztrid(ji,numeqmin(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
[5123]	597	ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
	598	ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1)
	599	zindterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)
[825]	600
[921]	601	!!case of only one layer in the ice (surface & ice equations are altered)
[825]	602
[5123]	603	IF ( nlay_i == 1 ) THEN
[921]	604	ztrid(ji,numeqmin(ji),1) = 0.0
[5123]	605	ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0) * 2.0
	606	ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0) * 2.0
	607	ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
	608	ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
[921]	609	ztrid(ji,numeqmin(ji)+1,3) = 0.0
[825]	610
[5123]	611	zindterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1) * &
	612	& ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
[921]	613	ENDIF
[825]	614
[921]	615	ELSE
[825]	616
[921]	617	!
	618	!------------------------------------------------------------------------------\|
	619	! case 4 : surface is melting - no snow \|
	620	!------------------------------------------------------------------------------\|
	621	!
	622	zdifcase(ji) = 4.0
	623	numeqmin(ji) = nlay_s + 2
	624	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	625
[921]	626	!!first layer of ice equation
	627	ztrid(ji,numeqmin(ji),1) = 0.0
[5123]	628	ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
[921]	629	ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)
[5123]	630	zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1) * &
	631	& ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) )
[825]	632
[921]	633	!!case of only one layer in the ice (surface & ice equations are altered)
[5123]	634	IF ( nlay_i == 1 ) THEN
[921]	635	ztrid(ji,numeqmin(ji),1) = 0.0
[5123]	636	ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
[921]	637	ztrid(ji,numeqmin(ji),3) = 0.0
[5123]	638	zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &
	639	& + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
[921]	640	ENDIF
[825]	641
[921]	642	ENDIF
	643	ENDIF
[825]	644
[921]	645	END DO
[825]	646
[921]	647	!
	648	!------------------------------------------------------------------------------\|
[5123]	649	! 10) tridiagonal system solving \|
[921]	650	!------------------------------------------------------------------------------\|
	651	!
[825]	652
[921]	653	! Solve the tridiagonal system with Gauss elimination method.
	654	! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
	655	! McGraw-Hill 1984.
[825]	656
[921]	657	maxnumeqmax = 0
[4873]	658	minnumeqmin = nlay_i+5
[825]	659
[921]	660	DO ji = kideb , kiut
[5123]	661	zindtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji))
[921]	662	zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2)
	663	minnumeqmin = MIN(numeqmin(ji),minnumeqmin)
	664	maxnumeqmax = MAX(numeqmax(ji),maxnumeqmax)
	665	END DO
	666
[4870]	667	DO jk = minnumeqmin+1, maxnumeqmax
[921]	668	DO ji = kideb , kiut
[4870]	669	numeq = min(max(numeqmin(ji)+1,jk),numeqmax(ji))
[5123]	670	zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3) / zdiagbis(ji,numeq-1)
	671	zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1)
[921]	672	END DO
	673	END DO
	674
	675	DO ji = kideb , kiut
	676	! ice temperatures
[5123]	677	t_i_1d(ji,nlay_i) = zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji))
[921]	678	END DO
	679
[5180]	680	DO numeq = nlay_i + nlay_s, nlay_s + 2, -1
[921]	681	DO ji = kideb , kiut
[4870]	682	jk = numeq - nlay_s - 1
[5123]	683	t_i_1d(ji,jk) = ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq)
[921]	684	END DO
	685	END DO
	686
	687	DO ji = kideb , kiut
[825]	688	! snow temperatures
[5123]	689	IF (ht_s_1d(ji) > 0._wp) &
	690	t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) &
	691	& / zdiagbis(ji,nlay_s+1) * MAX( 0.0, SIGN( 1.0, ht_s_1d(ji) ) )
[825]	692
	693	! surface temperature
[5123]	694	isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) )
[4872]	695	ztsubit(ji) = t_su_1d(ji)
[5123]	696	IF( t_su_1d(ji) < rt0 ) &
	697	t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) * &
	698	& ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))
[921]	699	END DO
	700	!
	701	!--------------------------------------------------------------------------
[5123]	702	! 11) Has the scheme converged ?, end of the iterative procedure \|
[921]	703	!--------------------------------------------------------------------------
	704	!
	705	! check that nowhere it has started to melt
[5123]	706	! zerrit(ji) is a measure of error, it has to be under terr_dif
[921]	707	DO ji = kideb , kiut
[5123]	708	t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , 190._wp )
	709	zerrit(ji) = ABS( t_su_1d(ji) - ztsubit(ji) )
[921]	710	END DO
[825]	711
[4870]	712	DO jk = 1, nlay_s
[921]	713	DO ji = kideb , kiut
[5123]	714	t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), 190._wp )
	715	zerrit(ji) = MAX( zerrit(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) )
[921]	716	END DO
	717	END DO
[825]	718
[4870]	719	DO jk = 1, nlay_i
[921]	720	DO ji = kideb , kiut
[5123]	721	ztmelt_i = -tmut * s_i_1d(ji,jk) + rt0
	722	t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp )
	723	zerrit(ji) = MAX( zerrit(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) )
[921]	724	END DO
	725	END DO
[825]	726
[921]	727	! Compute spatial maximum over all errors
[2715]	728	! note that this could be optimized substantially by iterating only the non-converging points
	729	zerritmax = 0._wp
	730	DO ji = kideb, kiut
	731	zerritmax = MAX( zerritmax, zerrit(ji) )
[921]	732	END DO
[2715]	733	IF( lk_mpp ) CALL mpp_max( zerritmax, kcom=ncomm_ice )
[825]	734
	735	END DO ! End of the do while iterative procedure
	736
[5128]	737	IF( ln_icectl .AND. lwp ) THEN
[1055]	738	WRITE(numout,*) ' zerritmax : ', zerritmax
	739	WRITE(numout,*) ' nconv : ', nconv
	740	ENDIF
[825]	741
[921]	742	!
[2715]	743	!-------------------------------------------------------------------------!
[5123]	744	! 12) Fluxes at the interfaces !
[2715]	745	!-------------------------------------------------------------------------!
[921]	746	DO ji = kideb, kiut
[2715]	747	! ! surface ice conduction flux
[5123]	748	isnow(ji) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) )
	749	fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) &
	750	& - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji))
[2715]	751	! ! bottom ice conduction flux
[4872]	752	fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
[921]	753	END DO
[825]	754
[4688]	755	! --- computes sea ice energy of melting compulsory for limthd_dh --- !
	756	CALL lim_thd_enmelt( kideb, kiut )
	757
[5146]	758	! --- diagnose the change in non-solar flux due to surface temperature change --- !
	759	IF ( ln_it_qnsice ) hfx_err_dif_1d(:) = hfx_err_dif_1d(:) - ( qns_ice_1d(:) - zqns_ice_b(:) ) * a_i_1d(:)
[5123]	760
[4990]	761	! --- diag conservation imbalance on heat diffusion - PART 2 --- !
[4688]	762	DO ji = kideb, kiut
[5123]	763	zdq(ji) = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + &
	764	& SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
	765	IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC
	766	zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice
	767	ELSE ! case T_su = 0degC
	768	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
	769	ENDIF
[4990]	770	hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
[5167]	771
	772	! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
	773	hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
[4990]	774	END DO
	775
[5123]	776	!-----------------------------------------
	777	! Heat flux used to warm/cool ice in W.m-2
	778	!-----------------------------------------
[4990]	779	DO ji = kideb, kiut
[5123]	780	IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC
	781	hfx_dif_1d(ji) = hfx_dif_1d(ji) + &
	782	& ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
	783	ELSE ! case T_su = 0degC
	784	hfx_dif_1d(ji) = hfx_dif_1d(ji) + &
	785	& ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
	786	ENDIF
[5167]	787	! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
	788	hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji)
	789	END DO
[4688]	790	!
[5123]	791	CALL wrk_dealloc( jpij, numeqmin, numeqmax )
	792	CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw )
	793	CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe )
[5167]	794	CALL wrk_dealloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 )
	795	CALL wrk_dealloc( jpij,nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 )
	796	CALL wrk_dealloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis )
	797	CALL wrk_dealloc( jpij,nlay_i+3,3, ztrid )
[4990]	798	CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err )
[4688]	799
	800	END SUBROUTINE lim_thd_dif
	801
	802	SUBROUTINE lim_thd_enmelt( kideb, kiut )
	803	!!-----------------------------------------------------------------------
	804	!! * ROUTINE lim_thd_enmelt *
	805	!!
	806	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature
	807	!!
	808	!! ** Method : Formula (Bitz and Lipscomb, 1999)
	809	!!-------------------------------------------------------------------
	810	INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop
	811	!
	812	INTEGER :: ji, jk ! dummy loop indices
[4990]	813	REAL(wp) :: ztmelts ! local scalar
[4688]	814	!!-------------------------------------------------------------------
	815	!
	816	DO jk = 1, nlay_i ! Sea ice energy of melting
[921]	817	DO ji = kideb, kiut
[5202]	818	ztmelts = - tmut * s_i_1d(ji,jk) + rt0
	819	t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point
	820	! (sometimes dif scheme produces abnormally high temperatures)
	821	q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) &
	822	& + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) ) &
	823	& - rcp * ( ztmelts-rt0 ) )
[921]	824	END DO
[4688]	825	END DO
	826	DO jk = 1, nlay_s ! Snow energy of melting
	827	DO ji = kideb, kiut
[5123]	828	q_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
[921]	829	END DO
[4688]	830	END DO
[2715]	831	!
[4688]	832	END SUBROUTINE lim_thd_enmelt
[825]	833
	834	#else
[2715]	835	!!----------------------------------------------------------------------
	836	!! Dummy Module No LIM-3 sea-ice model
	837	!!----------------------------------------------------------------------
[825]	838	CONTAINS
	839	SUBROUTINE lim_thd_dif ! Empty routine
	840	END SUBROUTINE lim_thd_dif
	841	#endif
[2528]	842	!!======================================================================
[921]	843	END MODULE limthd_dif

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

Download in other formats: