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

Context Navigation

source: branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 5225

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

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

Download in other formats: