1 | SUBROUTINE SISVAT_TS2 |
---|
2 | ! #ES. (ETSo_0,ETSo_1,ETSo_d) |
---|
3 | |
---|
4 | ! +------------------------------------------------------------------------+ |
---|
5 | ! | MAR SISVAT_TS2 Mon 16-08-2009 MAR | |
---|
6 | ! | SubRoutine SISVAT_TS2 computes the Soil/Snow temperature and fluxes | |
---|
7 | ! | using the same method as in LMDZ for consistency. | |
---|
8 | ! | The corresponding LMDZ routines are soil (soil.F90) and calcul_fluxs | |
---|
9 | ! | (calcul_fluxs_mod.F90). | |
---|
10 | ! +------------------------------------------------------------------------+ |
---|
11 | ! | | |
---|
12 | ! | | |
---|
13 | ! | PARAMETERS: klonv: Total Number of columns = | |
---|
14 | ! | ^^^^^^^^^^ = Total Number of grid boxes of surface type | |
---|
15 | ! | (land ice for now) | |
---|
16 | ! | | |
---|
17 | ! | INPUT: isnoSV = total Nb of Ice/Snow Layers | |
---|
18 | ! | ^^^^^ sol_SV : Downward Solar Radiation [W/m2] | |
---|
19 | ! | IRd_SV : Surface Downward Longwave Radiation [W/m2] | |
---|
20 | ! | VV__SV : SBL Top Wind Speed [m/s] | |
---|
21 | ! | TaT_SV : SBL Top Temperature [K] | |
---|
22 | ! | QaT_SV : SBL Top Specific Humidity [kg/kg] | |
---|
23 | ! | dzsnSV : Snow Layer Thickness [m] | |
---|
24 | ! | dt__SV : Time Step [s] | |
---|
25 | ! | | |
---|
26 | ! | SoSosv : Absorbed Solar Radiation by Surfac.(Normaliz)[-] | |
---|
27 | ! | Eso_sv : Soil+Snow Emissivity [-] | |
---|
28 | ! | ? rah_sv : Aerodynamic Resistance for Heat [s/m] | |
---|
29 | ! | | |
---|
30 | ! | dz1_SV : "inverse" layer thickness (centered) [1/m] | |
---|
31 | ! | dz2_SV : layer thickness (layer above (?)) [m] | |
---|
32 | ! | AcoHSV : coefficient for Enthalpy evolution, from atm. | |
---|
33 | ! | AcoHSV : coefficient for Enthalpy evolution, from atm. | |
---|
34 | ! | AcoQSV : coefficient for Humidity evolution, from atm. | |
---|
35 | ! | BcoQSV : coefficient for Humidity evolution, from atm. | |
---|
36 | ! | ps__SV : surface pressure [Pa] | |
---|
37 | ! | p1l_SV : 1st atmospheric layer pressure [Pa] | |
---|
38 | ! | cdH_SV : drag coeff Energy (?) | |
---|
39 | ! | rsolSV : Radiation balance surface [W/m2] | |
---|
40 | ! | lambSV : Coefficient for soil layer geometry [-] | |
---|
41 | ! | | |
---|
42 | ! | INPUT / TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)| |
---|
43 | ! | OUTPUT: & Snow Temperatures (layers 1,2,...,nsno) [K] | |
---|
44 | ! | ^^^^^^ rsolSV : Radiation balance surface [W/m2] | |
---|
45 | ! | | |
---|
46 | ! | OUTPUT: IRs_SV : Soil IR Radiation [W/m2] | |
---|
47 | ! | ^^^^^^ HSs_sv : Sensible Heat Flux [W/m2] | |
---|
48 | ! | HLs_sv : Latent Heat Flux [W/m2] | |
---|
49 | ! | TsfnSV : new surface temperature [K] | |
---|
50 | ! | Evp_sv : Evaporation [kg/m2] | |
---|
51 | ! | dSdTSV : Sensible Heat Flux temp. derivative [W/m2/K] | |
---|
52 | ! | dLdTSV : Latent Heat Flux temp. derivative [W/m2/K] | |
---|
53 | ! | | |
---|
54 | ! | ? ETSo_0 : Snow/Soil Energy Power, before Forcing [W/m2] | |
---|
55 | ! | ? ETSo_1 : Snow/Soil Energy Power, after Forcing [W/m2] | |
---|
56 | ! | ? ETSo_d : Snow/Soil Energy Power Forcing [W/m2] | |
---|
57 | ! | | |
---|
58 | ! |________________________________________________________________________| |
---|
59 | |
---|
60 | USE VAR_SV |
---|
61 | USE VARdSV |
---|
62 | |
---|
63 | USE VARySV |
---|
64 | USE VARtSV |
---|
65 | USE VARxSV |
---|
66 | USE VARphy |
---|
67 | USE indice_sol_mod |
---|
68 | USE lmdz_comsoil, ONLY: inertie_sol, inertie_sno, inertie_sic, inertie_lic, iflag_sic, iflag_inertie |
---|
69 | USE lmdz_yoethf |
---|
70 | |
---|
71 | USE lmdz_yomcst |
---|
72 | |
---|
73 | IMPLICIT NONE |
---|
74 | INCLUDE "FCTTRE.h" |
---|
75 | |
---|
76 | |
---|
77 | ! +--Global Variables |
---|
78 | ! + ================ |
---|
79 | |
---|
80 | ! INCLUDE "indicesol.h" |
---|
81 | ! include "LMDZphy.inc" |
---|
82 | |
---|
83 | ! +--OUTPUT for Stand Alone NetCDF File |
---|
84 | ! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
85 | ! #NC real*8 SOsoKL(klonv) ! Absorbed Solar Radiation |
---|
86 | ! #NC real*8 IRsoKL(klonv) ! Absorbed IR Radiation |
---|
87 | ! #NC real*8 HSsoKL(klonv) ! Absorbed Sensible Heat Flux |
---|
88 | ! #NC real*8 HLsoKL(klonv) ! Absorbed Latent Heat Flux |
---|
89 | ! #NC real*8 HLs_KL(klonv) ! Evaporation |
---|
90 | ! #NC real*8 HLv_KL(klonv) ! Transpiration |
---|
91 | ! #NC common/DumpNC/SOsoKL,IRsoKL |
---|
92 | ! #NC . ,HSsoKL,HLsoKL |
---|
93 | ! #NC . ,HLs_KL,HLv_KL |
---|
94 | |
---|
95 | ! +--Internal Variables |
---|
96 | ! + ================== |
---|
97 | |
---|
98 | INTEGER :: ig, jk, isl |
---|
99 | REAL :: mu |
---|
100 | REAL :: Tsrf(klonv) ! surface temperature as extrapolated from soil |
---|
101 | REAL :: mug(klonv) !hj coef top layers |
---|
102 | REAL :: ztherm_i(klonv), zdz2(klonv, -nsol:nsno), z1s |
---|
103 | REAL :: pfluxgrd(klonv), pcapcal(klonv), cal(klonv) |
---|
104 | REAL :: beta(klonv), dif_grnd(klonv) |
---|
105 | REAL :: C_coef(klonv, -nsol:nsno), D_coef(klonv, -nsol:nsno) |
---|
106 | |
---|
107 | REAL, DIMENSION(klonv) :: zx_mh, zx_nh, zx_oh |
---|
108 | REAL, DIMENSION(klonv) :: zx_mq, zx_nq, zx_oq |
---|
109 | REAL, DIMENSION(klonv) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef |
---|
110 | REAL, DIMENSION(klonv) :: zx_sl, zx_k1 |
---|
111 | REAL, DIMENSION(klonv) :: d_ts |
---|
112 | REAL :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh |
---|
113 | REAL :: qsat_new, q1_new |
---|
114 | ! REAL, PARAMETER :: t_grnd = 271.35, t_coup = 273.15 |
---|
115 | ! REAL, PARAMETER :: max_eau_sol = 150.0 |
---|
116 | REAL, DIMENSION(klonv) :: IRs__D, dIRsdT |
---|
117 | |
---|
118 | REAL :: t_grnd ! not used |
---|
119 | parameter(t_grnd = 271.35) ! |
---|
120 | REAL :: t_coup ! distinguish evap/sublimation |
---|
121 | parameter(t_coup = 273.15) ! |
---|
122 | REAL :: max_eau_sol |
---|
123 | parameter(max_eau_sol = 150.0) |
---|
124 | |
---|
125 | |
---|
126 | ! WRITE(*,*)'T check' |
---|
127 | |
---|
128 | ! DO ig = 1,knonv |
---|
129 | ! DO jk = 1,isnoSV(ig) !nsno |
---|
130 | ! IF (TsisSV(ig,jk) <= 1.) THEN !hj check |
---|
131 | ! TsisSV(ig,jk) = TsisSV(ig,isnoSV(ig)) |
---|
132 | ! ENDIF |
---|
133 | |
---|
134 | ! IF (TsisSV(ig,jk) <= 1.) THEN !hj check |
---|
135 | ! TsisSV(ig,jk) = 273.15 |
---|
136 | ! ENDIF |
---|
137 | ! END DO |
---|
138 | ! END DO |
---|
139 | |
---|
140 | !!======================================================================= |
---|
141 | !! I. First part: corresponds to soil.F90 in LMDZ |
---|
142 | !!======================================================================= |
---|
143 | |
---|
144 | DO ig = 1, knonv |
---|
145 | DO jk = 1, isnoSV(ig) |
---|
146 | dz2_SV(ig, jk) = dzsnSV(ig, jk) |
---|
147 | !! use arithmetic center between layers to derive dz1 for snow layers for simplicity: |
---|
148 | dz1_SV(ig, jk) = 2. / (dzsnSV(ig, jk) + dzsnSV(ig, jk - 1)) |
---|
149 | ENDDO |
---|
150 | ENDDO |
---|
151 | |
---|
152 | DO ig = 1, knonv |
---|
153 | ztherm_i(ig) = inertie_lic |
---|
154 | IF (isnoSV(ig) > 0) ztherm_i(ig) = inertie_sno |
---|
155 | ENDDO |
---|
156 | |
---|
157 | !!----------------------------------------------------------------------- |
---|
158 | !! 1) |
---|
159 | !! Calculation of Cgrf and Dgrd coefficients using soil temperature from |
---|
160 | !! previous time step. |
---|
161 | !! |
---|
162 | !! These variables are recalculated on the local compressed grid instead |
---|
163 | !! of saved in restart file. |
---|
164 | !!----------------------------------------------------------------------- |
---|
165 | DO ig = 1, knonv |
---|
166 | DO jk = -nsol, nsno |
---|
167 | zdz2(ig, jk) = dz2_SV(ig, jk) / dt__SV !ptimestep |
---|
168 | ENDDO |
---|
169 | ENDDO |
---|
170 | |
---|
171 | DO ig = 1, knonv |
---|
172 | z1s = zdz2(ig, -nsol) + dz1_SV(ig, -nsol + 1) |
---|
173 | C_coef(ig, -nsol + 1) = zdz2(ig, -nsol) * TsisSV(ig, -nsol) / z1s |
---|
174 | D_coef(ig, -nsol + 1) = dz1_SV(ig, -nsol + 1) / z1s |
---|
175 | ENDDO |
---|
176 | |
---|
177 | DO ig = 1, knonv |
---|
178 | DO jk = -nsol + 1, isnoSV(ig) - 1, 1 |
---|
179 | z1s = 1. / (zdz2(ig, jk) + dz1_SV(ig, jk + 1) + dz1_SV(ig, jk) & |
---|
180 | * (1. - D_coef(ig, jk))) |
---|
181 | C_coef(ig, jk + 1) = & |
---|
182 | (TsisSV(ig, jk) * zdz2(ig, jk) & |
---|
183 | + dz1_SV(ig, jk) * C_coef(ig, jk)) * z1s |
---|
184 | D_coef(ig, jk + 1) = dz1_SV(ig, jk + 1) * z1s |
---|
185 | ENDDO |
---|
186 | ENDDO |
---|
187 | |
---|
188 | !!----------------------------------------------------------------------- |
---|
189 | !! 2) |
---|
190 | !! Computation of the soil temperatures using the Cgrd and Dgrd |
---|
191 | !! coefficient computed above |
---|
192 | !! |
---|
193 | !!----------------------------------------------------------------------- |
---|
194 | !! Extrapolate surface Temperature !hj check |
---|
195 | mu = 1. / ((2.**1.5 - 1.) / (2.**(0.5) - 1.) - 1.) |
---|
196 | |
---|
197 | ! IF (knonv>0) THEN |
---|
198 | ! DO ig=1,8 |
---|
199 | ! WRITE(*,*)ig,'sisvat: Tsis ',TsisSV(ig,isnoSV(ig)) |
---|
200 | ! WRITE(*,*)'max-1 ',TsisSV(ig,isnoSV(ig)-1) |
---|
201 | ! WRITE(*,*)'max-2 ',TsisSV(ig,isnoSV(ig)-2) |
---|
202 | ! WRITE(*,*)'0 ',TsisSV(ig,0) |
---|
203 | !! WRITE(*,*)min(max(isnoSV(ig),0),1),max(1-isnoSV(ig),0) |
---|
204 | ! ENDDO |
---|
205 | ! END IF |
---|
206 | |
---|
207 | DO ig = 1, knonv |
---|
208 | IF (isnoSV(ig)>0) THEN |
---|
209 | IF (isnoSV(ig)>1) THEN |
---|
210 | mug(ig) = 1. / (1. + dzsnSV(ig, isnoSV(ig) - 1) / dzsnSV(ig, isnoSV(ig))) !mu |
---|
211 | ELSE |
---|
212 | mug(ig) = 1. / (1. + dzsnSV(ig, isnoSV(ig) - 1) / dz_dSV(0)) !mu |
---|
213 | ENDIF |
---|
214 | ELSE |
---|
215 | mug(ig) = lambSV |
---|
216 | ENDIF |
---|
217 | |
---|
218 | IF (mug(ig) <= 0.05) THEN |
---|
219 | WRITE(*, *)'Attention mu low', mug(ig) |
---|
220 | ENDIF |
---|
221 | IF (mug(ig) >= 0.98) THEN |
---|
222 | WRITE(*, *)'Attention mu high', mug(ig) |
---|
223 | ENDIF |
---|
224 | |
---|
225 | Tsrf(ig) = (1.5 * TsisSV(ig, isnoSV(ig)) - 0.5 * TsisSV(ig, isnoSV(ig) - 1))& |
---|
226 | * min(max(isnoSV(ig), 0), 1) + & |
---|
227 | ((mug(ig) + 1) * TsisSV(ig, 0) - mug(ig) * TsisSV(ig, -1)) & |
---|
228 | * max(1 - isnoSV(ig), 0) |
---|
229 | ENDDO |
---|
230 | |
---|
231 | |
---|
232 | |
---|
233 | !! Surface temperature |
---|
234 | DO ig = 1, knonv |
---|
235 | TsisSV(ig, isnoSV(ig)) = (mug(ig) * C_coef(ig, isnoSV(ig)) + Tsf_SV(ig)) / & |
---|
236 | (mug(ig) * (1. - D_coef(ig, isnoSV(ig))) + 1.) |
---|
237 | ENDDO |
---|
238 | |
---|
239 | !! Other temperatures |
---|
240 | DO ig = 1, knonv |
---|
241 | DO jk = isnoSV(ig), -nsol + 1, -1 |
---|
242 | TsisSV(ig, jk - 1) = C_coef(ig, jk) + D_coef(ig, jk) & |
---|
243 | * TsisSV(ig, jk) |
---|
244 | ENDDO |
---|
245 | ENDDO |
---|
246 | ! WRITE(*,*)ig,'Tsis',TsisSV(ig,0) |
---|
247 | |
---|
248 | ! IF (indice == is_sic) THEN |
---|
249 | ! DO ig = 1,knonv |
---|
250 | ! TsisSV(ig,-nsol) = RTT - 1.8 |
---|
251 | ! END DO |
---|
252 | ! ENDIF |
---|
253 | |
---|
254 | !C !hj new 11 03 2010 |
---|
255 | DO ig = 1, knonv |
---|
256 | isl = isnoSV(ig) |
---|
257 | ! dIRsdT(ig) = Eso_sv(ig)* SteBo * 4. & ! - d(IR)/d(T) |
---|
258 | ! & * Tsf_SV(ig) & !T TsisSV(ig,isl) ! |
---|
259 | ! & * Tsf_SV(ig) & !TsisSV(ig,isl) ! |
---|
260 | ! & * Tsf_SV(ig) !TsisSV(ig,isl) ! |
---|
261 | ! IRs__D(ig) = dIRsdT(ig)* Tsf_SV(ig) * 0.75 !TsisSV(ig,isl) * 0.75 !: |
---|
262 | dIRsdT(ig) = Eso_sv(ig) * StefBo * 4. & ! - d(IR)/d(T) |
---|
263 | * TsisSV(ig, isl) & ! |
---|
264 | * TsisSV(ig, isl) & ! |
---|
265 | * TsisSV(ig, isl) |
---|
266 | IRs__D(ig) = dIRsdT(ig) * TsisSV(ig, isl) * 0.75 !: |
---|
267 | END DO |
---|
268 | !hj |
---|
269 | !!----------------------------------------------------------------------- |
---|
270 | !! 3) |
---|
271 | !! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil |
---|
272 | !! temperature |
---|
273 | !!----------------------------------------------------------------------- |
---|
274 | DO ig = 1, knonv |
---|
275 | z1s = zdz2(ig, -nsol) + dz1_SV(ig, -nsol + 1) |
---|
276 | C_coef(ig, -nsol + 1) = zdz2(ig, -nsol) * TsisSV(ig, -nsol) / z1s |
---|
277 | D_coef(ig, -nsol + 1) = dz1_SV(ig, -nsol + 1) / z1s |
---|
278 | ENDDO |
---|
279 | |
---|
280 | DO ig = 1, knonv |
---|
281 | DO jk = -nsol + 1, isnoSV(ig) - 1, 1 |
---|
282 | z1s = 1. / (zdz2(ig, jk) + dz1_SV(ig, jk + 1) + dz1_SV(ig, jk) & |
---|
283 | * (1. - D_coef(ig, jk))) |
---|
284 | C_coef(ig, jk + 1) = (TsisSV(ig, jk) * zdz2(ig, jk) + & |
---|
285 | dz1_SV(ig, jk) * C_coef(ig, jk)) * z1s |
---|
286 | D_coef(ig, jk + 1) = dz1_SV(ig, jk + 1) * z1s |
---|
287 | ENDDO |
---|
288 | ENDDO |
---|
289 | |
---|
290 | !!----------------------------------------------------------------------- |
---|
291 | !! 4) |
---|
292 | !! Computation of the surface diffusive flux from ground and |
---|
293 | !! calorific capacity of the ground |
---|
294 | !!----------------------------------------------------------------------- |
---|
295 | DO ig = 1, knonv |
---|
296 | !! (pfluxgrd) |
---|
297 | pfluxgrd(ig) = ztherm_i(ig) * dz1_SV(ig, isnoSV(ig)) * & |
---|
298 | (C_coef(ig, isnoSV(ig)) + (D_coef(ig, isnoSV(ig)) - 1.) & |
---|
299 | * TsisSV(ig, isnoSV(ig))) |
---|
300 | !! (pcapcal) |
---|
301 | pcapcal(ig) = ztherm_i(ig) * & |
---|
302 | (dz2_SV(ig, isnoSV(ig)) + dt__SV * (1. - D_coef(ig, isnoSV(ig))) & |
---|
303 | * dz1_SV(ig, isnoSV(ig))) |
---|
304 | z1s = mug(ig) * (1. - D_coef(ig, isnoSV(ig))) + 1. |
---|
305 | pcapcal(ig) = pcapcal(ig) / z1s |
---|
306 | pfluxgrd(ig) = (pfluxgrd(ig) & |
---|
307 | + pcapcal(ig) * (TsisSV(ig, isnoSV(ig)) * z1s & |
---|
308 | - mug(ig) * C_coef(ig, isnoSV(ig)) & |
---|
309 | - Tsf_SV(ig)) / dt__SV) |
---|
310 | ENDDO |
---|
311 | |
---|
312 | cal(1:knonv) = RCPD / pcapcal(1:knonv) |
---|
313 | rsolSV(1:knonv) = rsolSV(1:knonv) + pfluxgrd(1:knonv) |
---|
314 | !!======================================================================= |
---|
315 | !! II. Second part: corresponds to calcul_fluxs_mod.F90 in LMDZ |
---|
316 | !!======================================================================= |
---|
317 | |
---|
318 | Evp_sv = 0. |
---|
319 | ! #NC HSsoKL=0. |
---|
320 | ! #NC HLsoKL=0. |
---|
321 | dSdTSV = 0. |
---|
322 | dLdTSV = 0. |
---|
323 | |
---|
324 | beta(:) = 1.0 |
---|
325 | dif_grnd(:) = 0.0 |
---|
326 | |
---|
327 | !! zx_qs = qsat en kg/kg |
---|
328 | !!**********************************************************************x*************** |
---|
329 | |
---|
330 | DO ig = 1, knonv |
---|
331 | IF (ps__SV(ig)<1.) THEN |
---|
332 | ! WRITE(*,*)'ig',ig,'ps',ps__SV(ig) |
---|
333 | ps__SV(ig) = max(ps__SV(ig), 1.e-8) |
---|
334 | ENDIF |
---|
335 | IF (p1l_SV(ig)<1.) THEN |
---|
336 | ! WRITE(*,*)'ig',ig,'p1l',p1l_SV(ig) |
---|
337 | p1l_SV(ig) = max(p1l_SV(ig), 1.e-8) |
---|
338 | ENDIF |
---|
339 | IF (TaT_SV(ig)<180.) THEN |
---|
340 | ! WRITE(*,*)'ig',ig,'TaT',TaT_SV(ig) |
---|
341 | TaT_SV(ig) = max(TaT_SV(ig), 180.) |
---|
342 | ENDIF |
---|
343 | IF (QaT_SV(ig)<1.e-8) THEN |
---|
344 | ! WRITE(*,*)'ig',ig,'QaT',QaT_SV(ig) |
---|
345 | QaT_SV(ig) = max(QaT_SV(ig), 1.e-8) |
---|
346 | ENDIF |
---|
347 | IF (Tsf_SV(ig)<100.) THEN |
---|
348 | ! WRITE(*,*)'ig',ig,'Tsf',Tsf_SV(ig) |
---|
349 | Tsf_SV(ig) = max(Tsf_SV(ig), 180.) |
---|
350 | ENDIF |
---|
351 | IF (Tsf_SV(ig)>500.) THEN |
---|
352 | ! WRITE(*,*)'ig',ig,'Tsf',Tsf_SV(ig) |
---|
353 | Tsf_SV(ig) = min(Tsf_SV(ig), 400.) |
---|
354 | ENDIF |
---|
355 | ! IF (Tsrf(ig).LT.1.) THEN |
---|
356 | !! WRITE(*,*)'ig',ig,'Tsrf',Tsrf(ig) |
---|
357 | ! Tsrf(ig)=max(Tsrf(ig),TaT_SV(ig)-20.) |
---|
358 | ! ENDIF |
---|
359 | IF (cdH_SV(ig)<1.e-10) THEN |
---|
360 | ! IF (ig.le.3) WRITE(*,*)'ig',ig,'cdH',cdH_SV(ig) |
---|
361 | cdH_SV(ig) = .5 |
---|
362 | ENDIF |
---|
363 | ENDDO |
---|
364 | |
---|
365 | DO ig = 1, knonv |
---|
366 | zx_pkh(ig) = 1. ! (ps__SV(ig)/ps__SV(ig))**RKAPPA |
---|
367 | IF (thermcep) THEN |
---|
368 | zdelta = MAX(0., SIGN(1., rtt - Tsf_SV(ig))) |
---|
369 | zcvm5 = R5LES * LhvH2O * (1. - zdelta) + R5IES * LhsH2O * zdelta |
---|
370 | zcvm5 = zcvm5 / RCPD / (1.0 + RVTMP2 * QaT_SV(ig)) |
---|
371 | zx_qs = r2es * FOEEW(Tsf_SV(ig), zdelta) / ps__SV(ig) |
---|
372 | zx_qs = MIN(0.5, zx_qs) |
---|
373 | !WRITE(*,*)'zcor',retv*zx_qs |
---|
374 | zcor = 1. / (1. - retv * zx_qs) |
---|
375 | zx_qs = zx_qs * zcor |
---|
376 | zx_dq_s_dh = FOEDE(Tsf_SV(ig), zdelta, zcvm5, zx_qs, zcor) & |
---|
377 | / LhvH2O / zx_pkh(ig) |
---|
378 | ELSE |
---|
379 | IF (Tsf_SV(ig)<t_coup) THEN |
---|
380 | zx_qs = qsats(Tsf_SV(ig)) / ps__SV(ig) |
---|
381 | zx_dq_s_dh = dqsats(Tsf_SV(ig), zx_qs) / LhvH2O & |
---|
382 | / zx_pkh(ig) |
---|
383 | ELSE |
---|
384 | zx_qs = qsatl(Tsf_SV(ig)) / ps__SV(ig) |
---|
385 | zx_dq_s_dh = dqsatl(Tsf_SV(ig), zx_qs) / LhvH2O & |
---|
386 | / zx_pkh(ig) |
---|
387 | ENDIF |
---|
388 | ENDIF |
---|
389 | zx_dq_s_dt(ig) = RCPD * zx_pkh(ig) * zx_dq_s_dh |
---|
390 | zx_qsat(ig) = zx_qs |
---|
391 | ! zx_coef(ig) = cdH_SV(ig) * & |
---|
392 | ! & (1.0+SQRT(u1lay(ig)**2+v1lay(ig)**2)) * & |
---|
393 | ! & p1l_SV(ig)/(RD*t1lay(ig)) |
---|
394 | zx_coef(ig) = cdH_SV(ig) * & |
---|
395 | (1.0 + VV__SV(ig)) * & |
---|
396 | p1l_SV(ig) / (RD * TaT_SV(ig)) |
---|
397 | |
---|
398 | ENDDO |
---|
399 | |
---|
400 | |
---|
401 | !! === Calcul de la temperature de surface === |
---|
402 | !! zx_sl = chaleur latente d'evaporation ou de sublimation |
---|
403 | !!**************************************************************************** |
---|
404 | |
---|
405 | DO ig = 1, knonv |
---|
406 | zx_sl(ig) = LhvH2O |
---|
407 | IF (Tsf_SV(ig) < RTT) zx_sl(ig) = LhsH2O |
---|
408 | zx_k1(ig) = zx_coef(ig) |
---|
409 | ENDDO |
---|
410 | |
---|
411 | DO ig = 1, knonv |
---|
412 | !! Q |
---|
413 | zx_oq(ig) = 1. - (beta(ig) * zx_k1(ig) * BcoQSV(ig) * dt__SV) |
---|
414 | zx_mq(ig) = beta(ig) * zx_k1(ig) * & |
---|
415 | (AcoQSV(ig) - zx_qsat(ig) + & |
---|
416 | zx_dq_s_dt(ig) * Tsf_SV(ig)) & |
---|
417 | / zx_oq(ig) |
---|
418 | zx_nq(ig) = beta(ig) * zx_k1(ig) * (-1. * zx_dq_s_dt(ig)) & |
---|
419 | / zx_oq(ig) |
---|
420 | |
---|
421 | !! H |
---|
422 | zx_oh(ig) = 1. - (zx_k1(ig) * BcoHSV(ig) * dt__SV) |
---|
423 | zx_mh(ig) = zx_k1(ig) * AcoHSV(ig) / zx_oh(ig) |
---|
424 | zx_nh(ig) = - (zx_k1(ig) * RCPD * zx_pkh(ig)) / zx_oh(ig) |
---|
425 | |
---|
426 | !! surface temperature |
---|
427 | TsfnSV(ig) = (Tsf_SV(ig) + cal(ig) / RCPD * zx_pkh(ig) * dt__SV * & |
---|
428 | (rsolSV(ig) + zx_mh(ig) + zx_sl(ig) * zx_mq(ig)) & |
---|
429 | + dif_grnd(ig) * t_grnd * dt__SV) / & |
---|
430 | (1. - dt__SV * cal(ig) / (RCPD * zx_pkh(ig)) * & |
---|
431 | (zx_nh(ig) + zx_sl(ig) * zx_nq(ig)) & |
---|
432 | + dt__SV * dif_grnd(ig)) |
---|
433 | |
---|
434 | !hj rajoute 22 11 2010 tuning... |
---|
435 | TsfnSV(ig) = min(RTT + 0.02, TsfnSV(ig)) |
---|
436 | |
---|
437 | d_ts(ig) = TsfnSV(ig) - Tsf_SV(ig) |
---|
438 | |
---|
439 | |
---|
440 | !!== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas |
---|
441 | !!== flux_t est le flux de cpt (energie sensible): j/(m**2 s) |
---|
442 | Evp_sv(ig) = - zx_mq(ig) - zx_nq(ig) * TsfnSV(ig) |
---|
443 | HLs_sv(ig) = - Evp_sv(ig) * zx_sl(ig) |
---|
444 | HSs_sv(ig) = zx_mh(ig) + zx_nh(ig) * TsfnSV(ig) |
---|
445 | |
---|
446 | !! Derives des flux dF/dTs (W m-2 K-1): |
---|
447 | dSdTSV(ig) = zx_nh(ig) |
---|
448 | dLdTSV(ig) = zx_sl(ig) * zx_nq(ig) |
---|
449 | |
---|
450 | |
---|
451 | !hj new 11 03 2010 |
---|
452 | isl = isnoSV(ig) |
---|
453 | ! TsisSV(ig,isl) = TsfnSV(ig) |
---|
454 | IRs_SV(ig) = IRs__D(ig) & ! |
---|
455 | - dIRsdT(ig) * TsfnSV(ig) !TsisSV(ig,isl)? ! |
---|
456 | |
---|
457 | ! hj |
---|
458 | ! #NC SOsoKL(ig) = sol_SV(ig) * SoSosv(ig) ! Absorbed Sol. |
---|
459 | ! #NC IRsoKL(ig) = IRs_SV(ig) & !Up Surf. IR |
---|
460 | ! #NC& + tau_sv(ig) *IRd_SV(ig)*Eso_sv(ig) & !Down Atm IR |
---|
461 | ! #NC& -(1.0-tau_sv(ig)) *0.5*IRv_sv(ig) ! Down Veg IR |
---|
462 | ! #NC HLsoKL(ig) = HLs_sv(ig) |
---|
463 | ! #NC HSsoKL(ig) = HSs_sv(ig) |
---|
464 | ! #NC HLs_KL(ig) = Evp_sv(ig) |
---|
465 | |
---|
466 | !! Nouvelle valeure de l'humidite au dessus du sol |
---|
467 | qsat_new = zx_qsat(ig) + zx_dq_s_dt(ig) * d_ts(ig) |
---|
468 | q1_new = AcoQSV(ig) - BcoQSV(ig) * Evp_sv(ig) * dt__SV |
---|
469 | QaT_SV(ig) = q1_new * (1. - beta(ig)) + beta(ig) * qsat_new |
---|
470 | |
---|
471 | ENDDO |
---|
472 | |
---|
473 | END SUBROUTINE sisvat_ts2 |
---|