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

Context Navigation

source: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90 @ 4506

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

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

Download in other formats: