Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

sisvat_tso.f90 @ 5447

Last change on this file since 5447 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: 26.3 KB

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

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

Download in other formats:

Original Format