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_ent.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90 @ 3319

Last change on this file since 3319 was 3294, checked in by rblod, 12 years ago
Merge of 3.4beta into the trunk
Property svn:keywords set to `Id`
File size: 27.0 KB

Rev	Line
[825]	1	MODULE limthd_ent
	2	!!======================================================================
	3	!! * MODULE limthd_ent *
	4	!! Redistribution of Enthalpy in the ice
	5	!! on the new vertical grid
	6	!! after vertical growth/decay
	7	!!======================================================================
[2715]	8	!! History : LIM ! 2003-05 (M. Vancoppenolle) Original code in 1D
	9	!! ! 2005-07 (M. Vancoppenolle) 3D version
	10	!! ! 2006-11 (X. Fettweis) Vectorized
	11	!! 3.0 ! 2008-03 (M. Vancoppenolle) Energy conservation and clean code
	12	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
	13	!!----------------------------------------------------------------------
[2528]	14	#if defined key_lim3
	15	!!----------------------------------------------------------------------
	16	!! 'key_lim3' LIM3 sea-ice model
	17	!!----------------------------------------------------------------------
[834]	18	!! lim_thd_ent : ice redistribution of enthalpy
[2528]	19	!!----------------------------------------------------------------------
[825]	20	USE par_oce ! ocean parameters
[2715]	21	USE dom_oce ! domain variables
	22	USE domain !
	23	USE phycst ! physical constants
	24	USE ice ! LIM variables
	25	USE par_ice ! LIM parameters
	26	USE thd_ice ! LIM thermodynamics
	27	USE limvar ! LIM variables
	28	USE in_out_manager ! I/O manager
	29	USE lib_mpp ! MPP library
[3294]	30	USE wrk_nemo ! work arrays
[825]	31
	32	IMPLICIT NONE
	33	PRIVATE
	34
[3294]	35	PUBLIC lim_thd_ent ! called by lim_thd
[825]	36
[2715]	37	REAL(wp) :: epsi20 = 1e-20_wp ! constant values
	38	REAL(wp) :: epsi13 = 1e-13_wp !
	39	REAL(wp) :: epsi10 = 1e-10_wp !
	40	REAL(wp) :: epsi06 = 1e-06_wp !
	41	REAL(wp) :: zzero = 0._wp !
	42	REAL(wp) :: zone = 1._wp !
	43
[825]	44	!!----------------------------------------------------------------------
[2715]	45	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	46	!! $Id$
[2528]	47	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	48	!!----------------------------------------------------------------------
	49	CONTAINS
[3294]	50
[2715]	51	SUBROUTINE lim_thd_ent( kideb, kiut, jl )
[825]	52	!!-------------------------------------------------------------------
	53	!! * ROUTINE lim_thd_ent *
	54	!!
	55	!! ** Purpose :
	56	!! This routine computes new vertical grids
	57	!! in the ice and in the snow, and consistently redistributes
	58	!! temperatures in the snow / ice.
	59	!! Redistribution is made so as to ensure to energy conservation
	60	!!
	61	!!
	62	!! ** Method : linear conservative remapping
	63	!!
[834]	64	!! ** Steps : 1) Grid
	65	!! 2) Switches
	66	!! 3) Snow redistribution
	67	!! 4) Ice enthalpy redistribution
	68	!! 5) Ice salinity, recover temperature
[825]	69	!!
[2715]	70	!! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005
	71	!!-------------------------------------------------------------------
	72	INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied
	73	INTEGER , INTENT(in) :: jl ! Thickness cateogry number
[825]	74
[2715]	75	INTEGER :: ji,jk ! dummy loop indices
	76	INTEGER :: zji, zjj , & ! dummy indices
[825]	77	ntop0 , & ! old layer top index
	78	nbot1 , & ! new layer bottom index
	79	ntop1 , & ! new layer top index
	80	limsum , & ! temporary loop index
	81	nlayi0,nlays0 , & ! old number of layers
	82	maxnbot0 , & ! old layer bottom index
	83	layer0, layer1 ! old/new layer indexes
	84
	85
	86	REAL(wp) :: &
	87	ztmelts , & ! ice melting point
	88	zqsnic , & ! enthalpy of snow ice layer
	89	zhsnow , & ! temporary snow thickness variable
	90	zswitch , & ! dummy switch argument
	91	zfac1 , & ! dummy factor
	92	zfac2 , & ! dummy factor
	93	ztform , & !: bottom formation temperature
	94	zaaa , & !: dummy factor
	95	zbbb , & !: dummy factor
	96	zccc , & !: dummy factor
	97	zdiscrim !: dummy factor
	98
[3294]	99	INTEGER, POINTER, DIMENSION(:) :: snswi ! snow switch
	100	INTEGER, POINTER, DIMENSION(:) :: nbot0 ! old layer bottom index
	101	INTEGER, POINTER, DIMENSION(:) :: icsuind ! ice surface index
	102	INTEGER, POINTER, DIMENSION(:) :: icsuswi ! ice surface switch
	103	INTEGER, POINTER, DIMENSION(:) :: icboind ! ice bottom index
	104	INTEGER, POINTER, DIMENSION(:) :: icboswi ! ice bottom switch
	105	INTEGER, POINTER, DIMENSION(:) :: snicind ! snow ice index
	106	INTEGER, POINTER, DIMENSION(:) :: snicswi ! snow ice switch
	107	INTEGER, POINTER, DIMENSION(:) :: snind ! snow index
[2715]	108	!
[3294]	109	REAL(wp), POINTER, DIMENSION(:) :: zh_i ! thickness of an ice layer
	110	REAL(wp), POINTER, DIMENSION(:) :: zh_s ! thickness of a snow layer
	111	REAL(wp), POINTER, DIMENSION(:) :: zqsnow ! enthalpy of the snow put in snow ice
	112	REAL(wp), POINTER, DIMENSION(:) :: zdeltah ! temporary variable
	113	REAL(wp), POINTER, DIMENSION(:) :: zqti_in, zqts_in
	114	REAL(wp), POINTER, DIMENSION(:) :: zqti_fin, zqts_fin
	115
	116	REAL(wp), POINTER, DIMENSION(:,:) :: zm0 ! old layer-system vertical cotes
	117	REAL(wp), POINTER, DIMENSION(:,:) :: qm0 ! old layer-system heat content
	118	REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! new snow system vertical cotes
	119	REAL(wp), POINTER, DIMENSION(:,:) :: z_i ! new ice system vertical cotes
	120	REAL(wp), POINTER, DIMENSION(:,:) :: zthick0 ! old ice thickness
	121	REAL(wp), POINTER, DIMENSION(:,:) :: zhl0 ! old and new layer thicknesses
	122	REAL(wp), POINTER, DIMENSION(:,:) :: zrl01
[2715]	123	!!-------------------------------------------------------------------
[825]	124
[3294]	125	CALL wrk_alloc( jpij, snswi, nbot0, icsuind, icsuswi, icboind, icboswi, snicind, snicswi, snind ) ! integer
	126	CALL wrk_alloc( jpij, zh_i, zh_s, zqsnow, zdeltah, zqti_in, zqts_in, zqti_fin, zqts_fin ) ! real
	127	CALL wrk_alloc( jpij,jkmax+4, zm0, qm0, z_s, z_i, zthick0, zhl0, kjstart = 0 )
	128	CALL wrk_alloc( jkmax+4,jkmax+4, zrl01, kistart = 0, kjstart = 0 )
[825]	129
[2715]	130	zthick0(:,:) = 0._wp
	131	zm0 (:,:) = 0._wp
	132	qm0 (:,:) = 0._wp
	133	zrl01 (:,:) = 0._wp
	134	zhl0 (:,:) = 0._wp
	135	z_i (:,:) = 0._wp
	136	z_s (:,:) = 0._wp
[825]	137
[921]	138	!
	139	!------------------------------------------------------------------------------\|
	140	! 1) Grid \|
	141	!------------------------------------------------------------------------------\|
[2715]	142	nlays0 = nlay_s
	143	nlayi0 = nlay_i
[825]	144
	145	DO ji = kideb, kiut
[2715]	146	zh_i(ji) = old_ht_i_b(ji) / nlay_i
	147	zh_s(ji) = old_ht_s_b(ji) / nlay_s
	148	END DO
[825]	149
[921]	150	!
	151	!------------------------------------------------------------------------------\|
	152	! 2) Switches \|
	153	!------------------------------------------------------------------------------\|
[825]	154	! 2.1 snind(ji), snswi(ji)
	155	! snow surface behaviour : computation of snind(ji)-snswi(ji)
	156	! snind(ji) : index which equals
	157	! 0 if snow is accumulating
	158	! 1 if 1st layer is melting
	159	! 2 if 2nd layer is melting ...
	160	DO ji = kideb, kiut
[2715]	161	snind (ji) = 0
	162	zdeltah(ji) = 0._wp
[825]	163	ENDDO !ji
	164
	165	DO jk = 1, nlays0
[921]	166	DO ji = kideb, kiut
[2715]	167	snind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-epsi20))) &
	168	+ snind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-epsi20))))
[921]	169	zdeltah(ji)= zdeltah(ji) + zh_s(ji)
	170	END DO ! ji
[2715]	171	END DO ! jk
[825]	172
	173	! snswi(ji) : switch which value equals 1 if snow melts
	174	! 0 if not
	175	DO ji = kideb, kiut
[2715]	176	snswi(ji) = MAX(0,INT(-dh_s_tot(ji)/MAX(epsi20,ABS(dh_s_tot(ji)))))
	177	END DO ! ji
[921]	178
[825]	179	! 2.2 icsuind(ji), icsuswi(ji)
	180	! ice surface behaviour : computation of icsuind(ji)-icsuswi(ji)
	181	! icsuind(ji) : index which equals
	182	! 0 if nothing happens at the surface
	183	! 1 if first layer is melting
	184	! 2 if 2nd layer is reached by melt ...
	185	DO ji = kideb, kiut
[2715]	186	icsuind(ji) = 0
	187	zdeltah(ji) = 0._wp
	188	END DO !ji
[825]	189	DO jk = 1, nlayi0
[921]	190	DO ji = kideb, kiut
[2715]	191	icsuind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-epsi20))) &
	192	+ icsuind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-epsi20))))
[921]	193	zdeltah(ji) = zdeltah(ji) + zh_i(ji)
	194	END DO ! ji
[825]	195	ENDDO !jk
	196
	197	! icsuswi(ji) : switch which equals
	198	! 1 if ice melts at the surface
	199	! 0 if not
	200	DO ji = kideb, kiut
[2715]	201	icsuswi(ji) = MAX(0,INT(-dh_i_surf(ji)/MAX(epsi20 , ABS(dh_i_surf(ji)) ) ) )
[825]	202	ENDDO
	203
	204	! 2.3 icboind(ji), icboswi(ji)
	205	! ice bottom behaviour : computation of icboind(ji)-icboswi(ji)
	206	! icboind(ji) : index which equals
	207	! 0 if accretion is on the way
	208	! 1 if last layer has started to melt
	209	! 2 if penultiem layer is melting ... and so on
	210	! N+1 if all layers melt and that snow transforms into ice
	211	DO ji = kideb, kiut
[2715]	212	icboind(ji) = 0
	213	zdeltah(ji) = 0._wp
	214	END DO
[825]	215	DO jk = nlayi0, 1, -1
[921]	216	DO ji = kideb, kiut
[2715]	217	icboind(ji) = (nlayi0+1-jk) * INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-epsi20))) &
	218	& + icboind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-epsi20))))
[921]	219	zdeltah(ji) = zdeltah(ji) + zh_i(ji)
	220	END DO
[2715]	221	END DO
[825]	222
	223	DO ji = kideb, kiut
	224	! case of total ablation with remaining snow
[2715]	225	IF ( ( ht_i_b(ji) .GT. epsi20 ) .AND. &
	226	( ht_i_b(ji) - dh_snowice(ji) .LT. epsi20 ) ) icboind(ji) = nlay_i + 1
[825]	227	END DO
	228
	229	! icboswi(ji) : switch which equals
	230	! 1 if ice accretion is on the way
	231	! 0 if ablation is on the way
	232	DO ji = kideb, kiut
[2715]	233	icboswi(ji) = MAX(0,INT(dh_i_bott(ji) / MAX(epsi20,ABS(dh_i_bott(ji)))))
	234	END DO
[825]	235
	236	! 2.4 snicind(ji), snicswi(ji)
	237	! snow ice formation : calcul de snicind(ji)-snicswi(ji)
	238	! snicind(ji) : index which equals
	239	! 0 if no snow-ice forms
	240	! 1 if last layer of snow has started to melt
	241	! 2 if penultiem layer ...
	242	DO ji = kideb, kiut
[2715]	243	snicind(ji) = 0
	244	zdeltah(ji) = 0._wp
	245	END DO
[825]	246	DO jk = nlays0, 1, -1
[921]	247	DO ji = kideb, kiut
	248	snicind(ji) = (nlays0+1-jk) &
[2715]	249	* INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-epsi20))) + snicind(ji) &
	250	* (1 - INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-epsi20))))
[921]	251	zdeltah(ji) = zdeltah(ji) + zh_s(ji)
	252	END DO
[2715]	253	END DO
[825]	254
	255	! snicswi(ji) : switch which equals
	256	! 1 if snow-ice forms
	257	! 0 if not
	258	DO ji = kideb, kiut
[2715]	259	snicswi(ji) = MAX(0,INT(dh_snowice(ji)/MAX(epsi20,ABS(dh_snowice(ji)))))
[825]	260	ENDDO
	261
[921]	262	!
	263	!------------------------------------------------------------------------------\|
	264	! 3) Snow redistribution \|
	265	!------------------------------------------------------------------------------\|
	266	!
[825]	267	!-------------
	268	! Old profile
	269	!-------------
	270
	271	! by 'old', it is meant that layers coming from accretion are included,
	272	! and that interfacial layers which were partly melted are reduced
	273
	274	! indexes of the vectors
	275	!------------------------
[2715]	276	ntop0 = 1
	277	maxnbot0 = 0
[825]	278
	279	DO ji = kideb, kiut
[2715]	280	nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1. - snicind(ji) ) * snicswi(ji)
[921]	281	! cotes of the top of the layers
[2715]	282	zm0(ji,0) = 0._wp
	283	maxnbot0 = MAX ( maxnbot0 , nbot0(ji) )
	284	END DO
	285	IF( lk_mpp ) CALL mpp_max( maxnbot0, kcom=ncomm_ice )
[825]	286
	287	DO jk = 1, maxnbot0
[921]	288	DO ji = kideb, kiut
	289	!change
[2715]	290	limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 )
	291	limsum = MIN( limsum , nlay_s )
	292	zm0(ji,jk) = dh_s_tot(ji) + zh_s(ji) * limsum
[921]	293	END DO
[2715]	294	END DO
[825]	295
	296	DO ji = kideb, kiut
[2715]	297	zm0(ji,nbot0(ji)) = dh_s_tot(ji) - snicswi(ji) * dh_snowice(ji) + zh_s(ji) * nlays0
	298	zm0(ji,1) = dh_s_tot(ji) * (1 -snswi(ji) ) + snswi(ji) * zm0(ji,1)
	299	END DO
[825]	300
	301	DO jk = ntop0, maxnbot0
[921]	302	DO ji = kideb, kiut
[2715]	303	zthick0(ji,jk) = zm0(ji,jk) - zm0(ji,jk-1) ! layer thickness
[921]	304	END DO
[2715]	305	END DO
[825]	306
[2715]	307	zqts_in(:) = 0._wp
[921]	308
[2715]	309	DO ji = kideb, kiut ! layer heat content
	310	qm0 (ji,1) = rhosn * ( cpic * ( rtt - ( 1. - snswi(ji) ) * tatm_ice_1d(ji) &
	311	& - snswi(ji) * t_s_b (ji,1) ) &
	312	& + lfus ) * zthick0(ji,1)
	313	zqts_in(ji) = zqts_in(ji) + qm0(ji,1)
	314	END DO
[825]	315
	316	DO jk = 2, maxnbot0
[921]	317	DO ji = kideb, kiut
[2715]	318	limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 )
[921]	319	limsum = MIN( limsum , nlay_s )
[2715]	320	qm0(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,limsum) ) + lfus ) * zthick0(ji,jk)
	321	zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, epsi20 - ht_s_b(ji) ) )
[921]	322	zqts_in(ji) = zqts_in(ji) + ( 1. - snswi(ji) ) * qm0(ji,jk) * zswitch
	323	END DO ! jk
[2715]	324	END DO ! ji
[825]	325
	326	!------------------------------------------------
	327	! Energy given by the snow in snow-ice formation
	328	!------------------------------------------------
	329	! zqsnow, enthalpy of the flooded snow
	330	DO ji = kideb, kiut
[2715]	331	zqsnow (ji) = rhosn * lfus
	332	zdeltah(ji) = 0._wp
	333	END DO
[825]	334
	335	DO jk = nlays0, 1, -1
[921]	336	DO ji = kideb, kiut
[2715]	337	zhsnow = MAX( 0._wp , dh_snowice(ji)-zdeltah(ji) )
	338	zqsnow (ji) = zqsnow (ji) + rhosncpic(rtt-t_s_b(ji,jk))
[921]	339	zdeltah(ji) = zdeltah(ji) + zh_s(ji)
	340	END DO
[2715]	341	END DO
[825]	342
	343	DO ji = kideb, kiut
	344	zqsnow(ji) = zqsnow(ji) * dh_snowice(ji)
	345	END DO
	346
	347	!------------------
[834]	348	! new snow profile
[825]	349	!------------------
	350
	351	!--------------
	352	! Vector index
	353	!--------------
[2715]	354	ntop1 = 1
	355	nbot1 = nlay_s
[825]	356
	357	!-------------------
	358	! Layer coordinates
	359	!-------------------
	360	DO ji = kideb, kiut
	361	zh_s(ji) = ht_s_b(ji) / nlay_s
[2715]	362	z_s(ji,0) = 0._wp
[825]	363	ENDDO
	364
	365	DO jk = 1, nlay_s
[921]	366	DO ji = kideb, kiut
	367	z_s(ji,jk) = zh_s(ji) * jk
	368	END DO
[2715]	369	END DO
[825]	370
	371	!-----------------
	372	! Layer thickness
	373	!-----------------
	374	DO layer0 = ntop0, maxnbot0
[921]	375	DO ji = kideb, kiut
	376	zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1)
	377	END DO
[2715]	378	END DO
[825]	379
	380	DO layer1 = ntop1, nbot1
[921]	381	DO ji = kideb, kiut
[2715]	382	q_s_b(ji,layer1) = 0._wp
[921]	383	END DO
[2715]	384	END DO
[825]	385
	386	!----------------
	387	! Weight factors
	388	!----------------
	389	DO layer0 = ntop0, maxnbot0
[921]	390	DO layer1 = ntop1, nbot1
	391	DO ji = kideb, kiut
[2715]	392	zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1)) &
	393	& - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10))
	394	q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0) &
	395	& * MAX(0.0,SIGN(1.0,nbot0(ji)-layer0+epsi20))
[921]	396	END DO
	397	END DO
[2715]	398	END DO
[825]	399
	400	! Heat conservation
[2715]	401	zqts_fin(:) = 0._wp
[825]	402	DO jk = 1, nlay_s
	403	DO ji = kideb, kiut
	404	zqts_fin(ji) = zqts_fin(ji) + q_s_b(ji,jk)
	405	END DO
	406	END DO
	407
	408	IF ( con_i ) THEN
[921]	409	DO ji = kideb, kiut
	410	IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) / rdt_ice .GT. 1.0e-6 ) THEN
	411	zji = MOD( npb(ji) - 1, jpi ) + 1
	412	zjj = ( npb(ji) - 1 ) / jpi + 1
	413	WRITE(numout,*) ' violation of heat conservation : ', &
	414	ABS ( zqts_in(ji) - zqts_fin(ji) ) / rdt_ice
	415	WRITE(numout,*) ' ji, jj : ', zji, zjj
	416	WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji)
	417	WRITE(numout,*) ' zqts_in : ', zqts_in(ji) / rdt_ice
	418	WRITE(numout,*) ' zqts_fin : ', zqts_fin(ji) / rdt_ice
	419	WRITE(numout,*) ' dh_snowice : ', dh_snowice(ji)
	420	WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji)
	421	WRITE(numout,*) ' snswi : ', snswi(ji)
	422	ENDIF
	423	END DO
[825]	424	ENDIF
	425
	426	!---------------------
	427	! Recover heat content
	428	!---------------------
[1724]	429	DO jk = 1, nlay_s
[825]	430	DO ji = kideb, kiut
[2715]	431	q_s_b(ji,jk) = q_s_b(ji,jk) / MAX( zh_s(ji) , epsi20 )
[825]	432	END DO !ji
[2715]	433	END DO !jk
[825]	434
	435	!---------------------
	436	! Recover temperature
	437	!---------------------
	438	zfac1 = 1. / ( rhosn * cpic )
	439	zfac2 = lfus / cpic
	440	DO jk = 1, nlay_s
[921]	441	DO ji = kideb, kiut
[2715]	442	zswitch = MAX ( 0.0 , SIGN ( 1.0, epsi20 - ht_s_b(ji) ) )
	443	t_s_b(ji,jk) = rtt + ( 1.0 - zswitch ) * ( - zfac1 * q_s_b(ji,jk) + zfac2 )
[921]	444	END DO
[2715]	445	END DO
[921]	446	!
	447	!------------------------------------------------------------------------------\|
	448	! 4) Ice redistribution \|
	449	!------------------------------------------------------------------------------\|
	450	!
[825]	451	!-------------
	452	! OLD PROFILE
	453	!-------------
	454
	455	!----------------
	456	! Vector indexes
	457	!----------------
[2715]	458	ntop0 = 1
	459	maxnbot0 = 0
[825]	460
	461	DO ji = kideb, kiut
[921]	462	! reference number of the bottommost layer
[2715]	463	nbot0(ji) = MAX( 1 , MIN( nlayi0 + ( 1 - icboind(ji) ) + &
	464	& ( 1 - icsuind(ji) ) * icsuswi(ji) + snicswi(ji) , nlay_i + 2 ) )
[825]	465	! maximum reference number of the bottommost layer over all domain
[2715]	466	maxnbot0 = MAX( maxnbot0 , nbot0(ji) )
	467	END DO
[825]	468
	469	!-------------------------
	470	! Cotes of old ice layers
	471	!-------------------------
[2777]	472	zm0(:,0) = 0._wp
[921]	473
[825]	474	DO jk = 1, maxnbot0
	475	DO ji = kideb, kiut
	476	! jk goes from 1 to nbot0
	477	! the ice layer number goes from 1 to nlay_i
	478	! limsum is the real ice layer number corresponding to present jk
[834]	479	limsum = ( (icsuswi(ji)*(icsuind(ji)+jk-1) + &
[921]	480	(1-icsuswi(ji))jk))(1-snicswi(ji)) + (jk-1)*snicswi(ji)
[825]	481	zm0(ji,jk)= icsuswi(ji)dh_i_surf(ji) + snicswi(ji)dh_snowice(ji) &
[921]	482	+ limsum * zh_i(ji)
[825]	483	END DO
[2715]	484	END DO
[825]	485
	486	DO ji = kideb, kiut
	487	zm0(ji,nbot0(ji)) = icsuswi(ji)dh_i_surf(ji) + snicswi(ji)dh_snowice(ji) + dh_i_bott(ji) &
[921]	488	+ zh_i(ji) * nlayi0
[825]	489	zm0(ji,1) = snicswi(ji)dh_snowice(ji) + (1-snicswi(ji))zm0(ji,1)
[2715]	490	END DO
[825]	491
	492	!-----------------------------
	493	! Thickness of old ice layers
	494	!-----------------------------
	495	DO jk = ntop0, maxnbot0
[921]	496	DO ji = kideb, kiut
	497	zthick0(ji,jk) = zm0(ji,jk) - zm0(ji,jk-1)
	498	END DO
[2715]	499	END DO
[825]	500
	501	!---------------------------
	502	! Inner layers heat content
	503	!---------------------------
	504	qm0(:,:) = 0.0
	505	zqti_in(:) = 0.0
	506
	507	DO jk = ntop0, maxnbot0
	508	DO ji = kideb, kiut
	509	limsum = MAX(1,MIN(snicswi(ji)(jk-1) + icsuswi(ji)(jk-1+icsuind(ji)) + &
[921]	510	(1-icsuswi(ji))(1-snicswi(ji))jk,nlay_i))
[825]	511	ztmelts = -tmut * s_i_b(ji,limsum) + rtt
	512	qm0(ji,jk) = rhoic * ( cpic * (ztmelts-t_i_b(ji,limsum)) + lfus * ( 1.0-(ztmelts-rtt)/ &
[2715]	513	MIN((t_i_b(ji,limsum)-rtt),-epsi20) ) - rcp*(ztmelts-rtt) ) &
[921]	514	* zthick0(ji,jk)
[825]	515	END DO
[2715]	516	END DO
[825]	517
	518	!----------------------------
	519	! Bottom layers heat content
	520	!----------------------------
	521	DO ji = kideb, kiut
[2715]	522	ztmelts = ( 1.0 - icboswi(ji) ) * (-tmut * s_i_b (ji,nlayi0) ) & ! case of melting ice
	523	& + icboswi(ji) * (-tmut * s_i_new(ji) ) & ! case of forming ice
	524	& + rtt ! in Kelvin
[825]	525
[921]	526	! bottom formation temperature
	527	ztform = t_i_b(ji,nlay_i)
	528	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) ztform = t_bo_b(ji)
[2715]	529	qm0(ji,nbot0(ji)) = ( 1.0 - icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice
	530	& + icboswi(ji) * rhoic * ( cpic*(ztmelts-ztform) & ! case of forming ice
	531	+ lfus *( 1.0-(ztmelts-rtt) / MIN ( (ztform-rtt) , - epsi10 ) ) &
	532	- rcp(ztmelts-rtt) ) zthick0(ji,nbot0(ji) )
	533	END DO
[825]	534
	535	!-----------------------------
	536	! Snow ice layer heat content
	537	!-----------------------------
	538	DO ji = kideb, kiut
	539	! energy of the flooding seawater
	540	zqsnic = rau0 * rcp * ( rtt - t_bo_b(ji) ) * dh_snowice(ji) * &
[921]	541	(rhoic - rhosn) / rhoic * snicswi(ji) ! generally positive
[825]	542	! Heat conservation diagnostic
	543	qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic
	544
	545	qldif_1d(ji) = qldif_1d(ji) + zqsnic * a_i_b(ji)
	546
	547	! enthalpy of the newly formed snow-ice layer
	548	! = enthalpy of snow + enthalpy of frozen water
	549	zqsnic = zqsnow(ji) + zqsnic
	550	qm0(ji,1) = snicswi(ji) * zqsnic + ( 1 - snicswi(ji) ) * qm0(ji,1)
	551
[2715]	552	END DO ! ji
[825]	553
	554	DO jk = ntop0, maxnbot0
[921]	555	DO ji = kideb, kiut
	556	! Heat conservation
[2715]	557	zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06+epsi20) ) &
	558	& * MAX( 0.0 , SIGN( 1. , nbot0(ji) - jk + epsi20 ) )
[921]	559	END DO
[2715]	560	END DO
[825]	561
	562	!-------------
	563	! NEW PROFILE
	564	!-------------
	565
	566	!---------------
	567	! Vectors index
	568	!---------------
[2715]	569	ntop1 = 1
	570	nbot1 = nlay_i
[825]	571
	572	!------------------
	573	! Layers thickness
	574	!------------------
	575	DO ji = kideb, kiut
[2715]	576	zh_i(ji) = ht_i_b(ji) / nlay_i
[825]	577	ENDDO
	578
	579	!-------------
	580	! Layer cotes
	581	!-------------
[2715]	582	z_i(:,0) = 0._wp
[825]	583	DO jk = 1, nlay_i
[921]	584	DO ji = kideb, kiut
	585	z_i(ji,jk) = zh_i(ji) * jk
	586	END DO
[2715]	587	END DO
[825]	588
	589	!--thicknesses of the layers
	590	DO layer0 = ntop0, maxnbot0
[921]	591	DO ji = kideb, kiut
[2715]	592	zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1) ! thicknesses of the layers
[921]	593	END DO
[2715]	594	END DO
[825]	595
	596	!------------------------
	597	! Weights for relayering
	598	!------------------------
[2715]	599	q_i_b(:,:) = 0._wp
[825]	600	DO layer0 = ntop0, maxnbot0
[921]	601	DO layer1 = ntop1, nbot1
	602	DO ji = kideb, kiut
	603	zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) &
	604	- MAX(zm0(ji,layer0-1), z_i(ji,layer1-1)))/MAX(zhl0(ji,layer0),epsi10))
	605	q_i_b(ji,layer1) = q_i_b(ji,layer1) &
	606	+ zrl01(layer1,layer0)*qm0(ji,layer0) &
[2715]	607	* MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06+epsi20)) &
	608	* MAX(0.0,SIGN(1.0,nbot0(ji)-layer0+epsi20))
[921]	609	END DO
	610	END DO
[2715]	611	END DO
[825]	612
	613	!-------------------------
	614	! Heat conservation check
	615	!-------------------------
[2715]	616	zqti_fin(:) = 0._wp
[825]	617	DO jk = 1, nlay_i
	618	DO ji = kideb, kiut
	619	zqti_fin(ji) = zqti_fin(ji) + q_i_b(ji,jk)
	620	END DO
	621	END DO
[921]	622	!
[825]	623	DO ji = kideb, kiut
	624	IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) / rdt_ice .GT. 1.0e-6 ) THEN
	625	zji = MOD( npb(ji) - 1, jpi ) + 1
	626	zjj = ( npb(ji) - 1 ) / jpi + 1
[2715]	627	WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) / rdt_ice
[825]	628	WRITE(numout,*) ' ji, jj : ', zji, zjj
	629	WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji)
	630	WRITE(numout,*) ' zqti_in : ', zqti_in(ji) / rdt_ice
	631	WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) / rdt_ice
	632	WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji)
	633	WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji)
	634	WRITE(numout,*) ' dh_snowice:', dh_snowice(ji)
	635	WRITE(numout,*) ' icsuswi : ', icsuswi(ji)
	636	WRITE(numout,*) ' icboswi : ', icboswi(ji)
	637	WRITE(numout,*) ' snicswi : ', snicswi(ji)
	638	ENDIF
	639	END DO
	640
	641	!----------------------
	642	! Recover heat content
	643	!----------------------
	644	DO jk = 1, nlay_i
	645	DO ji = kideb, kiut
[2715]	646	q_i_b(ji,jk) = q_i_b(ji,jk) / MAX( zh_i(ji) , epsi20 )
[825]	647	END DO !ji
[2715]	648	END DO !jk
[825]	649
	650	! Heat conservation
	651	zqti_fin(:) = 0.0
	652	DO jk = 1, nlay_i
	653	DO ji = kideb, kiut
	654	zqti_fin(ji) = zqti_fin(ji) + q_i_b(ji,jk) * zh_i(ji)
	655	END DO
	656	END DO
[834]	657
[921]	658	!
	659	!------------------------------------------------------------------------------\|
	660	! 5) Update salinity and recover temperature \|
	661	!------------------------------------------------------------------------------\|
	662	!
[834]	663	! Update salinity (basal entrapment, snow ice formation)
[825]	664	DO ji = kideb, kiut
[2715]	665	sm_i_b(ji) = sm_i_b(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji)
[825]	666	END DO !ji
	667
	668	! Recover temperature
	669	DO jk = 1, nlay_i
	670	DO ji = kideb, kiut
	671	ztmelts = -tmut*s_i_b(ji,jk) + rtt
	672	!Conversion q(S,T) -> T (second order equation)
	673	zaaa = cpic
[2715]	674	zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus
[825]	675	zccc = lfus * ( ztmelts - rtt )
	676	zdiscrim = SQRT( MAX(zbbbzbbb - 4.0zaaa*zccc,0.0) )
[2715]	677	t_i_b(ji,jk) = rtt - ( zbbb + zdiscrim ) / ( 2.0 *zaaa )
[825]	678	END DO !ji
	679
	680	END DO !jk
[2715]	681	!
[3294]	682	CALL wrk_dealloc( jpij, snswi, nbot0, icsuind, icsuswi, icboind, icboswi, snicind, snicswi, snind ) ! integer
	683	CALL wrk_dealloc( jpij, zh_i, zh_s, zqsnow, zdeltah, zqti_in, zqts_in, zqti_fin, zqts_fin ) ! real
	684	CALL wrk_dealloc( jpij,jkmax+4, zm0, qm0, z_s, z_i, zthick0, zhl0, kjstart = 0 )
	685	CALL wrk_dealloc( jkmax+4,jkmax+4, zrl01, kistart = 0, kjstart = 0 )
[2715]	686	!
[921]	687	END SUBROUTINE lim_thd_ent
[825]	688
	689	#else
[2715]	690	!!----------------------------------------------------------------------
	691	!! Default option NO LIM3 sea-ice model
	692	!!----------------------------------------------------------------------
[825]	693	CONTAINS
	694	SUBROUTINE lim_thd_ent ! Empty routine
	695	END SUBROUTINE lim_thd_ent
	696	#endif
[2715]	697
	698	!!======================================================================
[921]	699	END MODULE limthd_ent

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

Download in other formats: