Context Navigation

sisvat_zsn.F @ 4666

Last change on this file since 4666 was 3900, checked in by evignon, 3 years ago
Commission de la nouvelle interface LMDZ-SISVAT la prochaine commission consistera a supprimer l'ancien repertoire sisvat et a faire un peu de nettoyage.
File size: 33.6 KB

Rev	Line
[3792]	1
	2
	3	subroutine SISVAT_zSn
	4
	5	C +------------------------------------------------------------------------+
	6	C \| MAR SISVAT_zSn 12-07-2019 MAR \|
	7	C \| SubRoutine SISVAT_zSn manages the Snow Pack vertical Discretization \|
	8	C \| \|
	9	C +------------------------------------------------------------------------+
	10	C \| \|
	11	C \| PARAMETERS: knonv: Total Number of columns = \|
	12	C \| ^^^^^^^^^^ = Total Number of continental grid boxes \|
	13	C \| X Number of Mosaic Cell per grid box \|
	14	C \| \|
	15	C \| INPUT / NLaysv = New Snow Layer Switch \|
	16	C \| OUTPUT: isnoSV = total Nb of Ice/Snow Layers \|
	17	C \| ^^^^^^ ispiSV = 0,...,nsno: Uppermost Superimposed Ice Layer \|
	18	C \| iiceSV = total Nb of Ice Layers \|
	19	C \| istoSV = 0,...,5 : Snow History (see istdSV data) \|
	20	C \| \|
	21	C \| INPUT / TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)\|
	22	C \| OUTPUT: & Snow Temperatures (layers 1,2,...,nsno) [K] \|
	23	C \| ^^^^^^ ro__SV : Soil/Snow Volumic Mass [kg/m3] \|
	24	C \| eta_SV : Soil/Snow Water Content [m3/m3] \|
	25	C \| dzsnSV : Snow Layer Thickness [m] \|
	26	C \| G1snSV : Dendricity (<0) or Sphericity (>0) of Snow Layer \|
	27	C \| G2snSV : Sphericity (>0) or Size of Snow Layer \|
	28	C \| agsnSV : Snow Age [day] \|
	29	C \| \|
	30	C \| METHOD: 1) Agregate the thinest Snow Layer \|
	31	C \| ^^^^^^ if a new Snow Layer has been precipitated (NLaysv = 1) \|
	32	C \| 2) Divide a too thick Snow Layer except \|
	33	C \| if the maximum Number of Layer is reached \|
	34	C \| in this case forces NLay_s = 1 \|
	35	C \| 3) Agregate the thinest Snow Layer \|
	36	C \| in order to divide a too thick Snow Layer \|
	37	C \| at next Time Step when NLay_s = 1 \|
	38	C \| \|
	39	C +------------------------------------------------------------------------+
	40
	41
	42
	43
	44	C +--Global Variables
	45	C + ================
	46
	47
	48	use VARphy
	49	use VAR_SV
	50	use VARdSV
	51	use VAR0SV
	52	use VARxSV
	53	use VARySV
[3900]	54	use surface_data, only: ok_zsn_ii
[3792]	55
	56	IMPLICIT NONE
	57
	58
	59	C +--Internal Variables
	60	C + ==================
	61
	62	integer ikl ,isn ,i !
	63
	64
	65
	66
	67
	68
	69	integer NLay_s(knonv) ! Split Snow Layer Switch
	70	integer isagr1(knonv) ! 1st Layer History
	71	integer isagr2(knonv) ! 2nd Layer History
	72	integer LstLay ! 0 ====> isnoSV = 1
	73	integer isno_n ! Snow Normal.Profile
	74	integer iice_n ! Ice Normal.Profile
	75	integer iiceOK ! Ice Switch
	76	integer icemix ! 0 ====> Agregated Snow+Ice=Snow
	77	C + ! 1 Ice
	78	integer isn1 (knonv) ! 1st layer to stagger
	79	real staggr ! stagger Switch
	80
	81	real WEagre(knonv) ! Snow Water Equivalent Thickness
	82	real dzthin(knonv) ! Thickness of the thinest layer
	83	real OKthin ! Swich ON a new thinest layer
	84	real dz_dif ! difference from ideal discret.
	85	real thickL ! Thick Layer Indicator
	86	real OK_ICE ! Swich ON uppermost Ice Layer
	87
	88	real Agrege(knonv) ! 1. when Agregation constrained
	89	real dzepsi ! Min Single Snw Layer Thickness
	90	real dzxmin ! Min Acceptable Layer Thickness
	91	real dz_min ! Min Layer Thickness
	92	real dz_max ! Max Layer Thickness
	93	real dzagr1(knonv) ! 1st Layer Thickness
	94	real dzagr2(knonv) ! 2nd Layer Thickness
	95	real T_agr1(knonv) ! 1st Layer Temperature
	96	real T_agr2(knonv) ! 2nd Layer Temperature
	97	real roagr1(knonv) ! 1st Layer Density
	98	real roagr2(knonv) ! 2nd Layer Density
	99	real etagr1(knonv) ! 1st Layer Water Content
	100	real etagr2(knonv) ! 2nd Layer Water Content
	101	real G1agr1(knonv) ! 1st Layer Dendricity/Spher.
	102	real G1agr2(knonv) ! 2nd Layer Dendricity/Spher.
	103	real G2agr1(knonv) ! 1st Layer Sphericity/Size
	104	real G2agr2(knonv) ! 2nd Layer Sphericity/Size
	105	real agagr1(knonv) ! 1st Layer Age
	106	real agagr2(knonv) ! 2nd Layer Age
	107
	108
	109	C +--DATA
	110	C + ====
	111
	112	data icemix / 0 / ! 0 ====> Agregated Snow+Ice=Snow
	113	data dzepsi / 0.0020/ ! Min single Layer Thickness
	114	data dzxmin / 0.0025/ ! Min accept.Layer Thickness
	115	c #EU data dz_min / 0.0050/ ! Min Local Layer Thickness < SMn
	116	data dz_min / 0.0040/ ! Min Local Layer Thickness < SMn
	117	data dz_max / 0.0300/ ! Min Gener. Layer Thickness
	118	C + CAUTION: dz_max > dz_min*2 is required ! Otherwise re-agregation is
	119	C + ! activated after splitting
	120
	121
	122
	123
	124
	125	C +--Constrains Agregation of too thin Layers
	126	C + =================================================
	127
	128	C +--Search the thinest non-zero Layer
	129	C + ----------------------------------
	130
	131	DO ikl=1,knonv
	132	if(isnoSV(ikl)<=2) dz_min=max(0.0050,dz_min)
	133
	134	dzepsi=0.0015
	135	if(ro__SV(ikl,isnoSV(ikl))>920) dzepsi=0.0020
	136
	137	dzthin(ikl) = 0. ! Arbitrary unrealistic
	138	END DO ! Layer Thickness
	139	cXF
	140	DO ikl=1,knonv
	141	DO isn=1,isnoSV(ikl)-3 ! no agregation of 3 first snowlayers
	142	! XF 04/07/2019
	143
	144	isno_n = isnoSV(ikl)-isn+1 ! Snow Normal.Profile
	145	iice_n = iiceSV(ikl)-isn ! Ice Normal.Profile
	146	iiceOK = min(1,max(0,iice_n +1)) ! Ice Switch
	147	! #vz dz_ref(isn) = !
	148	! #vz. dz_min ((1-iiceOK)isno_n*isno_n ! Theoretical Profile
	149	! #vz. + iiceOK * 2**iice_n) !
	150	! #vz. /max(1,isnoSV(ikl)) !
	151	dz_dif = max(zero, ! Actual Profile
	152	. dz_min !
	153	. ((1-iiceOK)isno_n*isno_n ! Theoretical Profile
	154	. + iiceOK 2. *iice_n) !
	155	. - dzsnSV(ikl, isn) ) ! Actual Profile
	156	! #vz dzwdif(isn) = dz_dif !
	157	OKthin = max(zero, !
	158	. sign(unun, !
	159	. dz_dif-dzthin(ikl))) ! 1.=> New thinest Lay.
	160	. * max(0, ! 1 => .le. isnoSV
	161	. min(1, ! 1 => isn is in the
	162	. isnoSV(ikl)-isn +1 )) ! Snow Pack
	163	. * min(unun, !
	164	!
	165	! 1st additional Condition to accept OKthin
	166	. max(zero, ! combination
	167	. sign(unun,G1snSV(ikl, isn ) ! G1 with same
	168	. *G1snSV(ikl,max(1,isn-1)))) ! sign => OK
	169	!
	170	! 2nd additional Condition to accept OKthin
	171	. + max(zero, ! G1>0
	172	. sign(unun,G1snSV(ikl, isn ))) ! =>OK
	173	!
	174	! 3rd additional Condition to accept OKthin
	175	. + max(zero, ! dz too small
	176	. sign(unun,dzxmin ! =>OK
	177	. -dzsnSV(ikl, isn ))))!
	178
	179	i_thin(ikl) = (1. - OKthin) * i_thin(ikl) ! Update thinest Lay.
	180	. + OKthin * isn ! Index
	181	dzthin(ikl) = (1. - OKthin) * dzthin(ikl) !
	182	. + OKthin * dz_dif !
	183	END DO
	184	END DO
	185
	186
	187
	188	DO ikl=1,knonv
	189	DO isn=1,isnoSV(ikl)
	190	OKthin = max(zero, !
	191	. sign(unun, !
	192	. dz_min !
	193	. -dzsnSV(ikl,isn))) !
	194	. * max(zero, ! ON if dz > 0
	195	. sign(unun, !
	196	. dzsnSV(ikl,isn)-epsi)) !
	197	. *min(1,max(0, ! Multiple Snow Lay.
	198	. min (1, ! Switch = 1
	199	. isnoSV(ikl) ! if isno > iice + 1
	200	. -iiceSV(ikl)-1)) !
	201	C + !
	202	. +int(max(zero, !
	203	. sign(unun, !
	204	. dzepsi ! Minimum accepted for
	205	. -dzsnSV(ikl,isn)))) ! 1 Snow Layer over Ice
	206	. *int(max(zero, ! ON if dz > 0
	207	. sign(unun, !
	208	. dzsnSV(ikl,isn)-epsi)))!
	209	. *(1 -min (abs(isnoSV(ikl) ! Switch = 1
	210	. -iiceSV(ikl)-1),1)) ! if isno = iice + 1
	211	C + !
	212	. +max(0, ! Ice
	213	. min (1, ! Switch
	214	. iiceSV(ikl)+1-isn))) !
	215	. *min(unun, !
	216	. max(zero, ! combination
	217	. sign(unun,G1snSV(ikl, isn ) ! G1>0 + G1<0
	218	. *G1snSV(ikl,max(1,isn-1)))) ! NO
	219	. + max(zero, !
	220	. sign(unun,G1snSV(ikl, isn ))) !
	221	. + max(zero, !
	222	. sign(unun,dzxmin !
	223	. -dzsnSV(ikl, isn ))))!
	224	i_thin(ikl) = (1. - OKthin) * i_thin(ikl) ! Update thinest Lay.
	225	. + OKthin * isn ! Index
	226	END DO
	227	END DO
	228
	229
	230
	231
	232	C + ***************
	233	call SISVAT_zCr
	234	C + ***************
	235
	236
	237	C +--Assign the 2 Layers to agregate
	238	C + -------------------------------
	239
	240	DO ikl=1,knonv
	241	isn = i_thin(ikl)
	242	if(LIndsv(ikl)>0) isn=min(nsno-1,isn) ! cXF
	243	isagr1(ikl) = istoSV(ikl,isn)
	244	isagr2(ikl) = istoSV(ikl,isn+LIndsv(ikl))
	245	dzagr1(ikl) = dzsnSV(ikl,isn)
	246	dzagr2(ikl) = dzsnSV(ikl,isn+LIndsv(ikl))
	247	T_agr1(ikl) = TsisSV(ikl,isn)
	248	T_agr2(ikl) = TsisSV(ikl,isn+LIndsv(ikl))
	249	roagr1(ikl) = ro__SV(ikl,isn)
	250	roagr2(ikl) = ro__SV(ikl,isn+LIndsv(ikl))
	251	etagr1(ikl) = eta_SV(ikl,isn)
	252	etagr2(ikl) = eta_SV(ikl,isn+LIndsv(ikl))
	253	G1agr1(ikl) = G1snSV(ikl,isn)
	254	G1agr2(ikl) = G1snSV(ikl,isn+LIndsv(ikl))
	255	G2agr1(ikl) = G2snSV(ikl,isn)
	256	G2agr2(ikl) = G2snSV(ikl,isn+LIndsv(ikl))
	257	agagr1(ikl) = agsnSV(ikl,isn)
	258	agagr2(ikl) = agsnSV(ikl,isn+LIndsv(ikl))
	259	LstLay = min(1,max( 0,isnoSV(ikl) -1)) ! 0 if single Layer
	260	isnoSV(ikl) = isnoSV(ikl) ! decrement isnoSV
	261	. -(1-LstLay)* max(zero, ! if downmost Layer
	262	. sign(unun,eps_21 ! < 1.e-21 m
	263	. -dzsnSV(ikl,1))) !
	264	isnoSV(ikl) = max( 0, isnoSV(ikl) ) !
	265	Agrege(ikl) = max(zero, !
	266	. sign(unun,dz_min ! No Agregation
	267	. -dzagr1(ikl) )) ! if too thick Layer
	268	. *LstLay ! if a single Layer
	269	. * min( max(0 ,isnoSV(ikl)+1 ! if Agregation
	270	. -i_thin(ikl) ! with a Layer
	271	. -LIndsv(ikl) ),1) ! above the Pack
	272
	273	WEagre(ikl) = 0.
	274	END DO
	275
	276
	277	DO ikl=1,knonv
	278	DO isn=1,isnoSV(ikl)
	279	WEagre(ikl) = WEagre(ikl) + ro__SV(ikl,isn)*dzsnSV(ikl,isn)
	280	. *min(1,max(0,i_thin(ikl)+1-isn))
	281	ENDDO
	282	ENDDO
	283
	284
	285	C +--Agregates
	286	C + ---------
	287
	288	C + ***************
	289	call SISVAT_zAg
	290	. (isagr1,isagr2,WEagre
	291	. ,dzagr1,dzagr2,T_agr1,T_agr2
	292	. ,roagr1,roagr2,etagr1,etagr2
	293	. ,G1agr1,G1agr2,G2agr1,G2agr2
	294	. ,agagr1,agagr2,Agrege
	295	. )
	296	C + ***************
	297
	298
	299	C +--Rearranges the Layers
	300	C + ---------------------
	301
	302	C +--New (agregated) Snow layer
	303	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	304	DO ikl=1,knonv
	305	isn = i_thin(ikl)
	306	isn = min(isn,isn+LIndsv(ikl))
	307	isnoSV(ikl) = max(0.,isnoSV(ikl) -Agrege(ikl))
	308	iiceSV(ikl) = iiceSV(ikl)
	309	. -max(0,sign(1,iiceSV(ikl) -isn +icemix))
	310	. *Agrege(ikl)
	311	. *max(0,sign(1,iiceSV(ikl) -1 ))
	312	istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn)
	313	. + Agrege(ikl) *isagr1(ikl)
	314	dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn)
	315	. + Agrege(ikl) *dzagr1(ikl)
	316	TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn)
	317	. + Agrege(ikl) *T_agr1(ikl)
	318	ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn)
	319	. + Agrege(ikl) *roagr1(ikl)
	320	eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn)
	321	. + Agrege(ikl) *etagr1(ikl)
	322	G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn)
	323	. + Agrege(ikl) *G1agr1(ikl)
	324	G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn)
	325	. + Agrege(ikl) *G2agr1(ikl)
	326	agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn)
	327	. + Agrege(ikl) *agagr1(ikl)
	328	END DO
	329
	330	C +--Above
	331	C + ^^^^^
	332	DO ikl=1,knonv
	333	isn1(ikl)=max(i_thin(ikl),i_thin(ikl)+LIndsv(ikl))
	334	END DO
	335	DO i= 1,nsno-1
	336	DO ikl=1,knonv
	337	staggr = min(1,max(0,i +1 -isn1(ikl) ))
	338	istoSV(ikl,i) = (1.-staggr )*istoSV(ikl,i )
	339	. + staggr((1.-Agrege(ikl))istoSV(ikl,i )
	340	. + Agrege(ikl) *istoSV(ikl,i+1))
	341	dzsnSV(ikl,i) = (1.-staggr )*dzsnSV(ikl,i )
	342	. + staggr((1.-Agrege(ikl))dzsnSV(ikl,i )
	343	. + Agrege(ikl) *dzsnSV(ikl,i+1))
	344	TsisSV(ikl,i) = (1.-staggr )*TsisSV(ikl,i )
	345	. + staggr((1.-Agrege(ikl))TsisSV(ikl,i )
	346	. + Agrege(ikl) *TsisSV(ikl,i+1))
	347	ro__SV(ikl,i) = (1.-staggr )*ro__SV(ikl,i )
	348	. + staggr((1.-Agrege(ikl))ro__SV(ikl,i )
	349	. + Agrege(ikl) *ro__SV(ikl,i+1))
	350	eta_SV(ikl,i) = (1.-staggr )*eta_SV(ikl,i )
	351	. + staggr((1.-Agrege(ikl))eta_SV(ikl,i )
	352	. + Agrege(ikl) *eta_SV(ikl,i+1))
	353	G1snSV(ikl,i) = (1.-staggr )*G1snSV(ikl,i )
	354	. + staggr((1.-Agrege(ikl))G1snSV(ikl,i )
	355	. + Agrege(ikl) *G1snSV(ikl,i+1))
	356	G2snSV(ikl,i) = (1.-staggr )*G2snSV(ikl,i )
	357	. + staggr((1.-Agrege(ikl))G2snSV(ikl,i )
	358	. + Agrege(ikl) *G2snSV(ikl,i+1))
	359	agsnSV(ikl,i) = (1.-staggr )*agsnSV(ikl,i )
	360	. + staggr((1.-Agrege(ikl))agsnSV(ikl,i )
	361	. + Agrege(ikl) *agsnSV(ikl,i+1))
	362	END DO
	363	END DO
	364
	365	DO ikl=1,knonv
	366	isn = min(isnoSV(ikl) +1,nsno)
	367	istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn)
	368	dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn)
	369	TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn)
	370	ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn)
	371	eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn)
	372	G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn)
	373	G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn)
	374	agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn)
	375	END DO
	376
	377
	378
	379
	380	C +--Constrains Splitting of too thick Layers
	381	C + =================================================
	382
	383
	384	C +--Search the thickest non-zero Layer
	385	C + ----------------------------------
	386
	387	DO ikl=1,knonv
	388	dzthin(ikl) = 0. ! Arbitrary unrealistic
	389	END DO
	390	DO ikl=1,knonv
	391	DO isn=1,isnoSV(ikl)
	392	isno_n = isnoSV(ikl)-isn+1 ! Snow Normal.Profile
	393	iice_n = iiceSV(ikl)-isn ! Ice Normal.Profile
	394	iiceOK = min(1,max(0,iice_n +1)) ! Ice Switch
	395	dz_dif =( dzsnSV(ikl,isn) ! Actual Profile
	396	. - dz_max ((1-iiceOK)isno_n*isno_n ! Theoretical Profile
	397	. + iiceOK 2. *iice_n) ) !
	398	. /max(dzsnSV(ikl,isn),epsi) !
	399	OKthin = max(zero, !
	400	. sign(unun, !
	401	. dz_dif-dzthin(ikl))) ! 1.=>New thickest Lay.
	402	. * max(0, ! 1 =>.le. isnoSV
	403	. min(1, !
	404	. isnoSV(ikl)-isn +1 )) !
	405	i_thin(ikl) = (1. - OKthin) * i_thin(ikl) ! Update thickest Lay.
	406	. + OKthin * isn ! Index
	407	dzthin(ikl) = (1. - OKthin) * dzthin(ikl) !
	408	. + OKthin * dz_dif !
	409	END DO
	410
	411	isn=max(1,isnoSV(ikl)-3)
	412	if(dzsnSV(ikl,isn)>0.30) then ! surface layer > 30cm
	413	i_thin(ikl) = isn ! XF 04/07/2019
	414	dzthin(ikl) = dzsnSV(ikl,isn)
	415	endif
	416
	417	isn=max(1,isnoSV(ikl)-2)
	418	if(dzsnSV(ikl,isn)>0.10) then ! surface layer > 10cm
	419	i_thin(ikl) = isn ! XF 04/07/2019
	420	dzthin(ikl) = dzsnSV(ikl,isn)
	421	endif
	422
	423	isn=max(1,isnoSV(ikl)-1)
	424	if(dzsnSV(ikl,isn)>0.05) then ! surface layer > 5cm
	425	i_thin(ikl) = isn ! XF 04/07/2019
	426	dzthin(ikl) = dzsnSV(ikl,isn)
	427	endif
	428
	429	isn=max(1,isnoSV(ikl))
	430	if(dzsnSV(ikl,isn)>0.02) then ! surface layer > 2cm
	431	i_thin(ikl) = isn ! XF 04/07/2019
	432	dzthin(ikl) = dzsnSV(ikl,isn)
	433	endif
	434
	435	END DO
	436
	437	DO ikl=1,knonv
	438	ThickL = max(zero, ! 1. => a too thick
	439	. sign(unun,dzthin(ikl) ! Layer exists
	440	. -epsi )) !
	441	. * max(0,1-max(0 , isnoSV(ikl) ! No spliting allowed
	442	. -nsno+1 )) ! if isno > nsno - 1
	443	Agrege(ikl) = ThickL ! 1. => effective split
	444	. * max(0,1-max(0 , NLaysv(ikl) !
	445	. +isnoSV(ikl) !
	446	. -nsno+1 )) !
	447	NLay_s(ikl) = ThickL ! Agregation
	448	. * max(0,1-max(0 , NLaysv(ikl) ! to allow Splitting
	449	. +isnoSV(ikl) ! at next Time Step
	450	. -nsno )) !
	451	. -Agrege(ikl) !
	452	NLay_s(ikl) = max(0 , NLay_s(ikl)) ! Agregation effective
	453	END DO
	454
	455
	456	C +--Rearranges the Layers
	457	C + ---------------------
	458
	459	DO isn=nsno,2,-1
	460	DO ikl=1,knonv
	461	IF (Agrege(ikl).gt.0..AND.i_thin(ikl).lt.isnoSV(ikl)) THEN
	462	staggr = min(1,max(0,isn-i_thin(ikl) -1))
	463	. * min(1,max(0, isnoSV(ikl)-isn+2))
	464	istoSV(ikl,isn) = staggr * istoSV(ikl ,isn-1)
	465	. + (1. - staggr) * istoSV(ikl ,isn )
	466	dzsnSV(ikl,isn) = staggr * dzsnSV(ikl ,isn-1)
	467	. + (1. - staggr) * dzsnSV(ikl ,isn )
	468	TsisSV(ikl,isn) = staggr * TsisSV(ikl ,isn-1)
	469	. + (1. - staggr) * TsisSV(ikl ,isn )
	470	ro__SV(ikl,isn) = staggr * ro__SV(ikl ,isn-1)
	471	. + (1. - staggr) * ro__SV(ikl ,isn )
	472	eta_SV(ikl,isn) = staggr * eta_SV(ikl ,isn-1)
	473	. + (1. - staggr) * eta_SV(ikl ,isn )
	474	G1snSV(ikl,isn) = staggr * G1snSV(ikl ,isn-1)
	475	. + (1. - staggr) * G1snSV(ikl ,isn )
	476	G2snSV(ikl,isn) = staggr * G2snSV(ikl ,isn-1)
	477	. + (1. - staggr) * G2snSV(ikl ,isn )
	478	agsnSV(ikl,isn) = staggr * agsnSV(ikl ,isn-1)
	479	. + (1. - staggr) * agsnSV(ikl ,isn )
	480	END IF
	481	END DO
	482	END DO
	483
	484	DO ikl=1,knonv
	485	isn = i_thin(ikl)
	486	dzsnSV(ikl,isn) = 0.5Agrege(ikl) dzsnSV(ikl,isn)
	487	. + (1.-Agrege(ikl))*dzsnSV(ikl,isn)
	488
	489	isn = min(i_thin(ikl) +1,nsno)
	490	istoSV(ikl,isn) = Agrege(ikl) *istoSV(ikl,isn-1)
	491	. + (1.-Agrege(ikl))*istoSV(ikl,isn)
	492	dzsnSV(ikl,isn) = Agrege(ikl) *dzsnSV(ikl,isn-1)
	493	. + (1.-Agrege(ikl))*dzsnSV(ikl,isn)
	494	TsisSV(ikl,isn) = Agrege(ikl) *TsisSV(ikl,isn-1)
	495	. + (1.-Agrege(ikl))*TsisSV(ikl,isn)
	496	ro__SV(ikl,isn) = Agrege(ikl) *ro__SV(ikl,isn-1)
	497	. + (1.-Agrege(ikl))*ro__SV(ikl,isn)
	498	eta_SV(ikl,isn) = Agrege(ikl) *eta_SV(ikl,isn-1)
	499	. + (1.-Agrege(ikl))*eta_SV(ikl,isn)
	500	G1snSV(ikl,isn) = Agrege(ikl) *G1snSV(ikl,isn-1)
	501	. + (1.-Agrege(ikl))*G1snSV(ikl,isn)
	502	G2snSV(ikl,isn) = Agrege(ikl) *G2snSV(ikl,isn-1)
	503	. + (1.-Agrege(ikl))*G2snSV(ikl,isn)
	504	agsnSV(ikl,isn) = Agrege(ikl) *agsnSV(ikl,isn-1)
	505	. + (1.-Agrege(ikl))*agsnSV(ikl,isn)
	506	isnoSV(ikl) = min(Agrege(ikl) +isnoSV(ikl),real(nsno))
	507	iiceSV(ikl) = iiceSV(ikl)
	508	. + Agrege(ikl) *max(0,sign(1,iiceSV(ikl)
	509	. -isn +icemix))
	510	. *max(0,sign(1,iiceSV(ikl)
	511	. -1 ))
	512	END DO
	513
	514
	515	C +--Constrains Agregation in case of too much Layers
	516	C + =================================================
	517
	518	C +--Search the thinest non-zero Layer
	519	C + -----------------------------------
	520
	521
	522
	523	DO ikl=1,knonv
	524	dzthin(ikl) = 0. ! Arbitrary unrealistic
	525	END DO ! Layer Thickness
	526	DO ikl=1,knonv
	527	DO isn=1,isnoSV(ikl)-3 ! no agregation of 3 first snowlayers
	528	! XF 04/07/2019
	529
	530	isno_n = isnoSV(ikl)-isn+1 ! Snow Normal.Profile
	531	iice_n = iiceSV(ikl)-isn ! Ice Normal.Profile
	532	iiceOK = min(1,max(0,iice_n +1)) ! Ice Switch
	533	! #vz dz_ref(isn) = !
	534	! #vz. dz_min ((1-iiceOK)isno_n*isno_n ! Theoretical Profile
	535	! #vz. + iiceOK * 2**iice_n) !
	536	! #vz. /max(1,isnoSV(ikl)) !
	537	dz_dif = dz_min ! Actual Profile
	538	. - dzsnSV(ikl ,isn) !
	539	. /max(epsi,((1-iiceOK)isno_nisno_n ! Theoretical Profile
	540	. + iiceOK 2. *iice_n)) !
	541	! #vz dzwdif(isn) = dz_dif !
	542	OKthin = max(zero, !
	543	. sign(unun, !
	544	. dz_dif - dzthin(ikl)))! 1.=> New thinest Lay.
	545	. * max(0, ! 1 => .le. isnoSV
	546	. min(1, !
	547	. isnoSV(ikl)-isn +1 )) !
	548	i_thin(ikl) = (1. - OKthin) * i_thin(ikl) ! Update thinest Lay.
	549	. + OKthin * isn ! Index
	550	dzthin(ikl) = (1. - OKthin) * dzthin(ikl) !
	551	. + OKthin * dz_dif !
	552
	553
	554	END DO
	555	END DO
	556
	557
	558
	559
	560
	561	C +--Index of the contiguous Layer to agregate
	562	C + -----------------------------------------
	563
	564	C + ***************
	565	call SISVAT_zCr
	566	C + ***************
	567
	568
	569	C +--Assign the 2 Layers to agregate
	570	C + -------------------------------
	571
	572	DO ikl=1,knonv
	573	isn = i_thin(ikl)
	574	if(LIndsv(ikl)>0) isn=min(isn, nsno-1) !cXF
	575	isagr1(ikl) = istoSV(ikl,isn)
	576	isagr2(ikl) = istoSV(ikl,isn+LIndsv(ikl))
	577	dzagr1(ikl) = dzsnSV(ikl,isn)
	578	dzagr2(ikl) = dzsnSV(ikl,isn+LIndsv(ikl))
	579	T_agr1(ikl) = TsisSV(ikl,isn)
	580	T_agr2(ikl) = TsisSV(ikl,isn+LIndsv(ikl))
	581	roagr1(ikl) = ro__SV(ikl,isn)
	582	roagr2(ikl) = ro__SV(ikl,isn+LIndsv(ikl))
	583	etagr1(ikl) = eta_SV(ikl,isn)
	584	etagr2(ikl) = eta_SV(ikl,isn+LIndsv(ikl))
	585	G1agr1(ikl) = G1snSV(ikl,isn)
	586	G1agr2(ikl) = G1snSV(ikl,isn+LIndsv(ikl))
	587	G2agr1(ikl) = G2snSV(ikl,isn)
	588	G2agr2(ikl) = G2snSV(ikl,isn+LIndsv(ikl))
	589	agagr1(ikl) = agsnSV(ikl,isn)
	590	agagr2(ikl) = agsnSV(ikl,isn+LIndsv(ikl))
	591	LstLay = min(1,max( 0, isnoSV(ikl)-1 ))
	592	Agrege(ikl) = min(1,
	593	. max(0,
	594	. NLaysv(ikl) +isnoSV(ikl)-nsno
	595	. +NLay_s(ikl) )
	596	. *LstLay )
	597
	598	C + minimum uppermost layer thickness to guarantee a correct reproduction of the snow
	599	C + atmosphere coupling
	600	if(dzsnSV(ikl,max(1,isnoSV(ikl)-0))>0.02 .or. ! surface layers> 2-5-10
	601	. dzsnSV(ikl,max(1,isnoSV(ikl)-1))>0.05 .or. ! XF 04/07/2019
	602	. dzsnSV(ikl,max(1,isnoSV(ikl)-2))>0.10 .or.
	603	. dzsnSV(ikl,max(1,isnoSV(ikl)-3))>0.30 )then
	604	Agrege(ikl) = min(1,
	605	. max(0,
	606	. NLaysv(ikl) +isnoSV(ikl)+1-nsno ! nsno-1 layers ma
	607	. +NLay_s(ikl) )
	608	. *LstLay )
	609	endif
	610
	611	isnoSV(ikl) = isnoSV(ikl)
	612	. -(1-LstLay)*max(zero,
	613	. sign(unun, eps_21
	614	. -dzsnSV(ikl,1) ))
	615	isnoSV(ikl) =max( 0, isnoSV(ikl) )
	616
	617	WEagre(ikl) = 0.
	618	END DO
	619
	620	DO isn=1,nsno
	621	DO ikl=1,knonv
	622	WEagre(ikl) = WEagre(ikl) + ro__SV(ikl,isn)*dzsnSV(ikl,isn)
	623	. *min(1,max(0,i_thin(ikl)+1-isn))
	624	ENDDO
	625	ENDDO
	626
	627	C +--Agregates
	628	C + ---------
	629
	630	C + ***************
	631	call SISVAT_zAg
	632	. (isagr1,isagr2,WEagre
	633	. ,dzagr1,dzagr2,T_agr1,T_agr2
	634	. ,roagr1,roagr2,etagr1,etagr2
	635	. ,G1agr1,G1agr2,G2agr1,G2agr2
	636	. ,agagr1,agagr2,Agrege
	637	. )
	638	C + ***************
	639
	640
	641	C +--Rearranges the Layers
	642	C + ---------------------
	643
	644	C +--New (agregated) Snow layer
	645	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	646	DO ikl=1,knonv
	647	isn = i_thin(ikl)
	648	isn = min(isn,isn+LIndsv(ikl))
	649	isnoSV(ikl) = max(0.,isnoSV(ikl) -Agrege(ikl))
	650	iiceSV(ikl) = iiceSV(ikl)
	651	. -max(0,sign(1,iiceSV(ikl) -isn +icemix))
	652	. *Agrege(ikl)
	653	. *max(0,sign(1,iiceSV(ikl) -1 ))
	654	istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn)
	655	. + Agrege(ikl) *isagr1(ikl)
	656	dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn)
	657	. + Agrege(ikl) *dzagr1(ikl)
	658	TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn)
	659	. + Agrege(ikl) *T_agr1(ikl)
	660	ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn)
	661	. + Agrege(ikl) *roagr1(ikl)
	662	eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn)
	663	. + Agrege(ikl) *etagr1(ikl)
	664	G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn)
	665	. + Agrege(ikl) *G1agr1(ikl)
	666	G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn)
	667	. + Agrege(ikl) *G2agr1(ikl)
	668	agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn)
	669	. + Agrege(ikl) *agagr1(ikl)
	670	END DO
	671
	672	C +--Above
	673	C + ^^^^^
	674	DO ikl=1,knonv
	675	isn1(ikl)=max(i_thin(ikl),i_thin(ikl)+LIndsv(ikl))
	676	END DO
	677	DO i= 1,nsno-1
	678	DO ikl=1,knonv
	679	staggr = min(1,max(0,i +1 -isn1(ikl) ))
	680	istoSV(ikl,i) = (1.-staggr )*istoSV(ikl,i )
	681	. + staggr((1.-Agrege(ikl))istoSV(ikl,i )
	682	. + Agrege(ikl) *istoSV(ikl,i+1))
	683	dzsnSV(ikl,i) = (1.-staggr )*dzsnSV(ikl,i )
	684	. + staggr((1.-Agrege(ikl))dzsnSV(ikl,i )
	685	. + Agrege(ikl) *dzsnSV(ikl,i+1))
	686	TsisSV(ikl,i) = (1.-staggr )*TsisSV(ikl,i )
	687	. + staggr((1.-Agrege(ikl))TsisSV(ikl,i )
	688	. + Agrege(ikl) *TsisSV(ikl,i+1))
	689	ro__SV(ikl,i) = (1.-staggr )*ro__SV(ikl,i )
	690	. + staggr((1.-Agrege(ikl))ro__SV(ikl,i )
	691	. + Agrege(ikl) *ro__SV(ikl,i+1))
	692	eta_SV(ikl,i) = (1.-staggr )*eta_SV(ikl,i )
	693	. + staggr((1.-Agrege(ikl))eta_SV(ikl,i )
	694	. + Agrege(ikl) *eta_SV(ikl,i+1))
	695	G1snSV(ikl,i) = (1.-staggr )*G1snSV(ikl,i )
	696	. + staggr((1.-Agrege(ikl))G1snSV(ikl,i )
	697	. + Agrege(ikl) *G1snSV(ikl,i+1))
	698	G2snSV(ikl,i) = (1.-staggr )*G2snSV(ikl,i )
	699	. + staggr((1.-Agrege(ikl))G2snSV(ikl,i )
	700	. + Agrege(ikl) *G2snSV(ikl,i+1))
	701	agsnSV(ikl,i) = (1.-staggr )*agsnSV(ikl,i )
	702	. + staggr((1.-Agrege(ikl))agsnSV(ikl,i )
	703	. + Agrege(ikl) *agsnSV(ikl,i+1))
	704	END DO
	705	END DO
	706
	707	DO ikl=1,knonv
	708	isn = min(isnoSV(ikl) +1,nsno)
	709	istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn)
	710	dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn)
	711	TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn)
	712	ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn)
	713	eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn)
	714	G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn)
	715	G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn)
	716	agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn)
	717	END DO
	718
[3900]	719
	720	C +--Search new Ice/Snow Interface (option II in MAR)
	721	C + ===============================================
	722
	723	IF (ok_zsn_ii) THEN
	724
	725	DO ikl=1,knonv
	726	iiceSV(ikl) = 0
	727	END DO
[3792]	728
[3900]	729	DO ikl=1,knonv
	730	DO isn=1,isnoSV(ikl)
	731	OK_ICE = max(zero,sign(unun,ro__SV(ikl,isn)-ro_ice+20.))
	732	. * max(zero,sign(unun,dzsnSV(ikl,isn)-epsi))
	733	iiceSV(ikl) = (1.-OK_ICE) *iiceSV(ikl)
	734	. + OK_ICE *isn
	735	END DO
	736	END DO
	737
	738	END IF
	739
[3792]	740	return
	741	end

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

Download in other formats:

Original Format