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 @ 4444

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

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

Download in other formats: