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

Last change on this file since 5051 was 5051, checked in by clem, 9 years ago
LIM3 initialization is now called at the same time as other sbc fields
Property svn:keywords set to `Id`
File size: 42.7 KB

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

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

Download in other formats: