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

Context Navigation

source: trunk/NEMO/LIM_SRC_3/limthd_dif.F90 @ 1445

Last change on this file since 1445 was 1156, checked in by rblod, 16 years ago
Update Id and licence information, see ticket #210
Property svn:keywords set to `Id`
File size: 36.3 KB

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

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

Download in other formats: