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

Context Navigation

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 4332

Last change on this file since 4332 was 4332, checked in by clem, 11 years ago
update LIM3 to fix remaining bugs. Now working in global and regional config.
Property svn:keywords set to `Id`
File size: 37.0 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
[4045]	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)
[921]	27
[825]	28	IMPLICIT NONE
	29	PRIVATE
	30
[2528]	31	PUBLIC lim_thd_dif ! called by lim_thd
[825]	32
[4332]	33	REAL(wp) :: epsi10 = 1.e-10_wp !
[825]	34	!!----------------------------------------------------------------------
[4045]	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
	41	SUBROUTINE lim_thd_dif( kideb , kiut , jl )
[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)
	76	!! surface temperature : t_su_b
	77	!! ice/snow temperatures : t_i_b, t_s_b
	78	!! ice salinities : s_i_b
	79	!! number of layers in the ice/snow: nlay_i, nlay_s
	80	!! profile of the ice/snow layers : z_i, z_s
	81	!! total ice/snow thickness : ht_i_b, ht_s_b
	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	!!------------------------------------------------------------------
[3294]	93	INTEGER , INTENT (in) :: kideb ! Start point on which the the computation is applied
	94	INTEGER , INTENT (in) :: kiut ! End point on which the the computation is applied
	95	INTEGER , INTENT (in) :: jl ! Category number
[825]	96
[921]	97	!! * Local variables
[3294]	98	INTEGER :: ji ! spatial loop index
	99	INTEGER :: ii, ij ! temporary dummy loop index
	100	INTEGER :: numeq ! current reference number of equation
	101	INTEGER :: layer ! vertical dummy loop index
	102	INTEGER :: nconv ! number of iterations in iterative procedure
	103	INTEGER :: minnumeqmin, maxnumeqmax
[3610]	104	INTEGER, DIMENSION(kiut) :: numeqmin ! reference number of top equation
	105	INTEGER, DIMENSION(kiut) :: numeqmax ! reference number of bottom equation
	106	INTEGER, DIMENSION(kiut) :: isnow ! switch for presence (1) or absence (0) of snow
[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)
	111	REAL(wp) :: zraext_s = 1.e+8_wp ! extinction coefficient of radiation in the snow
	112	REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity
[2715]	113	REAL(wp) :: ztmelt_i ! ice melting temperature
	114	REAL(wp) :: zerritmax ! current maximal error on temperature
[3610]	115	REAL(wp), DIMENSION(kiut) :: ztfs ! ice melting point
	116	REAL(wp), DIMENSION(kiut) :: ztsuold ! old surface temperature (before the iterative procedure )
	117	REAL(wp), DIMENSION(kiut) :: ztsuoldit ! surface temperature at previous iteration
	118	REAL(wp), DIMENSION(kiut) :: zh_i ! ice layer thickness
	119	REAL(wp), DIMENSION(kiut) :: zh_s ! snow layer thickness
	120	REAL(wp), DIMENSION(kiut) :: zfsw ! solar radiation absorbed at the surface
	121	REAL(wp), DIMENSION(kiut) :: zf ! surface flux function
	122	REAL(wp), DIMENSION(kiut) :: dzf ! derivative of the surface flux function
	123	REAL(wp), DIMENSION(kiut) :: zerrit ! current error on temperature
	124	REAL(wp), DIMENSION(kiut) :: zdifcase ! case of the equation resolution (1->4)
	125	REAL(wp), DIMENSION(kiut) :: zftrice ! solar radiation transmitted through the ice
	126	REAL(wp), DIMENSION(kiut) :: zihic, zhsu
	127	REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztcond_i ! Ice thermal conductivity
	128	REAL(wp), DIMENSION(kiut,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice
	129	REAL(wp), DIMENSION(kiut,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice
	130	REAL(wp), DIMENSION(kiut,0:nlay_i) :: zkappa_i ! Kappa factor in the ice
	131	REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztiold ! Old temperature in the ice
	132	REAL(wp), DIMENSION(kiut,0:nlay_i) :: zeta_i ! Eta factor in the ice
	133	REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztitemp ! Temporary temperature in the ice to check the convergence
	134	REAL(wp), DIMENSION(kiut,0:nlay_i) :: zspeche_i ! Ice specific heat
	135	REAL(wp), DIMENSION(kiut,0:nlay_i) :: z_i ! Vertical cotes of the layers in the ice
	136	REAL(wp), DIMENSION(kiut,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow
	137	REAL(wp), DIMENSION(kiut,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow
	138	REAL(wp), DIMENSION(kiut,0:nlay_s) :: zkappa_s ! Kappa factor in the snow
	139	REAL(wp), DIMENSION(kiut,0:nlay_s) :: zeta_s ! Eta factor in the snow
	140	REAL(wp), DIMENSION(kiut,0:nlay_s) :: ztstemp ! Temporary temperature in the snow to check the convergence
	141	REAL(wp), DIMENSION(kiut,0:nlay_s) :: ztsold ! Temporary temperature in the snow
	142	REAL(wp), DIMENSION(kiut,0:nlay_s) :: z_s ! Vertical cotes of the layers in the snow
	143	REAL(wp), DIMENSION(kiut,jkmax+2) :: zindterm ! Independent term
	144	REAL(wp), DIMENSION(kiut,jkmax+2) :: zindtbis ! temporary independent term
	145	REAL(wp), DIMENSION(kiut,jkmax+2) :: zdiagbis
	146	REAL(wp), DIMENSION(kiut,jkmax+2,3) :: ztrid ! tridiagonal system terms
[3625]	147	!!------------------------------------------------------------------
[3610]	148	!
[921]	149	!------------------------------------------------------------------------------!
	150	! 1) Initialization !
	151	!------------------------------------------------------------------------------!
	152	!
	153	DO ji = kideb , kiut
	154	! is there snow or not
[4045]	155	isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) ) )
[921]	156	! surface temperature of fusion
[2715]	157	!!gm ??? ztfs(ji) = rtt !!!????
[4045]	158	ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt
[921]	159	! layer thickness
[4045]	160	zh_i(ji) = ht_i_b(ji) / REAL( nlay_i )
	161	zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
[921]	162	END DO
[825]	163
[921]	164	!--------------------
	165	! Ice / snow layers
	166	!--------------------
[825]	167
[2715]	168	z_s(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st snow layer
	169	z_i(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st ice layer
[825]	170
[2715]	171	DO layer = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer
[921]	172	DO ji = kideb , kiut
[4045]	173	z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s )
[921]	174	END DO
	175	END DO
[825]	176
[2715]	177	DO layer = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer
[921]	178	DO ji = kideb , kiut
[4045]	179	z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i )
[921]	180	END DO
	181	END DO
	182	!
	183	!------------------------------------------------------------------------------\|
	184	! 2) Radiations \|
	185	!------------------------------------------------------------------------------\|
	186	!
	187	!-------------------
	188	! Computation of i0
	189	!-------------------
	190	! i0 describes the fraction of solar radiation which does not contribute
	191	! to the surface energy budget but rather penetrates inside the ice.
	192	! We assume that no radiation is transmitted through the snow
	193	! If there is no no snow
	194	! zfsw = (1-i0).qsr_ice is absorbed at the surface
	195	! zftrice = io.qsr_ice is below the surface
	196	! fstbif = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
[825]	197
[921]	198	DO ji = kideb , kiut
	199	! switches
[4045]	200	isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) ) )
[921]	201	! hs > 0, isnow = 1
[2715]	202	zhsu (ji) = hnzst ! threshold for the computation of i0
	203	zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) )
[825]	204
[4045]	205	i0(ji) = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
[921]	206	!fr1_i0_1d = i0 for a thin ice surface
	207	!fr1_i0_2d = i0 for a thick ice surface
	208	! a function of the cloud cover
	209	!
	210	!i0(ji) = (1.0-FLOAT(isnow(ji)))3.0/(100ht_s_b(ji)+10.0)
	211	!formula used in Cice
	212	END DO
[825]	213
[921]	214	!-------------------------------------------------------
	215	! Solar radiation absorbed / transmitted at the surface
	216	! Derivative of the non solar flux
	217	!-------------------------------------------------------
	218	DO ji = kideb , kiut
[2715]	219	zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface
	220	zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer
	221	dzf (ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux
[921]	222	END DO
[825]	223
[921]	224	!---------------------------------------------------------
	225	! Transmission - absorption of solar radiation in the ice
	226	!---------------------------------------------------------
[825]	227
[2715]	228	DO ji = kideb, kiut ! snow initialization
	229	zradtr_s(ji,0) = zftrice(ji) ! radiation penetrating through snow
[921]	230	END DO
[825]	231
[2715]	232	DO layer = 1, nlay_s ! Radiation through snow
	233	DO ji = kideb, kiut
	234	! ! radiation transmitted below the layer-th snow layer
	235	zradtr_s(ji,layer) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,layer) ) ) )
	236	! ! radiation absorbed by the layer-th snow layer
[921]	237	zradab_s(ji,layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer)
	238	END DO
	239	END DO
[825]	240
[2715]	241	DO ji = kideb, kiut ! ice initialization
[4045]	242	zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * REAL( isnow(ji) ) + zftrice(ji) * REAL( 1 - isnow(ji) )
[921]	243	END DO
[825]	244
[2715]	245	DO layer = 1, nlay_i ! Radiation through ice
	246	DO ji = kideb, kiut
	247	! ! radiation transmitted below the layer-th ice layer
	248	zradtr_i(ji,layer) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,layer) ) ) )
	249	! ! radiation absorbed by the layer-th ice layer
[921]	250	zradab_i(ji,layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer)
	251	END DO
	252	END DO
[825]	253
[2715]	254	DO ji = kideb, kiut ! Radiation transmitted below the ice
[4045]	255	fstbif_1d(ji) = fstbif_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) ! clem modif
[921]	256	END DO
[834]	257
[921]	258	! +++++
	259	! just to check energy conservation
[2715]	260	DO ji = kideb, kiut
[3625]	261	ii = MOD( npb(ji) - 1 , jpi ) + 1
	262	ij = ( npb(ji) - 1 ) / jpi + 1
[4045]	263	fstroc(ii,ij,jl) = iatte_1d(ji) * zradtr_i(ji,nlay_i) ! clem modif
[921]	264	END DO
	265	! +++++
[825]	266
[921]	267	DO layer = 1, nlay_i
[2715]	268	DO ji = kideb, kiut
[921]	269	radab(ji,layer) = zradab_i(ji,layer)
	270	END DO
	271	END DO
[825]	272
[921]	273	!
	274	!------------------------------------------------------------------------------\|
	275	! 3) Iterative procedure begins \|
	276	!------------------------------------------------------------------------------\|
	277	!
[2715]	278	DO ji = kideb, kiut ! Old surface temperature
	279	ztsuold (ji) = t_su_b(ji) ! temperature at the beg of iter pr.
	280	ztsuoldit(ji) = t_su_b(ji) ! temperature at the previous iter
	281	t_su_b (ji) = MIN( t_su_b(ji), ztfs(ji)-0.00001 ) ! necessary
	282	zerrit (ji) = 1000._wp ! initial value of error
[921]	283	END DO
[825]	284
[2715]	285	DO layer = 1, nlay_s ! Old snow temperature
[921]	286	DO ji = kideb , kiut
[2715]	287	ztsold(ji,layer) = t_s_b(ji,layer)
[921]	288	END DO
	289	END DO
[825]	290
[2715]	291	DO layer = 1, nlay_i ! Old ice temperature
[921]	292	DO ji = kideb , kiut
[2715]	293	ztiold(ji,layer) = t_i_b(ji,layer)
[921]	294	END DO
	295	END DO
[825]	296
[2715]	297	nconv = 0 ! number of iterations
	298	zerritmax = 1000._wp ! maximal value of error on all points
