Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

sisvat_zsn.f90 @ 5428

Last change on this file since 5428 was 5246, checked in by abarral, 2 months ago
Convert fixed-form to free-form sources .F -> .{f,F}90 (WIP: some .F remain, will be handled in subsequent commits)
File size: 28.9 KB

Rev	Line
[3792]	1
	2
[5246]	3	subroutine SISVAT_zSn
[3792]	4
[5246]	5	! +------------------------------------------------------------------------+
	6	! \| MAR SISVAT_zSn 12-07-2019 MAR \|
	7	! \| SubRoutine SISVAT_zSn manages the Snow Pack vertical Discretization \|
	8	! \| \|
	9	! +------------------------------------------------------------------------+
	10	! \| \|
	11	! \| PARAMETERS: knonv: Total Number of columns = \|
	12	! \| ^^^^^^^^^^ = Total Number of continental grid boxes \|
	13	! \| X Number of Mosaic Cell per grid box \|
	14	! \| \|
	15	! \| INPUT / NLaysv = New Snow Layer Switch \|
	16	! \| OUTPUT: isnoSV = total Nb of Ice/Snow Layers \|
	17	! \| ^^^^^^ ispiSV = 0,...,nsno: Uppermost Superimposed Ice Layer \|
	18	! \| iiceSV = total Nb of Ice Layers \|
	19	! \| istoSV = 0,...,5 : Snow History (see istdSV data) \|
	20	! \| \|
	21	! \| INPUT / TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)\|
	22	! \| OUTPUT: & Snow Temperatures (layers 1,2,...,nsno) [K] \|
	23	! \| ^^^^^^ ro__SV : Soil/Snow Volumic Mass [kg/m3] \|
	24	! \| eta_SV : Soil/Snow Water Content [m3/m3] \|
	25	! \| dzsnSV : Snow Layer Thickness [m] \|
	26	! \| G1snSV : Dendricity (<0) or Sphericity (>0) of Snow Layer \|
	27	! \| G2snSV : Sphericity (>0) or Size of Snow Layer \|
	28	! \| agsnSV : Snow Age [day] \|
	29	! \| \|
	30	! \| METHOD: 1) Agregate the thinest Snow Layer \|
	31	! \| ^^^^^^ if a new Snow Layer has been precipitated (NLaysv = 1) \|
	32	! \| 2) Divide a too thick Snow Layer except \|
	33	! \| if the maximum Number of Layer is reached \|
	34	! \| in this case forces NLay_s = 1 \|
	35	! \| 3) Agregate the thinest Snow Layer \|
	36	! \| in order to divide a too thick Snow Layer \|
	37	! \| at next Time Step when NLay_s = 1 \|
	38	! \| \|
	39	! +------------------------------------------------------------------------+
[3792]	40
	41
	42
	43
[5246]	44	! +--Global Variables
	45	! + ================
[3792]	46
	47
[5246]	48	use VARphy
	49	use VAR_SV
	50	use VARdSV
	51	use VAR0SV
	52	use VARxSV
	53	use VARySV
	54	use surface_data, only: ok_zsn_ii
[3792]	55
[5246]	56	IMPLICIT NONE
[3792]	57
	58
[5246]	59	! +--Internal Variables
	60	! + ==================
[3792]	61
[5246]	62	integer :: ikl ,isn ,i !
[3792]	63
	64
	65
	66
	67
	68
[5246]	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	! + ! 1 Ice
	78	integer :: isn1 (knonv) ! 1st layer to stagger
	79	real :: staggr ! stagger Switch
[3792]	80
[5246]	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
[3792]	87
[5246]	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
[3792]	107
[3900]	108
[5246]	109	! +--DATA
	110	! + ====
[3900]	111
[5246]	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	! #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	! + CAUTION: dz_max > dz_min*2 is required ! Otherwise re-agregation is
	119	! + ! activated after splitting
[3900]	120
[5246]	121
	122
	123
	124
	125	! +--Constrains Agregation of too thin Layers
	126	! + =================================================
	127
	128	! +--Search the thinest non-zero Layer
	129	! + ----------------------------------
	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	!XF
	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	! + !
	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	! + !
	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	! + ***************
	233	call SISVAT_zCr
	234	! + ***************
	235
	236
	237	! +--Assign the 2 Layers to agregate
	238	! + -------------------------------
	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	! +--Agregates
	286	! + ---------
	287
	288	! + ***************
	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	! + ***************
	297
	298
	299	! +--Rearranges the Layers
	300	! + ---------------------
	301
	302	! +--New (agregated) Snow layer
	303	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	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	! +--Above
	331	! + ^^^^^
	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	! +--Constrains Splitting of too thick Layers
	381	! + =================================================
	382
	383
	384	! +--Search the thickest non-zero Layer
	385	! + ----------------------------------
	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	! +--Rearranges the Layers
	457	! + ---------------------
	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	! +--Constrains Agregation in case of too much Layers
	516	! + =================================================
	517
	518	! +--Search the thinest non-zero Layer
	519	! + -----------------------------------
	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	! +--Index of the contiguous Layer to agregate
	562	! + -----------------------------------------
	563
	564	! + ***************
	565	call SISVAT_zCr
	566	! + ***************
	567
	568
	569	! +--Assign the 2 Layers to agregate
	570	! + -------------------------------
	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	! + minimum uppermost layer thickness to guarantee a correct reproduction of the snow
	599	! + 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	! +--Agregates
	628	! + ---------
	629
	630	! + ***************
	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	! + ***************
	639
	640
	641	! +--Rearranges the Layers
	642	! + ---------------------
	643
	644	! +--New (agregated) Snow layer
	645	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	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	! +--Above
	673	! + ^^^^^
	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
	719
	720	! +--Search new Ice/Snow Interface (option II in MAR)
	721	! + ===============================================
	722
	723	IF (ok_zsn_ii) THEN
	724
	725	DO ikl=1,knonv
	726	iiceSV(ikl) = 0
	727	END DO
	728
	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
	740	return
	741	end subroutine sisvat_zsn

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

Download in other formats:

Original Format