Context Navigation

sisvat_tso.F @ 5453

Last change on this file since 5453 was 4013, checked in by lguez, 3 years ago
Sync latest trunk changes to Ocean_skin
File size: 29.3 KB

Line
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) ! 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	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	C +-- Initilialisation of local arrays
156	C + ================================
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	C +--Heat Conduction Coefficient (zero in the Layers over the highest one)
193	C + ===========================
194	C + ---------------- isl eta_SV, rho C (isl)
195	C +
196	C +--Soil ++++++++++++++++ etaMid, mu (isl)
197	C + ----
198	C + ---------------- 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	C + ! DR97 eq.3.31
224	mu_eta = ist__s mu_eta +ist__w vK_dSV ! Water Bodies
225	C + ! 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	C +--Soil/Snow Interface
238	C + -------------------
239
240	C +--Soil Contribution
241	C + ^^^^^^^^^^^^^^^^^
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	C + ! DR97 eq.3.31
251	mu_eta = ist__s mu_eta +ist__w vK_dSV ! Water Bodies
252
253	C +--Snow Contribution
254	C + ^^^^^^^^^^^^^^^^^
255	mu_sno(ikl) = CdidSV !
256	. (ro__SV(ikl,isl) /ro_Wat) * 1.88 !
257	mu_sno(ikl) = max(epsi,mu_sno(ikl)) !
258	C +... mu_sno : Snow Heat Conductivity Coefficient [Wm/K]
259	C + (Yen 1981, CRREL Rep., 81-10)
260
261	C +--Combined Heat Conductivity
262	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
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	C +--Inverted Heat Capacity
269	C + ^^^^^^^^^^^^^^^^^^^^^^
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	C +--Snow
276	C + ----
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	C +--Combined Heat Conductivity
291	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
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	C +--Inverted Heat Capacity
301	C + ^^^^^^^^^^^^^^^^^^^^^^
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	C +--Uppermost Effective Layer: NO conduction
309	C + ----------------------------------------
310
311	DO ikl=1,knonv
312	mu__dz(ikl,isnoSV(ikl)+1) = 0.0
313	END DO
314
315
316	C +--Energy Budget (IN)
317	C + ==================
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	C +--Tridiagonal Elimination: Set Up
334	C + ===============================
335
336	C +--Soil/Snow Interior
337	C + ^^^^^^^^^^^^^^^^^^
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	C +--Soil lowest Layer
356	C + ^^^^^^^^^^^^^^^^^^
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	C +--Snow highest Layer (dummy!)
372	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^
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	C +--Surface: UPwardIR Heat Flux
393	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^
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	C +--Surface: Richardson Number: T Derivative
403	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
404	c #RC dRidTs(ikl) =-gravit * za__SV(ikl)
405	c #RC. /(TaT_SV(ikl) * VV__SV(ikl)
406	c #RC. * VV__SV(ikl))
407
408	C +--Surface: Turbulent Heat Flux: Factors
409	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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	C +--Surface: Sensible Heat Flux: T Derivative
414	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
415	dSdTSV(ikl) = f_HSHL(ikl) * Cp !#- d(HS)/d(T)
416	c #RC. *(1.0 -(TsisSV(ikl,isl) -TaT_SV(ikl)) !#Richardson
417	c #RC. * dRidTs(ikl)*dFh_sv(ikl)/rah_sv(ikl)) ! Nb. Correct.
418	HS___D(ikl) = dSdTSV(ikl) * TaT_SV(ikl) !
419
420	C +--Surface: Latent Heat Flux: Saturation Specific Humidity
421	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
422	c den_qs = TsisSV(ikl,isl)- 35.8 !
423	c arg_qs = 17.27 *(TsisSV(ikl,isl)-273.16) !
424	c . / den_qs !
425	c 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	c 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	C +--Surface: Latent Heat Flux: Surface Relative Humidity
452	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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
459	DO ikl=1,knonv
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) !
466	END DO
467
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	c #VX faceta(ikl) = faceta(ikl) !
475	c #VX. /(1.+faceta(ikl)*dQa_SV(ikl)) ! Limitation
476	! by Atm.Conten
477	c #??. *max(0,sign(1.,qsatsg(ikl)-QaT_SV(ikl)))) ! NO Limitation
478	! of Downw.Flux
479	END DO !
480	DO itEuBk=1,2 !
481	DO ikl=1,knonv
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	c . /(Ro_Wat*dz_dSV(0)) !
503	. * dt_srf /(Ro_Wat*dz_dSV(0)) !
504	cXF 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
512
513
514	endif !
515
516	C +--Surface: Latent Heat Flux: Soil/Water Surface Contributions
517	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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
530
531	c #DL RHuSol(ikl) =(QaT_SV(ikl) !
532	c #DL. -HL___D(ikl) / f___HL(ikl)) !
533	c #DL. / qsatsg(ikl) !
534
535	C +--Surface: Latent Heat Flux: T Derivative
536	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
537	dLdTSV(ikl) = 0.
538	c #DL dLdTSV(ikl) = f___HL(ikl) * RHuSol(ikl) *dqs_dT(ikl) ! - d(HL)/d(T)
539	c #DL HL___D(ikl) = HL___D(ikl) !
540	c #DL. +dLdTSV(ikl) * TsisSV(ikl,isl) !
541	END DO !
542
543	C +--Surface: Tridiagonal Matrix Set Up
544	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
545	DO ikl=1,knonv
546	isl = isnoSV(ikl)
547	TSurf0(ikl) = TsisSV(ikl,isl)
548
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.
558
559	Term_D(ikl,isl) = Explic Elem_A TsisSV(ikl,isl-1)
560	. +(1.0d+0 -Explic Elem_A)TsisSV(ikl,isl)
561
562
563
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.
572
573	END DO
574
575
576	C +--Tridiagonal Elimination
577	C + =======================
578
579	C +--Forward Sweep
580	C + ^^^^^^^^^^^^^^
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
585
586	DO ikl= 1,knonv
587
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
594
595	DO ikl= 1,knonv
596	TsisSV(ikl,-nsol) = Term_D(ikl,-nsol)/Aux__P(ikl,-nsol)
597	END DO
598
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)
604
605
606	END DO
607	END DO
608
609	C +--Backward Sweep
610	C + ^^^^^^^^^^^^^^
611	DO ikl= 1,knonv
612	DO isl=min(nsno-1,isnoSV(ikl)+1),-nsol,-1
613
614
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
618
619	TsisSV(ikl,isl) = min(TaT_SV(ikl)+30,TsisSV(ikl,isl))
620	TsisSV(ikl,isl) = max(TaT_SV(ikl)-30,TsisSV(ikl,isl))
621
622
623	c #EU TsisSV(ikl,isl) = max(TaT_SV(ikl)-15.,TsisSV(ikl,isl))
624
625	!XF 18/11/2018 to avoid ST reaching 70�C!!
626	!It is an error compensation but does not work over tundra
627
628	endif
629
630
631
632	END DO
633
634	END DO
635
636
637
638	C +--Temperature Limits (avoids problems in case of no Snow Layers)
639	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
640	DO ikl= 1,knonv
641	isl = isnoSV(ikl)
642
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	C +--Update Surface Fluxes
660	C + ========================
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

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

Download in other formats:

Original Format