[825]	299
[2715]	300	DO WHILE ( zerritmax > maxer_i_thd .AND. nconv < nconv_i_thd )
[921]	301	!
[2715]	302	nconv = nconv + 1
	303	!
[921]	304	!------------------------------------------------------------------------------\|
	305	! 4) Sea ice thermal conductivity \|
	306	!------------------------------------------------------------------------------\|
	307	!
[2715]	308	IF( thcon_i_swi == 0 ) THEN ! Untersteiner (1964) formula
[921]	309	DO ji = kideb , kiut
[4332]	310	ztcond_i(ji,0) = rcdic + zbeta*s_i_b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt)
[921]	311	ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin)
	312	END DO
	313	DO layer = 1, nlay_i-1
	314	DO ji = kideb , kiut
[2777]	315	ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) / &
[4332]	316	MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)
[921]	317	ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin)
	318	END DO
	319	END DO
	320	ENDIF
[825]	321
[2715]	322	IF( thcon_i_swi == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
[921]	323	DO ji = kideb , kiut
[4332]	324	ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt ) &
[2715]	325	& - 0.011_wp * ( t_i_b(ji,1) - rtt )
	326	ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
[921]	327	END DO
[2715]	328	DO layer = 1, nlay_i-1
	329	DO ji = kideb , kiut
	330	ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) &
[4332]	331	& / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) &
[2715]	332	& - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )
	333	ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin )
	334	END DO
	335	END DO
	336	DO ji = kideb , kiut
[4332]	337	ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt) &
[2715]	338	& - 0.011_wp * ( t_bo_b(ji) - rtt )
	339	ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
	340	END DO
[921]	341	ENDIF
	342	!
	343	!------------------------------------------------------------------------------\|
	344	! 5) kappa factors \|
	345	!------------------------------------------------------------------------------\|
	346	!
	347	DO ji = kideb, kiut
[825]	348
[921]	349	!-- Snow kappa factors
[4332]	350	zkappa_s(ji,0) = rcdsn / MAX(epsi10,zh_s(ji))
	351	zkappa_s(ji,nlay_s) = rcdsn / MAX(epsi10,zh_s(ji))
[921]	352	END DO
[825]	353
[921]	354	DO layer = 1, nlay_s-1
	355	DO ji = kideb , kiut
	356	zkappa_s(ji,layer) = 2.0 * rcdsn / &
[4332]	357	MAX(epsi10,2.0*zh_s(ji))
[921]	358	END DO
	359	END DO
[825]	360
[921]	361	DO layer = 1, nlay_i-1
	362	DO ji = kideb , kiut
	363	!-- Ice kappa factors
	364	zkappa_i(ji,layer) = 2.0*ztcond_i(ji,layer)/ &
[4332]	365	MAX(epsi10,2.0*zh_i(ji))
[921]	366	END DO
	367	END DO
[825]	368
[921]	369	DO ji = kideb , kiut
[4332]	370	zkappa_i(ji,0) = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji))
	371	zkappa_i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji))
[921]	372	!-- Interface
[4332]	373	zkappa_s(ji,nlay_s) = 2.0rcdsnztcond_i(ji,0)/MAX(epsi10, &
[921]	374	(ztcond_i(ji,0)zh_s(ji) + rcdsnzh_i(ji)))
[4045]	375	zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) &
	376	+ zkappa_i(ji,0)*REAL( 1 - isnow(ji) )
[921]	377	END DO
	378	!
	379	!------------------------------------------------------------------------------\|
	380	! 6) Sea ice specific heat, eta factors \|
	381	!------------------------------------------------------------------------------\|
	382	!
	383	DO layer = 1, nlay_i
	384	DO ji = kideb , kiut
	385	ztitemp(ji,layer) = t_i_b(ji,layer)
	386	zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ &
[4332]	387	MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10)
[921]	388	zeta_i(ji,layer) = rdt_ice / MAX(rhoiczspeche_i(ji,layer)zh_i(ji), &
[4332]	389	epsi10)
[921]	390	END DO
	391	END DO
[825]	392
[921]	393	DO layer = 1, nlay_s
	394	DO ji = kideb , kiut
	395	ztstemp(ji,layer) = t_s_b(ji,layer)
