Context Navigation

icethd_zdf_bl99.F90 @ 13046

Last change on this file since 13046 was 13046, checked in by dancopsey, 4 years ago
Add snow_frac which is fraction of snow for coupling purposes.
File size: 44.6 KB

Rev	Line
[8984]	1	MODULE icethd_zdf_BL99
	2	!!======================================================================
	3	!! * MODULE icethd_zdf_BL99 *
	4	!! sea-ice: vertical heat diffusion in sea ice (computation of temperatures)
	5	!!======================================================================
[9656]	6	!! History : ! 2003-02 (M. Vancoppenolle) original 1D code
[9604]	7	!! ! 2005-06 (M. Vancoppenolle) 3d version
	8	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
[8984]	9	!!----------------------------------------------------------------------
[9570]	10	#if defined key_si3
[8984]	11	!!----------------------------------------------------------------------
[9570]	12	!! 'key_si3' SI3 sea-ice model
[8984]	13	!!----------------------------------------------------------------------
	14	!! ice_thd_zdf_BL99 : vertical diffusion computation
	15	!!----------------------------------------------------------------------
	16	USE dom_oce ! ocean space and time domain
	17	USE phycst ! physical constants (ocean directory)
	18	USE ice ! sea-ice: variables
	19	USE ice1D ! sea-ice: thermodynamics variables
	20	USE icevar ! sea-ice: operations
	21	!
	22	USE in_out_manager ! I/O manager
	23	USE lib_mpp ! MPP library
	24	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
	25
	26	IMPLICIT NONE
	27	PRIVATE
	28
	29	PUBLIC ice_thd_zdf_BL99 ! called by icethd_zdf
	30
	31	!!----------------------------------------------------------------------
[9598]	32	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	33	!! $Id$
[10068]	34	!! Software governed by the CeCILL license (see ./LICENSE)
[8984]	35	!!----------------------------------------------------------------------
	36	CONTAINS
	37
[10534]	38	SUBROUTINE ice_thd_zdf_BL99( k_cnd )
[8984]	39	!!-------------------------------------------------------------------
	40	!! * ROUTINE ice_thd_zdf_BL99 *
	41	!!
	42	!! ** Purpose : computes the time evolution of snow and sea-ice temperature
	43	!! profiles, using the original Bitz and Lipscomb (1999) algorithm
	44	!!
	45	!! ** Method : solves the heat equation diffusion with a Neumann boundary
	46	!! condition at the surface and a Dirichlet one at the bottom.
	47	!! Solar radiation is partially absorbed into the ice.
	48	!! The specific heat and thermal conductivities depend on ice
	49	!! salinity and temperature to take into account brine pocket
	50	!! melting. The numerical scheme is an iterative Crank-Nicolson
	51	!! on a non-uniform multilayer grid in the ice and snow system.
	52	!!
	53	!! The successive steps of this routine are
	54	!! 1. initialization of ice-snow layers thicknesses
	55	!! 2. Internal absorbed and transmitted radiation
	56	!! Then iterative procedure begins
	57	!! 3. Thermal conductivity
	58	!! 4. Kappa factors
	59	!! 5. specific heat in the ice
	60	!! 6. eta factors
	61	!! 7. surface flux computation
	62	!! 8. tridiagonal system terms
	63	!! 9. solving the tridiagonal system with Gauss elimination
	64	!! Iterative procedure ends according to a criterion on evolution
	65	!! of temperature
	66	!! 10. Fluxes at the interfaces
	67	!!
	68	!! ** Inputs / Ouputs : (global commons)
	69	!! surface temperature : t_su_1d
	70	!! ice/snow temperatures : t_i_1d, t_s_1d
	71	!! ice salinities : sz_i_1d
	72	!! number of layers in the ice/snow : nlay_i, nlay_s
	73	!! total ice/snow thickness : h_i_1d, h_s_1d
	74	!!-------------------------------------------------------------------
[10534]	75	INTEGER, INTENT(in) :: k_cnd ! conduction flux (off, on, emulated)
[8984]	76	!
	77	INTEGER :: ji, jk ! spatial loop index
	78	INTEGER :: jm ! current reference number of equation
	79	INTEGER :: jm_mint, jm_maxt
	80	INTEGER :: iconv ! number of iterations in iterative procedure
	81	INTEGER :: iconv_max = 50 ! max number of iterations in iterative procedure
	82	!
	83	INTEGER, DIMENSION(jpij) :: jm_min ! reference number of top equation
	84	INTEGER, DIMENSION(jpij) :: jm_max ! reference number of bottom equation
[10425]	85
	86	LOGICAL, DIMENSION(jpij) :: l_T_converged ! true when T converges (per grid point)
	87	!
[8984]	88	REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system
	89	REAL(wp) :: zg1 = 2._wp !
	90	REAL(wp) :: zgamma = 18009._wp ! for specific heat
	91	REAL(wp) :: zbeta = 0.117_wp ! for thermal conductivity (could be 0.13)
	92	REAL(wp) :: zraext_s = 10._wp ! extinction coefficient of radiation in the snow
	93	REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity
	94	REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C
	95	REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature
[9423]	96	REAL(wp) :: zhs_min = 0.01_wp ! minimum snow thickness for conductivity calculation
[9935]	97	REAL(wp) :: ztmelts ! ice melting temperature
[8984]	98	REAL(wp) :: zdti_max ! current maximal error on temperature
	99	REAL(wp) :: zcpi ! Ice specific heat
	100	REAL(wp) :: zhfx_err, zdq ! diag errors on heat
	101	REAL(wp) :: zfac ! dummy factor
	102	!
	103	REAL(wp), DIMENSION(jpij) :: isnow ! switch for presence (1) or absence (0) of snow
[13046]	104	REAL(wp), DIMENSION(jpij) :: snow_frac ! fraction of sea ice that is snow covered (for coupled variables only)
[8984]	105	REAL(wp), DIMENSION(jpij) :: ztsub ! surface temperature at previous iteration
	106	REAL(wp), DIMENSION(jpij) :: zh_i, z1_h_i ! ice layer thickness
	107	REAL(wp), DIMENSION(jpij) :: zh_s, z1_h_s ! snow layer thickness
	108	REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface
	109	REAL(wp), DIMENSION(jpij) :: zfnet ! surface flux function
	110	REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function
	111	!
	112	REAL(wp), DIMENSION(jpij ) :: ztsuold ! Old surface temperature in the ice
	113	REAL(wp), DIMENSION(jpij,nlay_i) :: ztiold ! Old temperature in the ice
	114	REAL(wp), DIMENSION(jpij,nlay_s) :: ztsold ! Old temperature in the snow
	115	REAL(wp), DIMENSION(jpij,nlay_i) :: ztib ! Temporary temperature in the ice to check the convergence
	116	REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow to check the convergence
	117	REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity
[10425]	118	REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i_cp ! copy
[8984]	119	REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice
	120	REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice
	121	REAL(wp), DIMENSION(jpij,0:nlay_i) :: zkappa_i ! Kappa factor in the ice
	122	REAL(wp), DIMENSION(jpij,0:nlay_i) :: zeta_i ! Eta factor in the ice
	123	REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow
	124	REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow
	125	REAL(wp), DIMENSION(jpij,0:nlay_s) :: zkappa_s ! Kappa factor in the snow
	126	REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeta_s ! Eta factor in the snow
[13046]	127	REAL(wp), DIMENSION(jpij) :: zkappa_combi ! Combined snow and ice surface conductivity
[8984]	128	REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindterm ! 'Ind'ependent term
	129	REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindtbis ! Temporary 'ind'ependent term
	130	REAL(wp), DIMENSION(jpij,nlay_i+3) :: zdiagbis ! Temporary 'dia'gonal term
	131	REAL(wp), DIMENSION(jpij,nlay_i+3,3) :: ztrid ! Tridiagonal system terms
	132	REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat
	133	REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat
	134	!
	135	! Mono-category
	136	REAL(wp) :: zepsilon ! determines thres. above which computation of G(h) is done
	137	REAL(wp) :: zhe ! dummy factor
	138	REAL(wp) :: zcnd_i ! mean sea ice thermal conductivity
	139	!!------------------------------------------------------------------
	140
	141	! --- diag error on heat diffusion - PART 1 --- !
	142	DO ji = 1, npti
	143	zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + &
	144	& SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
	145	END DO
	146
	147	!------------------
	148	! 1) Initialization
	149	!------------------
	150	DO ji = 1, npti
	151	isnow(ji) = 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) ) ! is there snow or not
[13046]	152	IF( h_s_1d(ji) < zhs_min ) THEN
	153	snow_frac(ji) = h_s_1d(ji) / zhs_min
	154	ELSE
	155	snow_frac(ji) = 1.0_wp
	156	END IF
[8984]	157	! layer thickness
	158	zh_i(ji) = h_i_1d(ji) * r1_nlay_i
	159	zh_s(ji) = h_s_1d(ji) * r1_nlay_s
	160	END DO
	161	!
	162	WHERE( zh_i(1:npti) >= epsi10 ) ; z1_h_i(1:npti) = 1._wp / zh_i(1:npti)
	163	ELSEWHERE ; z1_h_i(1:npti) = 0._wp
	164	END WHERE
	165	!
[9423]	166	WHERE( zh_s(1:npti) > 0._wp ) zh_s(1:npti) = MAX( zhs_min * r1_nlay_s, zh_s(1:npti) )
	167	!
	168	WHERE( zh_s(1:npti) > 0._wp ) ; z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
[8984]	169	ELSEWHERE ; z1_h_s(1:npti) = 0._wp
	170	END WHERE
	171	!
	172	! Store initial temperatures and non solar heat fluxes
[10534]	173	IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
[8984]	174	!
	175	ztsub (1:npti) = t_su_1d(1:npti) ! surface temperature at iteration n-1
	176	ztsuold (1:npti) = t_su_1d(1:npti) ! surface temperature initial value
	177	t_su_1d (1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err ) ! required to leave the choice between melting or not
	178	zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti) ! derivative of incoming nonsolar flux
	179	zqns_ice_b (1:npti) = qns_ice_1d(1:npti) ! store previous qns_ice_1d value
	180	!
	181	ENDIF
	182	!
	183	ztsold (1:npti,:) = t_s_1d(1:npti,:) ! Old snow temperature
	184	ztiold (1:npti,:) = t_i_1d(1:npti,:) ! Old ice temperature
	185
	186	!-------------
	187	! 2) Radiation
	188	!-------------
	189	! --- Transmission/absorption of solar radiation in the ice --- !
[9910]	190	zradtr_s(1:npti,0) = qtr_ice_top_1d(1:npti)
[8984]	191	DO jk = 1, nlay_s
	192	DO ji = 1, npti
	193	! ! radiation transmitted below the layer-th snow layer
[9423]	194	zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * h_s_1d(ji) * r1_nlay_s * REAL(jk) )
[8984]	195	! ! radiation absorbed by the layer-th snow layer
	196	zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
	197	END DO
	198	END DO
	199	!
[9910]	200	zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - isnow(1:npti) )
[8984]	201	DO jk = 1, nlay_i
	202	DO ji = 1, npti
	203	! ! radiation transmitted below the layer-th ice layer
	204	zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )
	205	! ! radiation absorbed by the layer-th ice layer
	206	zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
	207	END DO
	208	END DO
	209	!
[9910]	210	qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti,nlay_i) ! record radiation transmitted below the ice
[8984]	211	!
	212	iconv = 0 ! number of iterations
	213	!
[10425]	214	l_T_converged(:) = .FALSE.
	215	! Convergence calculated until all sub-domain grid points have converged
	216	! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation)
	217	! but values are not taken into account (results independant of MPI partitioning)
	218	!
[10926]	219	! !============================!
[10425]	220	DO WHILE ( ( .NOT. ALL (l_T_converged(1:npti)) ) .AND. iconv < iconv_max ) ! Iterative procedure begins !
[10926]	221	! !============================!
[8984]	222	iconv = iconv + 1
	223	!
	224	ztib(1:npti,:) = t_i_1d(1:npti,:)
	225	ztsb(1:npti,:) = t_s_1d(1:npti,:)
	226	!
	227	!--------------------------------
	228	! 3) Sea ice thermal conductivity
	229	!--------------------------------
	230	IF( ln_cndi_U64 ) THEN !-- Untersteiner (1964) formula: k = k0 + beta.S/T
	231	!
	232	DO ji = 1, npti
[10425]	233	ztcond_i_cp(ji,0) = rcnd_i + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
	234	ztcond_i_cp(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )
[8984]	235	END DO
	236	DO jk = 1, nlay_i-1
	237	DO ji = 1, npti
[10425]	238	ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &
	239	& MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
[8984]	240	END DO
	241	END DO
	242	!
	243	ELSEIF( ln_cndi_P07 ) THEN !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
	244	!
	245	DO ji = 1, npti
[10425]	246	ztcond_i_cp(ji,0) = rcnd_i + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) &
	247	& - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
	248	ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) &
	249	& - 0.011_wp * ( t_bo_1d(ji) - rt0 )
[8984]	250	END DO
	251	DO jk = 1, nlay_i-1
	252	DO ji = 1, npti
[10425]	253	ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &
	254	& MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) &
	255	& - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )
[8984]	256	END DO
	257	END DO
	258	!
	259	ENDIF
[10425]	260
	261	! Variable used after iterations
	262	! Value must be frozen after convergence for MPP independance reason
	263	DO ji = 1, npti
	264	IF ( .NOT. l_T_converged(ji) ) &
	265	ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) )
	266	END DO
[8984]	267	!
	268	!--- G(he) : enhancement of thermal conductivity in mono-category case
	269	! Computation of effective thermal conductivity G(h)
	270	! Used in mono-category case only to simulate an ITD implicitly
	271	! Fichefet and Morales Maqueda, JGR 1997
	272	zghe(1:npti) = 1._wp
	273	!
[10531]	274	IF( ln_virtual_itd ) THEN
[8984]	275	!
	276	zepsilon = 0.1_wp
	277	DO ji = 1, npti
	278	zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp ) ! Mean sea ice thermal conductivity
	279	zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i ) ! Effective thickness he (zhe)
	280	IF( zhe >= zepsilon * 0.5_wp * EXP(1._wp) ) &
	281	& zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) ) ! G(he)
	282	END DO
	283	!
[10531]	284	ENDIF
[8984]	285	!
	286	!-----------------
	287	! 4) kappa factors
	288	!-----------------
	289	!--- Snow
[10425]	290	! Variable used after iterations
	291	! Value must be frozen after convergence for MPP independance reason
[8984]	292	DO jk = 0, nlay_s-1
	293	DO ji = 1, npti
[10425]	294	IF ( .NOT. l_T_converged(ji) ) &
	295	zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
[8984]	296	END DO
	297	END DO
	298	DO ji = 1, npti ! Snow-ice interface
[10425]	299	IF ( .NOT. l_T_converged(ji) ) THEN
	300	zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
	301	IF( zfac > epsi10 ) THEN
	302	zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
	303	ELSE
	304	zkappa_s(ji,nlay_s) = 0._wp
	305	ENDIF
[8984]	306	ENDIF
	307	END DO
	308
	309	!--- Ice
[10425]	310	! Variable used after iterations
	311	! Value must be frozen after convergence for MPP independance reason
[8984]	312	DO jk = 0, nlay_i
	313	DO ji = 1, npti
[10425]	314	IF ( .NOT. l_T_converged(ji) ) &
	315	zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
[8984]	316	END DO
	317	END DO
	318	DO ji = 1, npti ! Snow-ice interface
[10425]	319	IF ( .NOT. l_T_converged(ji) ) &
[13046]	320	! Calculate combined surface snow and ice conductivity to pass through the coupler
	321	zkappa_combi(ji) = zkappa_s(ji,0) * snow_frac(ji) + zkappa_i(ji,0) * ( 1._wp - snow_frac(ji) )
	322
[10425]	323	zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
[8984]	324	END DO
	325	!
	326	!--------------------------------------
	327	! 5) Sea ice specific heat, eta factors
	328	!--------------------------------------
	329	DO jk = 1, nlay_i
	330	DO ji = 1, npti
[9935]	331	zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
	332	zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi )
[8984]	333	END DO
	334	END DO
	335
	336	DO jk = 1, nlay_s
	337	DO ji = 1, npti
[9935]	338	zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji)
[8984]	339	END DO
	340	END DO
	341	!
	342	!----------------------------------------!
	343	! !
[10534]	344	! Conduction flux is off or emulated !
[8984]	345	! !
	346	!----------------------------------------!
	347	!
[10534]	348	IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
[8984]	349	!
	350	! ==> The original BL99 temperature computation is used
	351	! (with qsr_ice, qns_ice and dqns_ice as inputs)
	352	!
	353	!----------------------------
	354	! 6) surface flux computation
	355	!----------------------------
	356	! update of the non solar flux according to the update in T_su
	357	DO ji = 1, npti
[10425]	358	! Variable used after iterations
	359	! Value must be frozen after convergence for MPP independance reason
	360	IF ( .NOT. l_T_converged(ji) ) &
	361	qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
[8984]	362	END DO
	363
	364	DO ji = 1, npti
[9910]	365	zfnet(ji) = qsr_ice_1d(ji) - qtr_ice_top_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar
[8984]	366	END DO
	367	!
	368	!----------------------------
	369	! 7) tridiagonal system terms
	370	!----------------------------
	371	! layer denotes the number of the layer in the snow or in the ice
	372	! jm denotes the reference number of the equation in the tridiagonal
	373	! system, terms of tridiagonal system are indexed as following :
	374	! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
	375
	376	! ice interior terms (top equation has the same form as the others)
	377	ztrid (1:npti,:,:) = 0._wp
	378	zindterm(1:npti,:) = 0._wp
	379	zindtbis(1:npti,:) = 0._wp
	380	zdiagbis(1:npti,:) = 0._wp
	381
	382	DO jm = nlay_s + 2, nlay_s + nlay_i
	383	DO ji = 1, npti
	384	jk = jm - nlay_s - 1
	385	ztrid (ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
	386	ztrid (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
	387	ztrid (ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk)
	388	zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
	389	END DO
	390	END DO
	391
	392	jm = nlay_s + nlay_i + 1
	393	DO ji = 1, npti
	394	! ice bottom term
	395	ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1)
	396	ztrid (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
	397	ztrid (ji,jm,3) = 0._wp
	398	zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * &
	399	& ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) )
	400	END DO
	401
	402	DO ji = 1, npti
	403	! !---------------------!
	404	IF( h_s_1d(ji) > 0._wp ) THEN ! snow-covered cells !
	405	! !---------------------!
	406	! snow interior terms (bottom equation has the same form as the others)
	407	DO jm = 3, nlay_s + 1
	408	jk = jm - 1
	409	ztrid (ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
	410	ztrid (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
	411	ztrid (ji,jm,3) = - zeta_s(ji,jk) * zkappa_s(ji,jk)
	412	zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
	413	END DO
	414
	415	! case of only one layer in the ice (ice equation is altered)
	416	IF( nlay_i == 1 ) THEN
	417	ztrid (ji,nlay_s+2,3) = 0._wp
[9068]	418	zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji)
[8984]	419	ENDIF
	420
	421	IF( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting
	422
	423	jm_min(ji) = 1
	424	jm_max(ji) = nlay_i + nlay_s + 1
	425
	426	! surface equation
	427	ztrid (ji,1,1) = 0._wp
	428	ztrid (ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
	429	ztrid (ji,1,3) = zg1s * zkappa_s(ji,0)
	430	zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
	431
	432	! first layer of snow equation
	433	ztrid (ji,2,1) = - zeta_s(ji,1) * zkappa_s(ji,0) * zg1s
	434	ztrid (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
	435	ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1)
	436	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
	437
	438	ELSE !-- case 2 : surface is melting
	439	!
	440	jm_min(ji) = 2
	441	jm_max(ji) = nlay_i + nlay_s + 1
	442
	443	! first layer of snow equation
	444	ztrid (ji,2,1) = 0._wp
	445	ztrid (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
	446	ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1)
	447	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )
	448	ENDIF
	449	! !---------------------!
	450	ELSE ! cells without snow !
	451	! !---------------------!
	452	!
	453	IF( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting
	454	!
	455	jm_min(ji) = nlay_s + 1
	456	jm_max(ji) = nlay_i + nlay_s + 1
	457
	458	! surface equation
	459	ztrid (ji,jm_min(ji),1) = 0._wp
	460	ztrid (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * zg1
	461	ztrid (ji,jm_min(ji),3) = zkappa_i(ji,0) * zg1
	462	zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
	463
	464	! first layer of ice equation
	465	ztrid (ji,jm_min(ji)+1,1) = - zeta_i(ji,1) * zkappa_i(ji,0) * zg1
	466	ztrid (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
	467	ztrid (ji,jm_min(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1)
	468	zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)
	469
	470	! case of only one layer in the ice (surface & ice equations are altered)
	471	IF( nlay_i == 1 ) THEN
	472	ztrid (ji,jm_min(ji),1) = 0._wp
	473	ztrid (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2._wp
	474	ztrid (ji,jm_min(ji),3) = zkappa_i(ji,0) * 2._wp
	475	ztrid (ji,jm_min(ji)+1,1) = - zeta_i(ji,1) * zkappa_i(ji,0) * 2._wp
	476	ztrid (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
	477	ztrid (ji,jm_min(ji)+1,3) = 0._wp
	478	zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji))
	479	ENDIF
	480
	481	ELSE !-- case 2 : surface is melting
	482
	483	jm_min(ji) = nlay_s + 2
	484	jm_max(ji) = nlay_i + nlay_s + 1
	485
	486	! first layer of ice equation
	487	ztrid (ji,jm_min(ji),1) = 0._wp
	488	ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
	489	ztrid (ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)
	490	zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji))
	491
	492	! case of only one layer in the ice (surface & ice equations are altered)
	493	IF( nlay_i == 1 ) THEN
	494	ztrid (ji,jm_min(ji),1) = 0._wp
	495	ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
	496	ztrid (ji,jm_min(ji),3) = 0._wp
	497	zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &
	498	& + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2._wp
	499	ENDIF
	500
	501	ENDIF
	502	ENDIF
	503	!
	504	zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
	505	zdiagbis(ji,jm_min(ji)) = ztrid (ji,jm_min(ji),2)
	506	!
	507	END DO
	508	!
	509	!------------------------------
	510	! 8) tridiagonal system solving
	511	!------------------------------
	512	! Solve the tridiagonal system with Gauss elimination method.
	513	! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
	514	jm_maxt = 0
	515	jm_mint = nlay_i+5
	516	DO ji = 1, npti
	517	jm_mint = MIN(jm_min(ji),jm_mint)
	518	jm_maxt = MAX(jm_max(ji),jm_maxt)
	519	END DO
	520
	521	DO jk = jm_mint+1, jm_maxt
	522	DO ji = 1, npti
	523	jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
	524	zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1)
	525	zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1)
	526	END DO
	527	END DO
	528
	529	! ice temperatures
	530	DO ji = 1, npti
[10425]	531	! Variable used after iterations
	532	! Value must be frozen after convergence for MPP independance reason
	533	IF ( .NOT. l_T_converged(ji) ) &
	534	t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
[8984]	535	END DO
	536
	537	DO jm = nlay_i + nlay_s, nlay_s + 2, -1
	538	DO ji = 1, npti
	539	jk = jm - nlay_s - 1
[10425]	540	IF ( .NOT. l_T_converged(ji) ) &
	541	t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
[8984]	542	END DO
	543	END DO
	544
	545	DO ji = 1, npti
[10425]	546	! Variables used after iterations
	547	! Value must be frozen after convergence for MPP independance reason
	548	IF ( .NOT. l_T_converged(ji) ) THEN
	549	! snow temperatures
	550	IF( h_s_1d(ji) > 0._wp ) THEN
	551	t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)
	552	ENDIF
	553	! surface temperature
	554	ztsub(ji) = t_su_1d(ji)
	555	IF( t_su_1d(ji) < rt0 ) THEN
	556	t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * &
	557	& ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
	558	ENDIF
[8984]	559	ENDIF
	560	END DO
	561	!
	562	!--------------------------------------------------------------
	563	! 9) Has the scheme converged?, end of the iterative procedure
	564	!--------------------------------------------------------------
	565	! check that nowhere it has started to melt
	566	! zdti_max is a measure of error, it has to be under zdti_bnd
[10425]	567
[8984]	568	DO ji = 1, npti
	569
[10425]	570	zdti_max = 0._wp
[8984]	571
[10425]	572	IF ( .NOT. l_T_converged(ji) ) THEN
	573	t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
	574	zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )
	575
	576	t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp )
	577	zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) )
	578
	579	DO jk = 1, nlay_i
	580	ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0
	581	t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
	582	zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
	583	END DO
	584
	585	IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE.
	586
	587	ENDIF
	588
[8984]	589	END DO
	590
	591	!----------------------------------------!
	592	! !
[10534]	593	! Conduction flux is on !
[8984]	594	! !
	595	!----------------------------------------!
	596	!
[10534]	597	ELSEIF( k_cnd == np_cnd_ON ) THEN
[8984]	598	!
	599	! ==> we use a modified BL99 solver with conduction flux (qcn_ice) as forcing term
	600	!
	601	!----------------------------
	602	! 7) tridiagonal system terms
	603	!----------------------------
	604	! layer denotes the number of the layer in the snow or in the ice
	605	! jm denotes the reference number of the equation in the tridiagonal
	606	! system, terms of tridiagonal system are indexed as following :
	607	! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
	608
	609	! ice interior terms (top equation has the same form as the others)
	610	ztrid (1:npti,:,:) = 0._wp
	611	zindterm(1:npti,:) = 0._wp
	612	zindtbis(1:npti,:) = 0._wp
	613	zdiagbis(1:npti,:) = 0._wp
	614
	615	DO jm = nlay_s + 2, nlay_s + nlay_i
	616	DO ji = 1, npti
	617	jk = jm - nlay_s - 1
	618	ztrid (ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
	619	ztrid (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
	620	ztrid (ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk)
	621	zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
	622	END DO
	623	ENDDO
	624
	625	jm = nlay_s + nlay_i + 1
	626	DO ji = 1, npti
	627	! ice bottom term
	628	ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1)
	629	ztrid (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
	630	ztrid (ji,jm,3) = 0._wp
	631	zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * &
	632	& ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) )
	633	ENDDO
	634
	635	DO ji = 1, npti
	636	! !---------------------!
	637	IF( h_s_1d(ji) > 0._wp ) THEN ! snow-covered cells !
	638	! !---------------------!
	639	! snow interior terms (bottom equation has the same form as the others)
	640	DO jm = 3, nlay_s + 1
	641	jk = jm - 1
	642	ztrid (ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
	643	ztrid (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
	644	ztrid (ji,jm,3) = - zeta_s(ji,jk) * zkappa_s(ji,jk)
	645	zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
	646	END DO
	647
	648	! case of only one layer in the ice (ice equation is altered)
	649	IF ( nlay_i == 1 ) THEN
	650	ztrid (ji,nlay_s+2,3) = 0._wp
[9068]	651	zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji)
[8984]	652	ENDIF
	653
	654	jm_min(ji) = 2
	655	jm_max(ji) = nlay_i + nlay_s + 1
	656
	657	! first layer of snow equation
	658	ztrid (ji,2,1) = 0._wp
	659	ztrid (ji,2,2) = 1._wp + zeta_s(ji,1) * zkappa_s(ji,1)
	660	ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1)
	661	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) )
	662
	663	! !---------------------!
	664	ELSE ! cells without snow !
	665	! !---------------------!
	666	jm_min(ji) = nlay_s + 2
	667	jm_max(ji) = nlay_i + nlay_s + 1
	668
	669	! first layer of ice equation
	670	ztrid (ji,jm_min(ji),1) = 0._wp
	671	ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)
	672	ztrid (ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)
	673	zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) )
	674
	675	! case of only one layer in the ice (surface & ice equations are altered)
	676	IF( nlay_i == 1 ) THEN
	677	ztrid (ji,jm_min(ji),1) = 0._wp
	678	ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)
	679	ztrid (ji,jm_min(ji),3) = 0._wp
	680	zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * &
	681	& ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) + qcn_ice_1d(ji) )
	682	ENDIF
	683
	684	ENDIF
	685	!
	686	zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
	687	zdiagbis(ji,jm_min(ji)) = ztrid (ji,jm_min(ji),2)
	688	!
	689	END DO
	690	!
	691	!------------------------------
	692	! 8) tridiagonal system solving
	693	!------------------------------
	694	! Solve the tridiagonal system with Gauss elimination method.
	695	! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
	696	jm_maxt = 0
	697	jm_mint = nlay_i+5
	698	DO ji = 1, npti
	699	jm_mint = MIN(jm_min(ji),jm_mint)
	700	jm_maxt = MAX(jm_max(ji),jm_maxt)
	701	END DO
	702
	703	DO jk = jm_mint+1, jm_maxt
	704	DO ji = 1, npti
	705	jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
	706	zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1)
	707	zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
	708	END DO
	709	END DO
	710
	711	! ice temperatures
[10425]	712	DO ji = 1, npti
	713	! Variable used after iterations
	714	! Value must be frozen after convergence for MPP independance reason
	715	IF ( .NOT. l_T_converged(ji) ) &
	716	t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
[8984]	717	END DO
	718
	719	DO jm = nlay_i + nlay_s, nlay_s + 2, -1
	720	DO ji = 1, npti
[10425]	721	IF ( .NOT. l_T_converged(ji) ) THEN
	722	jk = jm - nlay_s - 1
	723	t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
	724	ENDIF
[8984]	725	END DO
	726	END DO
	727
	728	! snow temperatures
	729	DO ji = 1, npti
[10425]	730	! Variable used after iterations
	731	! Value must be frozen after convergence for MPP independance reason
	732	IF ( .NOT. l_T_converged(ji) ) THEN
	733	IF( h_s_1d(ji) > 0._wp ) THEN
	734	t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)
	735	ENDIF
[8984]	736	ENDIF
	737	END DO
	738	!
	739	!--------------------------------------------------------------
	740	! 9) Has the scheme converged?, end of the iterative procedure
	741	!--------------------------------------------------------------
	742	! check that nowhere it has started to melt
	743	! zdti_max is a measure of error, it has to be under zdti_bnd
	744
[10425]	745	DO ji = 1, npti
	746
	747	zdti_max = 0._wp
	748
	749	IF ( .NOT. l_T_converged(ji) ) THEN
	750	! t_s
	751	t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp )
	752	zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) )
	753	! t_i
[10924]	754	DO jk = 1, nlay_i
[10425]	755	ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0
	756	t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
	757	zdti_max = MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
	758	END DO
	759
	760	IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE.
	761
	762	ENDIF
	763
[8984]	764	END DO
	765
[10534]	766	ENDIF ! k_cnd
[8984]	767
	768	END DO ! End of the do while iterative procedure
	769
	770	IF( ln_icectl .AND. lwp ) THEN
	771	WRITE(numout,*) ' zdti_max : ', zdti_max
	772	WRITE(numout,*) ' iconv : ', iconv
	773	ENDIF
	774
	775	!
	776	!-----------------------------
	777	! 10) Fluxes at the interfaces
	778	!-----------------------------
	779	!
[9916]	780	! --- calculate conduction fluxes (positive downward)
	781
[8984]	782	DO ji = 1, npti
	783	! ! surface ice conduction flux
[9916]	784	qcn_ice_top_1d(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) &
	785	& - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * ( t_i_1d(ji,1) - t_su_1d(ji) )
[8984]	786	! ! bottom ice conduction flux
[9916]	787	qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1 * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) )
[8984]	788	END DO
	789
	790	!
	791	! --- Diagnose the heat loss due to changing non-solar / conduction flux --- !
	792	!
[10534]	793	IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
[8984]	794	!
	795	DO ji = 1, npti
	796	hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji)
	797	END DO
	798	!
[10534]	799	ELSEIF( k_cnd == np_cnd_ON ) THEN
[8984]	800	!
	801	DO ji = 1, npti
[9916]	802	hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qcn_ice_top_1d(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji)
[8984]	803	END DO
	804	!
	805	ENDIF
	806
	807	!
	808	! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- !
	809	!
[10534]	810	IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_ON ) THEN
[8984]	811
	812	CALL ice_var_enthalpy
	813
	814	! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
	815	DO ji = 1, npti
	816	zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + &
	817	& SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
	818
[10534]	819	IF( k_cnd == np_cnd_OFF ) THEN
[8984]	820
	821	IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC
[9916]	822	zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji) &
	823	& + zdq * r1_rdtice ) * a_i_1d(ji)
[8984]	824	ELSE ! case T_su = 0degC
[9916]	825	zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji) &
	826	& + zdq * r1_rdtice ) * a_i_1d(ji)
[8984]	827	ENDIF
	828
[10534]	829	ELSEIF( k_cnd == np_cnd_ON ) THEN
[8984]	830
[9916]	831	zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji) &
	832	& + zdq * r1_rdtice ) * a_i_1d(ji)
[8984]	833
	834	ENDIF
	835	!
	836	! total heat sink to be sent to the ocean
	837	hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
	838	!
	839	! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2
	840	hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
	841	!
	842	END DO
	843	!
	844	ENDIF
	845	!
[10534]	846	!--------------------------------------------------------------------
	847	! 11) reset inner snow and ice temperatures, update conduction fluxes
	848	!--------------------------------------------------------------------
[8984]	849	! effective conductivity and 1st layer temperature (needed by Met Office)
	850	DO ji = 1, npti
[13043]	851	IF( h_i_1d(ji) > 0.1_wp ) THEN
[13046]	852	cnd_ice_1d(ji) = 2._wp * zkappa_combi(ji)
[8984]	853	ELSE
[13043]	854	cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp
[8984]	855	ENDIF
[13046]	856	t1_ice_1d(ji) = snow_frac(ji) * t_s_1d(ji,1) + ( 1._wp - snow_frac(ji) ) * t_i_1d(ji,1)
[8984]	857	END DO
	858	!
[10534]	859	IF( k_cnd == np_cnd_EMU ) THEN
[8984]	860	! Restore temperatures to their initial values
[9916]	861	t_s_1d (1:npti,:) = ztsold (1:npti,:)
	862	t_i_1d (1:npti,:) = ztiold (1:npti,:)
	863	qcn_ice_1d(1:npti) = qcn_ice_top_1d(1:npti)
[10926]	864
	865	!!clem
	866	! remettre t_su_1d, qns_ice_1d et dqns_ice_1d comme avant puisqu'on devrait faire comme si on avant conduction = input
	867	!clem
[8984]	868	ENDIF
	869	!
[9916]	870	! --- SIMIP diagnostics
	871	!
	872	DO ji = 1, npti
	873	!--- Snow-ice interfacial temperature (diagnostic SIMIP)
	874	zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
	875	IF( h_s_1d(ji) >= zhs_min ) THEN
	876	t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + &
	877	& ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac )
	878	ELSE
	879	t_si_1d(ji) = t_su_1d(ji)
	880	ENDIF
	881	END DO
	882	!
[8984]	883	END SUBROUTINE ice_thd_zdf_BL99
	884
	885	#else
	886	!!----------------------------------------------------------------------
[9570]	887	!! Default option Dummy Module No SI3 sea-ice model
[8984]	888	!!----------------------------------------------------------------------
	889	#endif
	890
	891	!!======================================================================
	892	END MODULE icethd_zdf_BL99

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.1_snow_frac_cpl/src/ICE/icethd_zdf_bl99.F90 @ 13046

Download in other formats: