Context Navigation

sisvat_tso.F @ 3817

Last change on this file since 3817 was 3792, checked in by evignon, 4 years ago
Ajout de INLANDSIS, nouvelle interface entre LMDZ et la neige de SISVAT Etienne, 04/01/2021
File size: 27.7 KB

Rev	Line
[3792]	1
	2
	3
	4
	5	subroutine SISVAT_TSo
	6	! #e1. (ETSo_0,ETSo_1,ETSo_d)
	7
	8	C +------------------------------------------------------------------------+
	9	C \| MAR SISVAT_TSo 06-10-2020 MAR \|
	10	C \| SubRoutine SISVAT_TSo computes the Soil/Snow Energy Balance \|
	11	C +------------------------------------------------------------------------+
	12	C \| \|
	13	C \| PARAMETERS: knonv: Total Number of columns = \|
	14	C \| ^^^^^^^^^^ = Total Number of continental grid boxes \|
	15	C \| X Number of Mosaic Cell per grid box \|
	16	C \| \|
	17	C \| INPUT: isotSV = 0,...,11: Soil Type \|
	18	C \| ^^^^^ 0: Water, Solid or Liquid \|
	19	C \| isnoSV = total Nb of Ice/Snow Layers \|
	20	C \| dQa_SV = Limitation of Water Vapor Turbulent Flux \|
	21	C \| \|
	22	C \| INPUT: sol_SV : Downward Solar Radiation [W/m2] \|
	23	C \| ^^^^^ IRd_SV : Surface Downward Longwave Radiation [W/m2] \|
	24	C \| za__SV : SBL Top Height [m] \|
	25	C \| VV__SV : SBL Top Wind Speed [m/s] \|
	26	C \| TaT_SV : SBL Top Temperature [K] \|
	27	C \| rhT_SV : SBL Top Air Density [kg/m3] \|
	28	C \| QaT_SV : SBL Top Specific Humidity [kg/kg] \|
	29	C \| LSdzsv : Vertical Discretization Factor [-] \|
	30	C \| = 1. Soil \|
	31	C \| = 1000. Ocean \|
	32	C \| dzsnSV : Snow Layer Thickness [m] \|
	33	C \| ro__SV : Snow/Soil Volumic Mass [kg/m3] \|
	34	C \| eta_SV : Soil Water Content [m3/m3] \|
	35	C \| dt__SV : Time Step [s] \|
	36	C \| \|
	37	C \| SoSosv : Absorbed Solar Radiation by Surfac.(Normaliz)[-] \|
	38	C \| Eso_sv : Soil+Snow Emissivity [-] \|
	39	C \| rah_sv : Aerodynamic Resistance for Heat [s/m] \|
	40	C \| Lx_H2O : Latent Heat of Vaporization/Sublimation [J/kg] \|
	41	C \| sEX_sv : Verticaly Integrated Extinction Coefficient [-] \|
	42	C \| \|
	43	C \| INPUT / TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)\|
	44	C \| OUTPUT: & Snow Temperatures (layers 1,2,...,nsno) [K] \|
	45	C \| ^^^^^^ \|
	46	C \| \|
	47	C \| OUTPUT: IRs_SV : Soil IR Radiation [W/m2] \|
	48	C \| ^^^^^^ HSs_sv : Sensible Heat Flux [W/m2] \|
	49	C \| HLs_sv : Latent Heat Flux [W/m2] \|
	50	C \| ETSo_0 : Snow/Soil Energy Power, before Forcing [W/m2] \|
	51	C \| ETSo_1 : Snow/Soil Energy Power, after Forcing [W/m2] \|
	52	C \| ETSo_d : Snow/Soil Energy Power Forcing [W/m2] \|
	53	C \| \|
	54	C \| Internal Variables: \|
	55	C \| ^^^^^^^^^^^^^^^^^^ \|
	56	C \| \|
	57	C \| METHOD: NO Skin Surface Temperature \|
	58	C \| ^^^^^^ Semi-Implicit Crank Nicholson Scheme \|
	59	C \| \|
	60	C \| # OPTIONS: #E0: Energy Budget Verification \|
	61	C \| # ^^^^^^^ #kd: KDsvat Option:NO Flux Limitor on HL \|
	62	C \| # #KD: KDsvat Option:Explicit Formulation of HL \|
	63	C \| # #NC: OUTPUT for Stand Alone NetCDF File \|
	64	C \| \|
	65	C +------------------------------------------------------------------------+
	66
	67
	68
	69
	70	C +--Global Variables
	71	C + ================
	72
	73	use VARphy
	74	use VAR_SV
	75	use VARdSV
	76	use VARxSV
	77	use VARySV
	78	use VARtSV
	79	use VAR0SV
	80
	81
	82	IMPLICIT NONE
	83
	84
	85	C +--OUTPUT
	86	C + ------
	87
	88	! #e1 real ETSo_0(knonv) ! Soil/Snow Power, before Forcing
	89	! #e1 real ETSo_1(knonv) ! Soil/Snow Power, after Forcing
	90	! #e1 real ETSo_d(knonv) ! Soil/Snow Power, Forcing
	91
	92
	93	C +--Internal Variables
	94	C + ==================
	95
	96	integer ikl ,isl ,jsl ,ist !
	97	integer ist__s,ist__w ! Soil/Water Body Identifier
	98	integer islsgn ! Soil/Snow Surfac.Identifier
	99	real eps__3 ! Arbitrary Low Number
	100	real etaMid,psiMid ! Layer Interface's Humidity
	101	real mu_eta ! Soil thermal Conductivity
	102	real mu_exp ! arg Soil thermal Conductivity
	103	real mu_min ! Min Soil thermal Conductivity
	104	real mu_max ! Max Soil thermal Conductivity
	105	real mu_sno(knonv),mu_aux ! Snow thermal Conductivity
	106	real mu__dz(knonv,-nsol:nsno+1) ! mu_(eta,sno) / dz
	107	real dtC_sv(knonv,-nsol:nsno) ! dt / C
	108	real IRs__D(knonv) ! UpwardIR Previous Iter.Contr.
	109	real dIRsdT(knonv) ! UpwardIR T Derivat.
	110	real f_HSHL(knonv) ! Factor common to HS and HL
	111	real dRidTs(knonv) ! d(Rib)/d(Ts)
	112	real HS___D(knonv) ! Sensible Heat Flux Atm.Contr.
	113	real f___HL(knonv) !
	114	real HL___D(knonv) ! Latent Heat Flux Atm.Contr.
	115	REAL TSurf0(knonv),dTSurf ! Previous Surface Temperature
	116	real qsatsg(knonv) !,den_qs,arg_qs ! Soil Saturat. Spec. Humidity
	117	real dqs_dT(knonv) ! d(qsatsg)/dTv
	118	real Psi( knonv) ! 1st Soil Layer Water Potential
	119	real RHuSol(knonv) ! Soil Surface Relative Humidity
	120	real etaSol ! Soil Surface Humidity
	121	real d__eta ! Soil Surface Humidity Increm.
	122	real Elem_A,Elem_C ! Diagonal Coefficients
	123	real Diag_A(knonv,-nsol:nsno) ! A Diagonal
	124	real Diag_B(knonv,-nsol:nsno) ! B Diagonal
	125	real Diag_C(knonv,-nsol:nsno) ! C Diagonal
	126	real Term_D(knonv,-nsol:nsno) ! Independant Term
	127	real Aux__P(knonv,-nsol:nsno) ! P Auxiliary Variable
	128	real Aux__Q(knonv,-nsol:nsno) ! Q Auxiliary Variable
	129	real Ts_Min,Ts_Max ! Temperature Limits
	130	! #e1 real Exist0 ! Existing Layer Switch
	131	real psat_wat, psat_ice, sp ! computation of qsat
	132
	133	integer nt_srf,it_srf,itEuBk ! HL: Surface Scheme
	134	parameter(nt_srf=10) !
	135	real agpsrf,xgpsrf,dt_srf,dt_ver !
	136	real etaBAK(knonv) !
	137	real etaNEW(knonv) !
	138	real etEuBk(knonv) !
	139	real fac_dt(knonv),faceta(knonv) !
	140	real PsiArg(knonv),SHuSol(knonv) !
	141
	142
	143
	144	C +--Internal DATA
	145	C + =============
	146
	147	data eps__3 / 1.e-3 / ! Arbitrary Low Number
	148	data mu_exp / -0.4343 / ! Soil Thermal Conductivity
	149	data mu_min / 0.172 / ! Min Soil Thermal Conductivity
	150	data mu_max / 2.000 / ! Max Soil Thermal Conductivity
	151	data Ts_Min / 175. / ! Temperature Minimum
	152	data Ts_Max / 300. / ! Temperature Acceptable Maximum
	153	C + ! including Snow Melt Energy
	154
	155
	156
	157	C +--Heat Conduction Coefficient (zero in the Layers over the highest one)
	158	C + ===========================
	159	C + ---------------- isl eta_SV, rho C (isl)
	160	C +
	161	C +--Soil ++++++++++++++++ etaMid, mu (isl)
	162	C + ----
	163	C + ---------------- isl-1 eta_SV, rho C (isl-1)
	164	isl=-nsol
	165	DO ikl=1,knonv
	166
	167	mu__dz(ikl,isl) = 0.
	168
	169	dtC_sv(ikl,isl) = dtz_SV2(isl) * dt__SV ! dt / (dz X rho C)
	170	. /((rocsSV(isotSV(ikl)) ! [s / (m.J/m3/K)]
	171	. +rcwdSV*eta_SV(ikl,isl)) !
	172	. *LSdzsv(ikl) ) !
	173	END DO
	174	DO isl=-nsol+1,0
	175	DO ikl=1,knonv
	176	ist = isotSV(ikl) ! Soil Type
	177	ist__s = min(ist, 1) ! 1 => Soil
	178	ist__w = 1 - ist__s ! 1 => Water Body
	179
	180	etaMid = 0.5(dz_dSV(isl-1)eta_SV(ikl,isl-1) ! eta at layers
	181	. +dz_dSV(isl) *eta_SV(ikl,isl) ) ! interface
	182	. /dzmiSV(isl) ! LSdzsv implicit !
	183	etaMid = max(etaMid,epsi)
	184	psiMid = psidSV(ist)
	185	. (etadSV(ist)/etaMid)*bCHdSV(ist)
	186	mu_eta = 3.82 (psiMid)*mu_exp ! Soil Thermal
	187	mu_eta = min(max(mu_eta, mu_min), mu_max) ! Conductivity
	188	C + ! DR97 eq.3.31
	189	mu_eta = ist__s mu_eta +ist__w vK_dSV ! Water Bodies
	190	C + ! Correction
	191	mu__dz(ikl,isl) = mu_eta/(dzmiSV(isl) !
	192	. *LSdzsv(ikl)) !
	193
	194	dtC_sv(ikl,isl) = dtz_SV2(isl)* dt__SV ! dt / (dz X rho C)
	195	. /((rocsSV(isotSV(ikl)) !
	196	. +rcwdSV*eta_SV(ikl,isl)) !
	197	. *LSdzsv(ikl) ) !
	198	END DO
	199	END DO
	200
	201
	202	C +--Soil/Snow Interface
	203	C + -------------------
	204
	205	C +--Soil Contribution
	206	C + ^^^^^^^^^^^^^^^^^
	207	isl=1
	208	DO ikl=1,knonv
	209	ist = isotSV(ikl) ! Soil Type
	210	ist__s = min(ist, 1) ! 1 => Soil
	211	ist__w = 1 - ist__s ! 1 => Water Body
	212	psiMid = psidSV(ist) ! Snow => Saturation
	213	mu_eta = 3.82 (psiMid)*mu_exp ! Soil Thermal
	214	mu_eta = min(max(mu_eta, mu_min), mu_max) ! Conductivity
	215	C + ! DR97 eq.3.31
	216	mu_eta = ist__s mu_eta +ist__w vK_dSV ! Water Bodies
	217
	218	C +--Snow Contribution
	219	C + ^^^^^^^^^^^^^^^^^
	220	mu_sno(ikl) = CdidSV !
	221	. (ro__SV(ikl,isl) /ro_Wat) * 1.88 !
	222	mu_sno(ikl) = max(epsi,mu_sno(ikl)) !
	223	C +... mu_sno : Snow Heat Conductivity Coefficient [Wm/K]
	224	C + (Yen 1981, CRREL Rep., 81-10)
	225
	226	C +--Combined Heat Conductivity
	227	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	228	mu__dz(ikl,isl) = 2./(dzsnSV(ikl,isl ) ! Combined Heat
	229	. /mu_sno(ikl) ! Conductivity
	230	. +LSdzsv(ikl) !
	231	. *dz_dSV( isl-1)/mu_eta) ! Coefficient
	232
	233	C +--Inverted Heat Capacity
	234	C + ^^^^^^^^^^^^^^^^^^^^^^
	235	dtC_sv(ikl,isl) = dt__SV/max(epsi, ! dt / (dz X rho C)
	236	. dzsnSV(ikl,isl) * ro__SV(ikl,isl) *Cn_dSV) !
	237	END DO
	238
	239
	240	C +--Snow
	241	C + ----
	242
	243	DO ikl=1,knonv
	244	DO isl=1,min(nsno,isnoSV(ikl)+1)
	245	ro__SV(ikl,isl) = !
	246	. ro__SV(ikl ,isl) !
	247	. * max(0,min(isnoSV(ikl)-isl+1,1)) !
	248
	249	END DO
	250	END DO
	251
	252	DO ikl=1,knonv
	253	DO isl=1,min(nsno,isnoSV(ikl)+1)
	254
	255	C +--Combined Heat Conductivity
	256	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	257	mu_aux = CdidSV !
	258	. (ro__SV(ikl,isl) /ro_Wat) * 1.88 !
	259	mu__dz(ikl,isl) = !
	260	. 2. mu_auxmu_sno(ikl) ! Combined Heat
	261	. /max(epsi,dzsnSV(ikl,isl )*mu_sno(ikl) ! Conductivity
	262	. +dzsnSV(ikl,isl-1)*mu_aux ) ! For upper Layer
	263	mu_sno(ikl) = mu_aux !
	264
	265	C +--Inverted Heat Capacity
	266	C + ^^^^^^^^^^^^^^^^^^^^^^
	267	dtC_sv(ikl,isl) = dt__SV/max(eps__3, ! dt / (dz X rho C)
	268	. dzsnSV(ikl,isl) * ro__SV(ikl,isl) *Cn_dSV) !
	269	END DO
	270	END DO
	271
	272
	273	C +--Uppermost Effective Layer: NO conduction
	274	C + ----------------------------------------
	275
	276	DO ikl=1,knonv
	277	mu__dz(ikl,isnoSV(ikl)+1) = 0.0
	278	END DO
	279
	280
	281	C +--Energy Budget (IN)
	282	C + ==================
	283
	284	! #e1 DO ikl=1,knonv
	285	! #e1 ETSo_0(ikl) = 0.
	286	! #e1 END DO
	287	! #e1 DO isl= -nsol,nsno
	288	! #e1 DO ikl=1,knonv
	289	! #e1 Exist0 = isl - isnoSV(ikl)
	290	! #e1 Exist0 = 1. - max(zero,min(unun,Exist0))
	291	! #e1 ETSo_0(ikl) = ETSo_0(ikl)
	292	! #e1. +(TsisSV(ikl,isl)-TfSnow)*Exist0
	293	! #e1. /dtC_sv(ikl,isl)
	294	! #e1 END DO
	295	! #e1 END DO
	296
	297
	298	C +--Tridiagonal Elimination: Set Up
	299	C + ===============================
	300
	301	C +--Soil/Snow Interior
	302	C + ^^^^^^^^^^^^^^^^^^
	303	DO ikl=1,knonv
	304	DO isl=-nsol+1,min(nsno-1,isnoSV(ikl)+1)
	305
	306	Elem_A = dtC_sv(ikl,isl) *mu__dz(ikl,isl)
	307	Elem_C = dtC_sv(ikl,isl) *mu__dz(ikl,isl+1)
	308	Diag_A(ikl,isl) = -Elem_A *Implic
	309	Diag_C(ikl,isl) = -Elem_C *Implic
	310	Diag_B(ikl,isl) = 1.0d+0 -Diag_A(ikl,isl)-Diag_C(ikl,isl)
	311	Term_D(ikl,isl) = Explic (Elem_A TsisSV(ikl,isl-1)
	312	. +Elem_C *TsisSV(ikl,isl+1))
	313	. +(1.0d+0 -Explic (Elem_A+Elem_C))TsisSV(ikl,isl)
	314	. + dtC_sv(ikl,isl) * sol_SV(ikl) *SoSosv(ikl)
	315	. *(sEX_sv(ikl,isl+1)
	316	. -sEX_sv(ikl,isl ))
	317	END DO
	318	END DO
	319
	320	C +--Soil lowest Layer
	321	C + ^^^^^^^^^^^^^^^^^^
	322	isl= -nsol
	323	DO ikl=1,knonv
	324	Elem_A = 0.
	325	Elem_C = dtC_sv(ikl,isl) *mu__dz(ikl,isl+1)
	326	Diag_A(ikl,isl) = 0.
	327	Diag_C(ikl,isl) = -Elem_C *Implic
	328	Diag_B(ikl,isl) = 1.0d+0 -Diag_A(ikl,isl)-Diag_C(ikl,isl)
	329	Term_D(ikl,isl) = Explic * Elem_C *TsisSV(ikl,isl+1)
	330	. +(1.0d+0 -Explic * Elem_C) *TsisSV(ikl,isl)
	331	. + dtC_sv(ikl,isl) * sol_SV(ikl) *SoSosv(ikl)
	332	. *(sEX_sv(ikl,isl+1)
	333	. -sEX_sv(ikl,isl ))
	334	END DO
	335
	336	C +--Snow highest Layer (dummy!)
	337	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^
	338	isl= min(isnoSV(1)+1,nsno)
	339	DO ikl=1,knonv
	340	Elem_A = dtC_sv(ikl,isl) *mu__dz(ikl,isl)
	341	Elem_C = 0.
	342	Diag_A(ikl,isl) = -Elem_A *Implic
	343	Diag_C(ikl,isl) = 0.
	344	Diag_B(ikl,isl) = 1.0d+0 -Diag_A(ikl,isl)
	345	Term_D(ikl,isl) = Explic * Elem_A *TsisSV(ikl,isl-1)
	346	. +(1.0d+0 -Explic * Elem_A) *TsisSV(ikl,isl)
	347	. + dtC_sv(ikl,isl) * (sol_SV(ikl) *SoSosv(ikl)
	348	. *(sEX_sv(ikl,isl+1)
	349	. -sEX_sv(ikl,isl )))
	350	END DO
	351
	352	C +--Surface: UPwardIR Heat Flux
	353	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^
	354	DO ikl=1,knonv
	355	isl = isnoSV(ikl)
	356	dIRsdT(ikl) = Eso_sv(ikl)* StefBo * 4. ! - d(IR)/d(T)
	357	. * TsisSV(ikl,isl) !
	358	. * TsisSV(ikl,isl) !
	359	. * TsisSV(ikl,isl) !
	360	IRs__D(ikl) = dIRsdT(ikl)* TsisSV(ikl,isl) * 0.75 !
	361
	362	C +--Surface: Richardson Number: T Derivative
	363	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	364	c #RC dRidTs(ikl) =-gravit * za__SV(ikl)
	365	c #RC. /(TaT_SV(ikl) * VV__SV(ikl)
	366	c #RC. * VV__SV(ikl))
	367
	368	C +--Surface: Turbulent Heat Flux: Factors
	369	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	370	f_HSHL(ikl) = rhT_SV(ikl) / rah_sv(ikl) ! to HS, HL
	371	f___HL(ikl) = f_HSHL(ikl) * Lx_H2O(ikl)
	372
	373	C +--Surface: Sensible Heat Flux: T Derivative
	374	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	375	dSdTSV(ikl) = f_HSHL(ikl) * Cp !#- d(HS)/d(T)
	376	c #RC. *(1.0 -(TsisSV(ikl,isl) -TaT_SV(ikl)) !#Richardson
	377	c #RC. * dRidTs(ikl)*dFh_sv(ikl)/rah_sv(ikl)) ! Nb. Correct.
	378	HS___D(ikl) = dSdTSV(ikl) * TaT_SV(ikl) !
	379
	380	C +--Surface: Latent Heat Flux: Saturation Specific Humidity
	381	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	382	c den_qs = TsisSV(ikl,isl)- 35.8 !
	383	c arg_qs = 17.27 *(TsisSV(ikl,isl)-273.16) !
	384	c . / den_qs !
	385	c qsatsg(ikl) = .0038 * exp(arg_qs) !
	386
	387	! sp = (pst_SV(ikl) + ptopSV) * 10.
	388
	389	sp=ps__SV(ikl)
	390	psat_ice = 6.1070 * exp(6150. *(1./273.16 -
	391	. 1./TsisSV(ikl,isl)))
	392
	393	psat_wat = 6.1078 * exp (5.138*log(273.16 /TsisSV(ikl,isl)))
	394	. * exp (6827.*(1. /273.16-1./TsisSV(ikl,isl)))
	395
	396	if(TsisSV(ikl,isl)<=273.16) then
	397	qsatsg(ikl) = 0.622 * psat_ice / (sp - 0.378 * psat_ice)
	398	else
	399	qsatsg(ikl) = 0.622 * psat_wat / (sp - 0.378 * psat_wat)
	400	endif
	401
	402	c dqs_dT(ikl) = qsatsg(ikl)* 4099.2 /(den_qs *den_qs)!
	403	fac_dt(ikl) = f_HSHL(ikl)/(ro_Wat * dz_dSV(0)) !
	404	END DO
	405
	406	C +--Surface: Latent Heat Flux: Surface Relative Humidity
	407	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	408	xgpsrf = 1.05 !
	409	agpsrf = dt__SV*( 1.0-xgpsrf ) !
	410	. /( 1.0-xgpsrf**nt_srf) !
	411	dt_srf = agpsrf !
	412	dt_ver = 0. !
	413	DO ikl=1,knonv
	414	isl = isnoSV(ikl) !
	415	etaBAK(ikl) = max(epsi,eta_SV(ikl ,isl)) !
	416	etaNEW(ikl) = etaBAK(ikl) !
	417	etEuBk(ikl) = etaNEW(ikl) !
	418	END DO !
	419	DO it_srf=1,nt_srf !
	420	dt_ver = dt_ver +dt_srf !
	421	DO ikl=1,knonv !
	422	faceta(ikl) = fac_dt(ikl)*dt_srf !
	423	c #VX faceta(ikl) = faceta(ikl) !
	424	c #VX. /(1.+faceta(ikl)*dQa_SV(ikl)) ! Limitation
	425	! by Atm.Conten
	426	c #??. *max(0,sign(1.,qsatsg(ikl)-QaT_SV(ikl)))) ! NO Limitation
	427	! of Downw.Flux
	428	END DO !
	429	DO itEuBk=1,2 !
	430	DO ikl=1,knonv
	431	ist = max(0,isotSV(ikl)-100*isnoSV(ikl)) ! 0 if H2O
	432	!
	433	Psi(ikl) = !
	434	. psidSV(ist) ! DR97, Eqn 3.34
	435	. *(etadSV(ist) !
	436	. /max(etEuBk(ikl),epsi)) !
	437	. **bCHdSV(ist) !
	438	PsiArg(ikl) = 7.2E-5*Psi(ikl) !
	439	RHuSol(ikl) = exp(-min(0.,PsiArg(ikl))) !
	440	SHuSol(ikl) = qsatsg(ikl) *RHuSol(ikl) ! DR97, Eqn 3.15
	441	etEuBk(ikl) = !
	442	. (etaNEW(ikl) + faceta(ikl)*(QaT_SV(ikl) !
	443	. -SHuSol(ikl) !
	444	. *(1. -bCHdSV(ist) !
	445	. *PsiArg(ikl)) )) !
	446	. /(1. + faceta(ikl)* SHuSol(ikl) !
	447	. *bCHdSV(ist) !
	448	. *PsiArg(ikl) !
	449	. /etaNEW(ikl)) !
	450	etEuBk(ikl) = etEuBk(ikl) !
	451	c . /(Ro_Wat*dz_dSV(0)) !
	452	. * dt_srf /(Ro_Wat*dz_dSV(0)) !
	453	cXF 15/05/2017 BUG
	454	END DO !
	455	END DO !
	456	DO ikl=1,knonv !
	457	etaNEW(ikl) = max(etEuBk(ikl),epsi) !
	458	END DO !
	459	dt_srf = dt_srf * xgpsrf !
	460	END DO !
	461
	462	C +--Surface: Latent Heat Flux: Soil/Water Surface Contributions
	463	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	464	DO ikl=1,knonv !
	465	isl = isnoSV(ikl) !
	466	ist = max(0,isotSV(ikl)-100*isnoSV(ikl)) ! 0 if H2O
	467	ist__s= min(1,ist) ! 1 if no H2O
	468	ist__w= 1-ist__s ! 1 if H2O
	469	d__eta = eta_SV(ikl,isl)-etaNEW(ikl) !
	470	! latent heat flux computation
	471	HL___D(ikl)=( ist__s ro_Wat dz_dSV(0) ! Soil Contrib.
	472	. *(etaNEW(ikl) -etaBAK(ikl)) / dt__SV !
	473	. +ist__w *f_HSHL(ikl) ! H2O Contrib.
	474	. *(QaT_SV(ikl) - qsatsg(ikl)) ) !
	475	. * Lx_H2O(ikl) ! common factor
	476
	477	c #DL RHuSol(ikl) =(QaT_SV(ikl) !
	478	c #DL. -HL___D(ikl) / f___HL(ikl)) !
	479	c #DL. / qsatsg(ikl) !
	480
	481	C +--Surface: Latent Heat Flux: T Derivative
	482	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	483	dLdTSV(ikl) = 0.
	484	c #DL dLdTSV(ikl) = f___HL(ikl) * RHuSol(ikl) *dqs_dT(ikl) ! - d(HL)/d(T)
	485	c #DL HL___D(ikl) = HL___D(ikl) !
	486	c #DL. +dLdTSV(ikl) * TsisSV(ikl,isl) !
	487	END DO !
	488
	489	C +--Surface: Tridiagonal Matrix Set Up
	490	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	491	DO ikl=1,knonv
	492	isl = isnoSV(ikl)
	493	TSurf0(ikl) = TsisSV(ikl,isl)
	494
	495	Elem_A = dtC_sv(ikl,isl)*mu__dz(ikl,isl)
	496	Elem_C = 0.
	497	Diag_A(ikl,isl) = -Elem_A *Implic
	498	Diag_C(ikl,isl) = 0.
	499	Diag_B(ikl,isl) = 1.0d+0 -Diag_A(ikl,isl)
	500	Diag_B(ikl,isl) = Diag_B(ikl,isl)
	501	. + dtC_sv(ikl,isl) * (dIRsdT(ikl) ! Upw. Sol IR
	502	. +dSdTSV(ikl) ! HS/Surf.Contr.
	503	. +dLdTSV(ikl)) ! HL/Surf.Contr.
	504
	505	Term_D(ikl,isl) = Explic Elem_A TsisSV(ikl,isl-1)
	506	. +(1.0d+0 -Explic Elem_A)TsisSV(ikl,isl)
	507
	508
	509
	510	Term_D(ikl,isl) = Term_D(ikl,isl)
	511	. + dtC_sv(ikl,isl) * (sol_SV(ikl) *SoSosv(ikl) ! Absorbed
	512	. *(sEX_sv(ikl,isl+1) ! Solar
	513	. -sEX_sv(ikl,isl ))!
	514	. + IRd_SV(ikl)*Eso_sv(ikl) ! Down Atm IR
	515	. +IRs__D(ikl) ! Upw. Sol IR
	516	. +HS___D(ikl) ! HS/Atmo.Contr.
	517	. +HL___D(ikl) )! HL/Atmo.Contr.
	518
	519	END DO
	520
	521
	522	C +--Tridiagonal Elimination
	523	C + =======================
	524
	525	C +--Forward Sweep
	526	C + ^^^^^^^^^^^^^^
	527	DO ikl= 1,knonv
	528	Aux__P(ikl,-nsol) = Diag_B(ikl,-nsol)
	529	Aux__Q(ikl,-nsol) =-Diag_C(ikl,-nsol)/Aux__P(ikl,-nsol)
	530	END DO
	531
	532	DO ikl= 1,knonv
	533
	534	DO isl=-nsol+1,min(nsno,isnoSV(ikl)+1)
	535	Aux__P(ikl,isl) = Diag_A(ikl,isl) *Aux__Q(ikl,isl-1)
	536	. +Diag_B(ikl,isl)
	537	Aux__Q(ikl,isl) =-Diag_C(ikl,isl) /Aux__P(ikl,isl)
	538	END DO
	539	END DO
	540
	541	DO ikl= 1,knonv
	542	TsisSV(ikl,-nsol) = Term_D(ikl,-nsol)/Aux__P(ikl,-nsol)
	543	END DO
	544
	545	DO ikl= 1,knonv
	546	DO isl=-nsol+1,min(nsno,isnoSV(ikl)+1)
	547	TsisSV(ikl,isl) =(Term_D(ikl,isl)
	548	. -Diag_A(ikl,isl) *TsisSV(ikl,isl-1))
	549	. /Aux__P(ikl,isl)
	550
	551
	552	END DO
	553	END DO
	554
	555	C +--Backward Sweep
	556	C + ^^^^^^^^^^^^^^
	557	DO ikl= 1,knonv
	558	DO isl=min(nsno-1,isnoSV(ikl)+1),-nsol,-1
	559
	560
	561	TsisSV(ikl,isl) = Aux__Q(ikl,isl) *TsisSV(ikl,isl+1)
	562	. +TsisSV(ikl,isl)
	563	if(isl==0.and.isnoSV(ikl)==0) then
	564
	565	TsisSV(ikl,isl) = min(TaT_SV(ikl)+30,TsisSV(ikl,isl))
	566	TsisSV(ikl,isl) = max(TaT_SV(ikl)-30,TsisSV(ikl,isl))
	567
	568
	569	c #EU TsisSV(ikl,isl) = max(TaT_SV(ikl)-15.,TsisSV(ikl,isl))
	570
	571	!XF 18/11/2018 to avoid ST reaching 70�C!!
	572	!It is an error compensation but does not work over tundra
	573
	574	endif
	575
	576
	577
	578	END DO
	579
	580	END DO
	581
	582	C +--Temperature Limits (avoids problems in case of no Snow Layers)
	583	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	584	DO ikl= 1,knonv
	585	isl = isnoSV(ikl)
	586	dTSurf = TsisSV(ikl,isl) - TSurf0(ikl)
	587	TsisSV(ikl,isl) = TSurf0(ikl) + sign(1.,dTSurf) ! 180.0 dgC/hr
	588	. * min(abs(dTSurf),5.e-2*dt__SV) ! =0.05 dgC/s
	589
	590
	591
	592	END DO
	593
	594	DO ikl= 1,knonv
	595	DO isl=min(nsno,isnoSV(ikl)+1),1 ,-1
	596	TsisSV(ikl,isl) = max(Ts_Min, TsisSV(ikl,isl))
	597	TsisSV(ikl,isl) = min(Ts_Max, TsisSV(ikl,isl))
	598	END DO
	599
	600	END DO
	601
	602	C +--Update Surface Fluxes
	603	C + ========================
	604
	605	DO ikl= 1,knonv
	606	isl = isnoSV(ikl)
	607	IRs_SV(ikl) = IRs__D(ikl) !
	608	. - dIRsdT(ikl) * TsisSV(ikl,isl) !
	609	HSs_sv(ikl) = HS___D(ikl) ! Sensible Heat
	610	. - dSdTSV(ikl) * TsisSV(ikl,isl) ! Downward > 0
	611	HLs_sv(ikl) = HL___D(ikl) ! Latent Heat
	612	. - dLdTSV(ikl) * TsisSV(ikl,isl) ! Downward > 0
	613	END DO
	614
	615
	616
	617	return
	618	end

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

Download in other formats:

Original Format