[4332]	396	zeta_s(ji,layer) = rdt_ice / MAX(rhosncpiczh_s(ji),epsi10)
[921]	397	END DO
	398	END DO
	399	!
	400	!------------------------------------------------------------------------------\|
	401	! 7) surface flux computation \|
	402	!------------------------------------------------------------------------------\|
	403	!
	404	DO ji = kideb , kiut
[825]	405
[921]	406	! update of the non solar flux according to the update in T_su
	407	qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * &
	408	( t_su_b(ji) - ztsuoldit(ji) )
[825]	409
[921]	410	! update incoming flux
	411	zf(ji) = zfsw(ji) & ! net absorbed solar radiation
	412	+ qnsr_ice_1d(ji) ! non solar total flux
	413	! (LWup, LWdw, SH, LH)
[825]	414
[921]	415	END DO
[825]	416
[921]	417	!
	418	!------------------------------------------------------------------------------\|
	419	! 8) tridiagonal system terms \|
	420	!------------------------------------------------------------------------------\|
	421	!
	422	!!layer denotes the number of the layer in the snow or in the ice
	423	!!numeq denotes the reference number of the equation in the tridiagonal
	424	!!system, terms of tridiagonal system are indexed as following :
	425	!!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
[825]	426
[921]	427	!!ice interior terms (top equation has the same form as the others)
	428
	429	DO numeq=1,jkmax+2
	430	DO ji = kideb , kiut
	431	ztrid(ji,numeq,1) = 0.
	432	ztrid(ji,numeq,2) = 0.
	433	ztrid(ji,numeq,3) = 0.
	434	zindterm(ji,numeq)= 0.
	435	zindtbis(ji,numeq)= 0.
	436	zdiagbis(ji,numeq)= 0.
	437	ENDDO
	438	ENDDO
	439
	440	DO numeq = nlay_s + 2, nlay_s + nlay_i
	441	DO ji = kideb , kiut
	442	layer = numeq - nlay_s - 1
	443	ztrid(ji,numeq,1) = - zeta_i(ji,layer)*zkappa_i(ji,layer-1)
	444	ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,layer)*(zkappa_i(ji,layer-1) + &
	445	zkappa_i(ji,layer))
	446	ztrid(ji,numeq,3) = - zeta_i(ji,layer)*zkappa_i(ji,layer)
	447	zindterm(ji,numeq) = ztiold(ji,layer) + zeta_i(ji,layer)* &
	448	zradab_i(ji,layer)
	449	END DO
	450	ENDDO
	451
	452	numeq = nlay_s + nlay_i + 1
	453	DO ji = kideb , kiut
[825]	454	!!ice bottom term
	455	ztrid(ji,numeq,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)
	456	ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,nlay_i)( zkappa_i(ji,nlay_i)zg1 &
[921]	457	+ zkappa_i(ji,nlay_i-1) )
[825]	458	ztrid(ji,numeq,3) = 0.0
	459	zindterm(ji,numeq) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i)* &
[921]	460	( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 &
	461	* t_bo_b(ji) )
	462	ENDDO
[825]	463
	464
[921]	465	DO ji = kideb , kiut
[825]	466	IF ( ht_s_b(ji).gt.0.0 ) THEN
[921]	467	!
	468	!------------------------------------------------------------------------------\|
	469	! snow-covered cells \|
	470	!------------------------------------------------------------------------------\|
	471	!
	472	!!snow interior terms (bottom equation has the same form as the others)
	473	DO numeq = 3, nlay_s + 1
	474	layer = numeq - 1
	475	ztrid(ji,numeq,1) = - zeta_s(ji,layer)*zkappa_s(ji,layer-1)
	476	ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,layer)*( zkappa_s(ji,layer-1) + &
	477	zkappa_s(ji,layer) )
	478	ztrid(ji,numeq,3) = - zeta_s(ji,layer)*zkappa_s(ji,layer)
	479	zindterm(ji,numeq) = ztsold(ji,layer) + zeta_s(ji,layer)* &
	480	zradab_s(ji,layer)
	481	END DO
[825]	482
[921]	483	!!case of only one layer in the ice (ice equation is altered)
	484	IF ( nlay_i.eq.1 ) THEN
	485	ztrid(ji,nlay_s+2,3) = 0.0
	486	zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* &
	487	t_bo_b(ji)
	488	ENDIF
[834]	489
[921]	490	IF ( t_su_b(ji) .LT. rtt ) THEN
[825]	491
[921]	492	!------------------------------------------------------------------------------\|
	493	! case 1 : no surface melting - snow present \|
	494	!------------------------------------------------------------------------------\|
	495	zdifcase(ji) = 1.0
	496	numeqmin(ji) = 1
	497	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	498
[921]	499	!!surface equation
	500	ztrid(ji,1,1) = 0.0
	501	ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0)
	502	ztrid(ji,1,3) = zg1s*zkappa_s(ji,0)
	503	zindterm(ji,1) = dzf(ji)*t_su_b(ji) - zf(ji)
[825]	504
[921]	505	!!first layer of snow equation
	506	ztrid(ji,2,1) = - zkappa_s(ji,0)zg1szeta_s(ji,1)
	507	ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)(zkappa_s(ji,1) + zkappa_s(ji,0)zg1s)
	508	ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1)
	509	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)
[825]	510
[921]	511	ELSE
	512	!
	513	!------------------------------------------------------------------------------\|
	514	! case 2 : surface is melting - snow present \|
	515	!------------------------------------------------------------------------------\|
	516	!
	517	zdifcase(ji) = 2.0
	518	numeqmin(ji) = 2
	519	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	520
[921]	521	!!first layer of snow equation
	522	ztrid(ji,2,1) = 0.0
	523	ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + &
	524	zkappa_s(ji,0) * zg1s )
	525	ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1)
	526	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * &
	527	( zradab_s(ji,1) + &
	528	zkappa_s(ji,0) * zg1s * t_su_b(ji) )
	529	ENDIF
	530	ELSE
	531	!
	532	!------------------------------------------------------------------------------\|
	533	! cells without snow \|
	534	!------------------------------------------------------------------------------\|
	535	!
	536	IF (t_su_b(ji) .LT. rtt) THEN
	537	!
	538	!------------------------------------------------------------------------------\|
	539	! case 3 : no surface melting - no snow \|
	540	!------------------------------------------------------------------------------\|
	541	!
	542	zdifcase(ji) = 3.0
	543	numeqmin(ji) = nlay_s + 1
	544	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	545
[921]	546	!!surface equation
	547	ztrid(ji,numeqmin(ji),1) = 0.0
	548	ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1
	549	ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1
	550	zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_b(ji) - zf(ji)
[825]	551
[921]	552	!!first layer of ice equation
	553	ztrid(ji,numeqmin(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
	554	ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) &
	555	+ zkappa_i(ji,0) * zg1 )
	556	ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1)
	557	zindterm(ji,numeqmin(ji)+1)= ztiold(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)
[825]	558
[921]	559	!!case of only one layer in the ice (surface & ice equations are altered)
[825]	560
[921]	561	IF (nlay_i.eq.1) THEN
	562	ztrid(ji,numeqmin(ji),1) = 0.0
	563	ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*2.0
	564	ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*2.0
	565	ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0)2.0zeta_i(ji,1)
	566	ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,0)2.0 + &
	567	zkappa_i(ji,1))
	568	ztrid(ji,numeqmin(ji)+1,3) = 0.0
[825]	569
[921]	570	zindterm(ji,numeqmin(ji)+1) = ztiold(ji,1) + zeta_i(ji,1)* &
	571	( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji) )
	572	ENDIF
[825]	573
[921]	574	ELSE
[825]	575
[921]	576	!
	577	!------------------------------------------------------------------------------\|
	578	! case 4 : surface is melting - no snow \|
	579	!------------------------------------------------------------------------------\|
	580	!
	581	zdifcase(ji) = 4.0
	582	numeqmin(ji) = nlay_s + 2
	583	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	584
[921]	585	!!first layer of ice equation
	586	ztrid(ji,numeqmin(ji),1) = 0.0
	587	ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,1) + zkappa_i(ji,0) &
	588	zg1)
	589	ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)
	590	zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &
	591	zkappa_i(ji,0) * zg1 * t_su_b(ji) )
[825]	592
[921]	593	!!case of only one layer in the ice (surface & ice equations are altered)
	594	IF (nlay_i.eq.1) THEN
	595	ztrid(ji,numeqmin(ji),1) = 0.0
	596	ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,0)2.0 + &
	597	zkappa_i(ji,1))
	598	ztrid(ji,numeqmin(ji),3) = 0.0
	599	zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)* &
	600	(zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji)) &
	601	+ t_su_b(ji)zeta_i(ji,1)zkappa_i(ji,0)*2.0
	602	ENDIF
[825]	603
[921]	604	ENDIF
	605	ENDIF
[825]	606
[921]	607	END DO
[825]	608
[921]	609	!
	610	!------------------------------------------------------------------------------\|
	611	! 9) tridiagonal system solving \|
	612	!------------------------------------------------------------------------------\|
	613	!
[825]	614
[921]	615	! Solve the tridiagonal system with Gauss elimination method.
	616	! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
	617	! McGraw-Hill 1984.
[825]	618
[921]	619	maxnumeqmax = 0
	620	minnumeqmin = jkmax+4
[825]	621
[921]	622	DO ji = kideb , kiut
	623	zindtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji))
	624	zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2)
	625	minnumeqmin = MIN(numeqmin(ji),minnumeqmin)
	626	maxnumeqmax = MAX(numeqmax(ji),maxnumeqmax)
	627	END DO
	628
	629	DO layer = minnumeqmin+1, maxnumeqmax
	630	DO ji = kideb , kiut
	631	numeq = min(max(numeqmin(ji)+1,layer),numeqmax(ji))
	632	zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* &
	633	ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1)
	634	zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1)* &
	635	zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1)
	636	END DO
	637	END DO
	638
	639	DO ji = kideb , kiut
	640	! ice temperatures
	641	t_i_b(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))
	642	END DO
	643
	644	DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1
	645	DO ji = kideb , kiut
	646	layer = numeq - nlay_s - 1
	647	t_i_b(ji,layer) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &
	648	t_i_b(ji,layer+1))/zdiagbis(ji,numeq)
	649	END DO
	650	END DO
	651
	652	DO ji = kideb , kiut
[825]	653	! snow temperatures
[4220]	654	IF (ht_s_b(ji).GT.0._wp) &
[921]	655	t_s_b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
	656	* t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) &
[4045]	657	* MAX(0.0,SIGN(1.0,ht_s_b(ji)))
[825]	658
	659	! surface temperature
[4045]	660	isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) ) ) )
[2715]	661	ztsuoldit(ji) = t_su_b(ji)
[4045]	662	IF( t_su_b(ji) < ztfs(ji) ) &
	663	t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1) &
	664	& + REAL( 1 - isnow(ji) )*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))
[921]	665	END DO
	666	!
	667	!--------------------------------------------------------------------------
	668	! 10) Has the scheme converged ?, end of the iterative procedure \|
	669	!--------------------------------------------------------------------------
	670	!
	671	! check that nowhere it has started to melt
	672	! zerrit(ji) is a measure of error, it has to be under maxer_i_thd
	673	DO ji = kideb , kiut
[2715]	674	t_su_b(ji) = MAX( MIN( t_su_b(ji) , ztfs(ji) ) , 190._wp )
	675	zerrit(ji) = ABS( t_su_b(ji) - ztsuoldit(ji) )
[921]	676	END DO
[825]	677
[921]	678	DO layer = 1, nlay_s
	679	DO ji = kideb , kiut
[2715]	680	ii = MOD( npb(ji) - 1, jpi ) + 1
	681	ij = ( npb(ji) - 1 ) / jpi + 1
	682	t_s_b(ji,layer) = MAX( MIN( t_s_b(ji,layer), rtt ), 190._wp )
	683	zerrit(ji) = MAX(zerrit(ji),ABS(t_s_b(ji,layer) - ztstemp(ji,layer)))
[921]	684	END DO
	685	END DO
[825]	686
[921]	687	DO layer = 1, nlay_i
	688	DO ji = kideb , kiut
[2715]	689	ztmelt_i = -tmut * s_i_b(ji,layer) + rtt
[4220]	690	t_i_b(ji,layer) = MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp)
[2715]	691	zerrit(ji) = MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer)))
[921]	692	END DO
	693	END DO
[825]	694
[921]	695	! Compute spatial maximum over all errors
[2715]	696	! note that this could be optimized substantially by iterating only the non-converging points
	697	zerritmax = 0._wp
	698	DO ji = kideb, kiut
	699	zerritmax = MAX( zerritmax, zerrit(ji) )
[921]	700	END DO
[2715]	701	IF( lk_mpp ) CALL mpp_max( zerritmax, kcom=ncomm_ice )
[825]	702
	703	END DO ! End of the do while iterative procedure
	704
[4332]	705	IF( ln_nicep .AND. lwp ) THEN
[1055]	706	WRITE(numout,*) ' zerritmax : ', zerritmax
	707	WRITE(numout,*) ' nconv : ', nconv
	708	ENDIF
[825]	709
[921]	710	!
[2715]	711	!-------------------------------------------------------------------------!
	712	! 11) Fluxes at the interfaces !
	713	!-------------------------------------------------------------------------!
[921]	714	DO ji = kideb, kiut
[3808]	715	#if ! defined key_coupled
	716	! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)
	717	qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) )
	718	#endif
[2715]	719	! ! surface ice conduction flux
[4045]	720	isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) ) )
	721	fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji)) &
	722	& - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_b(ji,1) - t_su_b(ji))
[2715]	723	! ! bottom ice conduction flux
	724	fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
[921]	725	END DO
[825]	726
[921]	727	!-------------------------!
	728	! Heat conservation !
	729	!-------------------------!
[4332]	730	IF( con_i .AND. jiindex_1d > 0 ) THEN
[921]	731	DO ji = kideb, kiut
	732	! Upper snow value
[4045]	733	fc_s(ji,0) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) )
[921]	734	! Bott. snow value
[4045]	735	fc_s(ji,1) = - REAL( isnow(ji) ) * zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) )
[921]	736	END DO
[2715]	737	DO ji = kideb, kiut ! Upper ice layer
[4045]	738	fc_i(ji,0) = - REAL( isnow(ji) ) * & ! interface flux if there is snow
[921]	739	( zkappa_i(ji,0) * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) &
[4045]	740	- REAL( 1 - isnow(ji) ) * ( zkappa_i(ji,0) * &
[921]	741	zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not
	742	END DO
[2715]	743	DO layer = 1, nlay_i - 1 ! Internal ice layers
[921]	744	DO ji = kideb, kiut
[2715]	745	fc_i(ji,layer) = - zkappa_i(ji,layer) * ( t_i_b(ji,layer+1) - t_i_b(ji,layer) )
	746	ii = MOD( npb(ji) - 1, jpi ) + 1
	747	ij = ( npb(ji) - 1 ) / jpi + 1
[921]	748	END DO
	749	END DO
[2715]	750	DO ji = kideb, kiut ! Bottom ice layers
	751	fc_i(ji,nlay_i) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
[921]	752	END DO
	753	ENDIF
[2715]	754	!
[921]	755	END SUBROUTINE lim_thd_dif
[825]	756
	757	#else
[2715]	758	!!----------------------------------------------------------------------
	759	!! Dummy Module No LIM-3 sea-ice model
	760	!!----------------------------------------------------------------------
[825]	761	CONTAINS
	762	SUBROUTINE lim_thd_dif ! Empty routine
	763	END SUBROUTINE lim_thd_dif
	764	#endif
[2528]	765	!!======================================================================
[921]	766	END MODULE limthd_dif

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

Download in other formats: