1 | !Bruno Ringeval; routines de calcul des densites de flux de CH4 a l interface surface/atmo |
---|
2 | !routines inspirées TRES fortement du modèle de Walter et al. |
---|
3 | !j indique quand c'est possible les routines du modele initial de Walter a laquelle appartiennent les differents "bouts" ci-dessous |
---|
4 | !modifications principales: modif du calcul du substrat de la methanogenèse |
---|
5 | !couplage avec ORCHIDEE |
---|
6 | |
---|
7 | !Cette routine calcule les densites de flux de CH4 pour un wetlands dont la water table depth est constante dans le temps et situee à pwater_wt2 cm EN-DESSOUS de la surface du sol |
---|
8 | !=> simplification par rapport au modele original (WTD variable) ; les parties inutiles ont ete mises en commentaires quand cela a ete possible (ces parties sont differentes dans stomate_*_ter_0.f90 et stomate_*_ter_weti.f90) |
---|
9 | !Dans cette routine, il y a un calcul d oxydation (sur la partie 0; pwater_wt2 cm) |
---|
10 | |
---|
11 | MODULE stomate_wet_ch4_pt_ter_wet2 |
---|
12 | ! modules used: |
---|
13 | USE ioipsl |
---|
14 | !pss+ |
---|
15 | ! USE stomate_cste_wetlands |
---|
16 | ! USE stomate_constants |
---|
17 | USE constantes |
---|
18 | USE constantes_soil |
---|
19 | USE pft_parameters |
---|
20 | !pss- |
---|
21 | |
---|
22 | IMPLICIT NONE |
---|
23 | |
---|
24 | ! private & public routines |
---|
25 | |
---|
26 | PRIVATE |
---|
27 | PUBLIC ch4_wet_flux_density_wet2,ch4_wet_flux_density_clear_wet2 |
---|
28 | |
---|
29 | ! first call |
---|
30 | LOGICAL, SAVE :: firstcall = .TRUE. |
---|
31 | |
---|
32 | CONTAINS |
---|
33 | |
---|
34 | SUBROUTINE ch4_wet_flux_density_clear_wet2 |
---|
35 | firstcall=.TRUE. |
---|
36 | END SUBROUTINE ch4_wet_flux_density_clear_wet2 |
---|
37 | |
---|
38 | |
---|
39 | SUBROUTINE ch4_wet_flux_density_wet2 ( npts,stempdiag,tsurf,tsurf_year,veget_max,veget,& |
---|
40 | & carbon_surf,lai,uo,uold2,ch4_flux_density_tot, ch4_flux_density_dif, & |
---|
41 | & ch4_flux_density_bub, ch4_flux_density_pla, catm) |
---|
42 | |
---|
43 | |
---|
44 | !BR: certaines variables du modele de Walter ne sont pas "consistent" avec |
---|
45 | !celles d ORCHIDEE: rprof (root depth, rprof), |
---|
46 | !pr le moment conserve le calcul fait au sein du modele de Walter |
---|
47 | !le travail de coéhrence entre les 2 modèles a été fait pour le LAI (j ai enlevé |
---|
48 | !le calcul du LAI au sein du modele de Walter) |
---|
49 | |
---|
50 | ! |
---|
51 | ! 0 declarations |
---|
52 | ! |
---|
53 | |
---|
54 | ! 0.1 input |
---|
55 | |
---|
56 | ! Domain size |
---|
57 | INTEGER(i_std), INTENT(in) :: npts |
---|
58 | ! natural space |
---|
59 | !REAL(r_std), DIMENSION(npts), INTENT(in) :: space_nat |
---|
60 | ! Soil temperature (K) |
---|
61 | ! REAL(r_std),DIMENSION (npts,nbdl), INTENT (in) :: tsoil |
---|
62 | |
---|
63 | REAL(r_std),DIMENSION (npts,nslm), INTENT (in) :: stempdiag |
---|
64 | |
---|
65 | ! "maximal" coverage fraction of a PFT (LAI -> infinity) on nat/agri ground |
---|
66 | REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max |
---|
67 | REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget |
---|
68 | ! carbon pool: active, slow, or passive, natural and agricultural (gC/m**2 of |
---|
69 | ! natural or agricultural ground) |
---|
70 | !pss+++++ |
---|
71 | REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(in) :: carbon_surf |
---|
72 | ! REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(in) :: carbon |
---|
73 | !pss----- |
---|
74 | |
---|
75 | ! temperature (K) at the surface |
---|
76 | REAL(r_std), DIMENSION(npts), INTENT(in) :: tsurf |
---|
77 | |
---|
78 | !modif bruno |
---|
79 | ! temperature (K) at the surface |
---|
80 | REAL(r_std), DIMENSION(npts), INTENT(in) :: tsurf_year |
---|
81 | !end modif bruno |
---|
82 | |
---|
83 | ! LAI |
---|
84 | REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lai |
---|
85 | |
---|
86 | |
---|
87 | REAL(r_std),INTENT (in) :: catm |
---|
88 | |
---|
89 | ! 0.2 modified fields |
---|
90 | |
---|
91 | !BR: a faire: rajouter le carbone du sol dans les variables modifiees |
---|
92 | |
---|
93 | REAL(r_std), DIMENSION(npts,nvert), INTENT(inout) :: uo |
---|
94 | REAL(r_std), DIMENSION(npts,nvert), INTENT(inout) :: uold2 |
---|
95 | |
---|
96 | ! 0.3 output |
---|
97 | |
---|
98 | ! density flux of methane calculated for entire pixel (gCH4/dt/m**2) |
---|
99 | REAL(r_std), DIMENSION(npts), INTENT(out) :: ch4_flux_density_tot |
---|
100 | REAL(r_std), DIMENSION(npts), INTENT(out) :: ch4_flux_density_dif |
---|
101 | REAL(r_std), DIMENSION(npts), INTENT(out) :: ch4_flux_density_bub |
---|
102 | REAL(r_std), DIMENSION(npts), INTENT(out) :: ch4_flux_density_pla |
---|
103 | |
---|
104 | ! 0.4 local |
---|
105 | |
---|
106 | ! ajout |
---|
107 | REAL(r_std), DIMENSION(npts) :: substrat |
---|
108 | REAL(r_std) :: mwater |
---|
109 | |
---|
110 | ! Index |
---|
111 | INTEGER(i_std) :: i |
---|
112 | INTEGER(i_std) :: j |
---|
113 | INTEGER(i_std), DIMENSION(npts) :: j_point |
---|
114 | INTEGER(i_std) :: ipoint |
---|
115 | |
---|
116 | !variables appelees dans l ancien Scalc |
---|
117 | !Vcalc.var |
---|
118 | INTEGER(i_std) :: nsoil |
---|
119 | INTEGER(i_std) :: nroot |
---|
120 | ! INTEGER(i_std) :: ihelp |
---|
121 | !Cmodel.cb |
---|
122 | !new dimension |
---|
123 | REAL(r_std), DIMENSION(npts) :: rcmax |
---|
124 | REAL(r_std),DIMENSION(npts, ns) :: forg |
---|
125 | !on definit uo dans modified fields: doit etre conserve d une iteration a l autre |
---|
126 | !Cread_carbon.cb |
---|
127 | !new dimension |
---|
128 | INTEGER(i_std),DIMENSION(npts) :: ibare |
---|
129 | INTEGER(i_std),DIMENSION(npts) :: iroot |
---|
130 | INTEGER(i_std),DIMENSION(npts) :: isoil |
---|
131 | !new dimension |
---|
132 | REAL(r_std),DIMENSION(npts) :: tveg |
---|
133 | !new dimension |
---|
134 | REAL(r_std),DIMENSION(npts) :: tmean |
---|
135 | ! REAL(r_std),DIMENSION(npts) :: std |
---|
136 | ! REAL(r_std),DIMENSION(npts) :: std3 |
---|
137 | ! REAL(r_std),DIMENSION(npts) :: std4 |
---|
138 | ! REAL(r_std),DIMENSION(npts) :: std5 |
---|
139 | !new dimension |
---|
140 | REAL(r_std),DIMENSION(npts) :: quotient |
---|
141 | REAL(r_std) :: ihelp |
---|
142 | !redondance et pas meme dimension |
---|
143 | !new dimension |
---|
144 | REAL(r_std),DIMENSION(npts,nvert) :: root |
---|
145 | !new dimension |
---|
146 | REAL(r_std),DIMENSION(npts,ns) :: sou |
---|
147 | !ancien Soutput |
---|
148 | !new dimension |
---|
149 | REAL(r_std),DIMENSION(npts) :: fdiffday |
---|
150 | REAL(r_std),DIMENSION(npts) :: fluxbub |
---|
151 | REAL(r_std),DIMENSION(npts) :: fluxplant |
---|
152 | REAL(r_std),DIMENSION(npts) :: fluxtot |
---|
153 | !ancien Vmodel.var --> s y rapporter pour voir les variables que l on a supprime |
---|
154 | !new dimension |
---|
155 | INTEGER(i_std),DIMENSION(npts) :: nw |
---|
156 | INTEGER(i_std),DIMENSION(npts) :: n0 |
---|
157 | INTEGER(i_std),DIMENSION(npts) :: n1 |
---|
158 | INTEGER(i_std) :: n2 |
---|
159 | INTEGER(i_std) :: n3 |
---|
160 | !new dimension |
---|
161 | INTEGER(i_std),DIMENSION(npts) :: nwater |
---|
162 | INTEGER(i_std) :: nthaw |
---|
163 | !redondance INTEGER(i_std) :: nroot |
---|
164 | INTEGER(i_std),DIMENSION(npts) :: j1 |
---|
165 | INTEGER(i_std),DIMENSION(npts) :: j2 |
---|
166 | INTEGER(i_std),DIMENSION(npts) :: j3 |
---|
167 | !new dimension |
---|
168 | INTEGER(i_std),DIMENSION(npts) :: ijump |
---|
169 | INTEGER(i_std),DIMENSION(npts) :: ndim |
---|
170 | INTEGER(i_std) :: iday |
---|
171 | !new dimension |
---|
172 | REAL(r_std),DIMENSION(npts) :: rwater |
---|
173 | REAL(r_std) :: diffsw |
---|
174 | REAL(r_std) :: diffsa |
---|
175 | REAL(r_std) :: diffwa |
---|
176 | !new dimension |
---|
177 | REAL(r_std),DIMENSION(npts) :: diffup |
---|
178 | REAL(r_std),DIMENSION(npts) :: diffdo |
---|
179 | REAL(r_std) :: amhalf |
---|
180 | REAL(r_std) :: aphalf |
---|
181 | REAL(r_std) :: source |
---|
182 | REAL(r_std) :: rexp |
---|
183 | REAL(r_std) :: q |
---|
184 | REAL(r_std) :: p |
---|
185 | !new dimension |
---|
186 | REAL(r_std),DIMENSION(npts) :: rvmax |
---|
187 | REAL(r_std) :: aold1 |
---|
188 | REAL(r_std) :: aold2 |
---|
189 | REAL(r_std) :: aold3 |
---|
190 | REAL(r_std) :: aold4 |
---|
191 | REAL(r_std) :: diffgrdo |
---|
192 | REAL(r_std) :: diffgrup |
---|
193 | REAL(r_std) :: aoldm2 |
---|
194 | REAL(r_std) :: aoldm1 |
---|
195 | REAL(r_std) :: aoldn |
---|
196 | !new dimension |
---|
197 | REAL(r_std),DIMENSION(npts) :: negconc |
---|
198 | REAL(r_std) :: cdiff |
---|
199 | !new dimension |
---|
200 | REAL(r_std),DIMENSION(npts) :: tmax |
---|
201 | REAL(r_std) :: cplant |
---|
202 | !new dimension |
---|
203 | REAL(r_std),DIMENSION(npts) :: tplant |
---|
204 | !new dimension |
---|
205 | REAL(r_std),DIMENSION(npts) :: fgrow |
---|
206 | REAL(r_std) :: cbub |
---|
207 | !new dimension |
---|
208 | REAL(r_std),DIMENSION(npts) :: wpro |
---|
209 | REAL(r_std),DIMENSION(npts) :: tout |
---|
210 | REAL(r_std),DIMENSION(npts) :: goxid |
---|
211 | REAL(r_std),DIMENSION(npts) :: dplant |
---|
212 | REAL(r_std) :: fdiffgrad |
---|
213 | REAL(r_std) :: wtdiff |
---|
214 | REAL(r_std) :: uwp0 |
---|
215 | REAL(r_std) :: uwp1 |
---|
216 | REAL(r_std) :: uwp2 |
---|
217 | REAL(r_std) :: xgrow |
---|
218 | !new dimension |
---|
219 | REAL(r_std),DIMENSION(npts,nvert,nvert) :: a |
---|
220 | REAL(r_std),DIMENSION(npts,nvert) :: b |
---|
221 | REAL(r_std),DIMENSION(npts,nvert) :: uout |
---|
222 | !new dimension |
---|
223 | ! REAL(r_std),DIMENSION(npts,nvert) :: uold2 |
---|
224 | REAL(r_std),DIMENSION(npts,nvert) :: tria |
---|
225 | REAL(r_std),DIMENSION(npts,nvert) :: trib |
---|
226 | REAL(r_std),DIMENSION(npts,nvert) :: tric |
---|
227 | REAL(r_std),DIMENSION(npts,nvert) :: trir |
---|
228 | REAL(r_std),DIMENSION(npts,nvert) :: triu |
---|
229 | !new dimension |
---|
230 | REAL(r_std),DIMENSION(npts,nvert) :: fpart |
---|
231 | REAL(r_std),DIMENSION(npts,nday) :: flplant |
---|
232 | !new dimension |
---|
233 | REAL(r_std),DIMENSION(npts,nday) :: flbub |
---|
234 | REAL(r_std),DIMENSION(npts,nday) :: fdifftime |
---|
235 | !new dimension |
---|
236 | REAL(r_std),DIMENSION(npts,nday) :: flbubdiff |
---|
237 | REAL(r_std),DIMENSION(npts) :: fluxbubdiff |
---|
238 | REAL(r_std),DIMENSION(npts) :: fbubdiff |
---|
239 | REAL(r_std),DIMENSION(npts) :: fluxdiff |
---|
240 | ! REAL(r_std) :: pwater |
---|
241 | |
---|
242 | !variable necessaire depuis l incoporation de la subrooutine tridag (cf. fin du programme) au code |
---|
243 | REAL(r_std),DIMENSION(npts,500) :: gam |
---|
244 | REAL(r_std),DIMENSION(npts) :: bet |
---|
245 | ! REAL(r_std),DIMENSION(npts) :: sum_veget |
---|
246 | REAL(r_std),DIMENSION(npts) :: sum_veget_nat |
---|
247 | REAL(r_std),DIMENSION(npts) :: sum_veget_nat_ss_nu |
---|
248 | !alpha est le pourcentage de carbone labile pouvant subir la methanogenese |
---|
249 | REAL(r_std),DIMENSION(npts) :: alpha_PFT |
---|
250 | REAL(r_std),DIMENSION(3) :: repart_PFT !1:boreaux, 2:temperes, 3:tropciaux |
---|
251 | |
---|
252 | |
---|
253 | !!si je veux rajouter le cas ou je peux avoir de l eau au-dessus du sol |
---|
254 | !=> je dois modifier la valeur max d indices des boucles n2,n3=f(ipoint) |
---|
255 | |
---|
256 | |
---|
257 | |
---|
258 | !2 Ancien Sread |
---|
259 | !attribue les variables d'ORCHIDEE a celles de Walter |
---|
260 | sum_veget_nat(:)=0. |
---|
261 | DO ipoint = 1,npts |
---|
262 | DO j=1,nvm |
---|
263 | IF ( natural(j) ) THEN ! for natural PFTs |
---|
264 | sum_veget_nat(ipoint) = sum_veget_nat(ipoint) + veget_max(ipoint,j) |
---|
265 | ENDIF |
---|
266 | END DO |
---|
267 | END DO |
---|
268 | sum_veget_nat_ss_nu(:)=0. |
---|
269 | DO ipoint = 1,npts |
---|
270 | DO j=2,nvm |
---|
271 | IF ( natural(j) ) THEN ! for natural PFTs |
---|
272 | sum_veget_nat_ss_nu(ipoint) = sum_veget_nat_ss_nu(ipoint) + veget_max(ipoint,j) |
---|
273 | ENDIF |
---|
274 | END DO |
---|
275 | END DO |
---|
276 | !l indice de boucle va jusque 11 pour ne prendre en compte que les PFT naturels |
---|
277 | !! a changer en etiquette "naturels" |
---|
278 | |
---|
279 | !arrondi par defaut (pas de round) |
---|
280 | ibare(:) = INT(100*veget_max(:,1)) |
---|
281 | |
---|
282 | !introduit une variation du alpha selon le type de PFT |
---|
283 | !!!!pss adapt pft_to_mtc style, need to know whether PFT is boreal, temperate or |
---|
284 | !!!!tropical PFT |
---|
285 | repart_PFT(:) = zero |
---|
286 | DO ipoint = 1,npts |
---|
287 | DO j = 2,nvm |
---|
288 | IF ( pft_to_mtc(j)==4 .OR. pft_to_mtc(j)==5 .OR. pft_to_mtc(j)==6 ) THEN |
---|
289 | !temperate |
---|
290 | repart_PFT(2) = repart_PFT(2) + veget_max(ipoint,j) |
---|
291 | ELSEIF ( is_tropical(j) ) THEN |
---|
292 | !tropical |
---|
293 | repart_PFT(3) = repart_PFT(3) + veget_max(ipoint,j) |
---|
294 | ELSE |
---|
295 | !boreal |
---|
296 | repart_PFT(1) = repart_PFT(1) + veget_max(ipoint,j) |
---|
297 | ENDIF |
---|
298 | END DO |
---|
299 | IF ((repart_PFT(1) .GT. repart_PFT(2)) .AND. (repart_PFT(1) .GE. repart_PFT(3))) THEN |
---|
300 | alpha_PFT(ipoint)=alpha_CH4(1) |
---|
301 | ELSEIF ((repart_PFT(2) .GT. repart_PFT(1)) .AND. (repart_PFT(2) .GE. repart_PFT(3))) THEN |
---|
302 | alpha_PFT(ipoint)=alpha_CH4(2) |
---|
303 | ELSEIF ((repart_PFT(3) .GE. repart_PFT(1)) .AND. (repart_PFT(3) .GE. repart_PFT(2))) THEN |
---|
304 | alpha_PFT(ipoint)=alpha_CH4(3) |
---|
305 | elseif ((repart_PFT(2) .EQ. repart_PFT(1)) .AND. & |
---|
306 | & ((repart_PFT(1) .GE. repart_PFT(3)) .OR. (repart_PFT(2) .GE. repart_PFT(3)))) THEN |
---|
307 | alpha_PFT(ipoint)=alpha_CH4(1) |
---|
308 | ENDIF |
---|
309 | |
---|
310 | |
---|
311 | END DO |
---|
312 | |
---|
313 | iroot(:) = 0 |
---|
314 | quotient(:) = 0 |
---|
315 | DO i=1,11 |
---|
316 | DO ipoint = 1,npts |
---|
317 | iroot(ipoint) = iroot(ipoint) + veget_max(ipoint,i)*rdepth_v(i)*tveg_v(i) |
---|
318 | quotient(ipoint) = quotient(ipoint) + veget_max(ipoint,i)*tveg_v(i) |
---|
319 | ENDDO |
---|
320 | END DO |
---|
321 | |
---|
322 | DO ipoint = 1,npts |
---|
323 | IF ((veget_max(ipoint,1) .LE. 0.95) .AND. (sum_veget_nat_ss_nu(ipoint) .NE. 0.)) THEN |
---|
324 | iroot(ipoint) = iroot(ipoint)/quotient(ipoint) |
---|
325 | ELSE |
---|
326 | iroot(ipoint) = 0.0 |
---|
327 | ENDIF |
---|
328 | ENDDO |
---|
329 | !non necessaire de diviser numerateur et denominateur par sum_veget_nat |
---|
330 | |
---|
331 | isoil(:)=0 |
---|
332 | DO i=1,11 |
---|
333 | DO ipoint = 1,npts |
---|
334 | isoil(ipoint)=isoil(ipoint)+veget_max(ipoint,i)*sdepth_v(i) |
---|
335 | ENDDO |
---|
336 | ENDDO |
---|
337 | DO ipoint = 1,npts |
---|
338 | IF (sum_veget_nat(ipoint) .NE. 0.) THEN |
---|
339 | isoil(ipoint) = isoil(ipoint)/sum_veget_nat(ipoint) |
---|
340 | ENDIF |
---|
341 | ENDDO |
---|
342 | |
---|
343 | |
---|
344 | tveg(:)=0 |
---|
345 | DO i=1,11 |
---|
346 | DO ipoint = 1,npts |
---|
347 | tveg(ipoint) = tveg(ipoint)+veget_max(ipoint,i)*tveg_v(i) |
---|
348 | ENDDO |
---|
349 | ENDDO |
---|
350 | DO ipoint = 1,npts |
---|
351 | IF (sum_veget_nat(ipoint) .NE. 0.) THEN |
---|
352 | tveg(ipoint) = tveg(ipoint)/sum_veget_nat(ipoint) |
---|
353 | ENDIF |
---|
354 | ENDDO |
---|
355 | tveg(:) = 0.5*tveg(:)/15. |
---|
356 | |
---|
357 | !ATTENTION: TRES IMPORTANT |
---|
358 | !il n est pas clair dans la publi de Walter et al. si Tmean varie au cours du temps ou non |
---|
359 | !=>cf. ma these |
---|
360 | !cf. stomate_season.f90 |
---|
361 | DO ipoint = 1,npts |
---|
362 | tmean(ipoint) = tsurf_year(ipoint)-273.16 |
---|
363 | ENDDO |
---|
364 | |
---|
365 | !std,std3,std4,std5: supprimees lors du couplage avec ORCHIDEE => prend directement les temperatures |
---|
366 | !du sol d'ORCHIDEE |
---|
367 | |
---|
368 | ! pwater=-9 |
---|
369 | rwater(:) = pwater_wet2 |
---|
370 | !nwater = NINT(rwater) |
---|
371 | DO ipoint = 1,npts |
---|
372 | nwater(ipoint) = NINT(rwater(ipoint)) |
---|
373 | ENDDO |
---|
374 | |
---|
375 | |
---|
376 | substrat(:) = 0.0 |
---|
377 | DO i=1,11 |
---|
378 | DO ipoint = 1,npts |
---|
379 | substrat(ipoint) = substrat(ipoint) + veget_max(ipoint,i)*carbon_surf(ipoint,1,i) |
---|
380 | ENDDO |
---|
381 | ENDDO |
---|
382 | DO ipoint = 1,npts |
---|
383 | IF (sum_veget_nat(ipoint) .NE. 0.) THEN |
---|
384 | substrat(ipoint) = substrat(ipoint)/sum_veget_nat(ipoint) |
---|
385 | ENDIF |
---|
386 | ENDDO |
---|
387 | |
---|
388 | |
---|
389 | |
---|
390 | !impose valeurs par defaut pour eviter que le modele s arrete dans le cas d un grid-cell |
---|
391 | !où iroot,isoil et tveg seraient nuls |
---|
392 | WHERE ((veget_max(:,1) .GT. 0.95) .OR. (sum_veget_nat_ss_nu(:) .LE. 0.1)) |
---|
393 | iroot(:)=1*64*1/1 |
---|
394 | isoil(:)=129 |
---|
395 | tveg(:)=0.5*1/15. |
---|
396 | ENDWHERE |
---|
397 | |
---|
398 | !************************************ |
---|
399 | !**ANCIEN Scalc********************** |
---|
400 | !************************************ |
---|
401 | |
---|
402 | rcmax(:)=scmax+(ibare(:)/100.)*scmax |
---|
403 | |
---|
404 | forg(:,:)=0. |
---|
405 | |
---|
406 | DO i=1,ns |
---|
407 | DO ipoint = 1,npts |
---|
408 | nsoil=ns-isoil(ipoint) |
---|
409 | nroot=ns-iroot(ipoint) |
---|
410 | ihelp=(nroot-nsoil)+1 |
---|
411 | IF ((i .GE. nsoil) .AND. (i .LE. nroot)) THEN |
---|
412 | ihelp=ihelp-1 |
---|
413 | forg(ipoint,i)=EXP((-(ihelp*1.))/10.) |
---|
414 | ELSEIF (i .GE. nroot+1) THEN |
---|
415 | forg(ipoint,i)=1. |
---|
416 | ENDIF |
---|
417 | ENDDO |
---|
418 | ENDDO |
---|
419 | |
---|
420 | |
---|
421 | |
---|
422 | !**************ANCIEN Smodel |
---|
423 | |
---|
424 | !1 Initialisation |
---|
425 | ! |
---|
426 | ! |
---|
427 | ! |
---|
428 | |
---|
429 | IF ( firstcall ) THEN |
---|
430 | |
---|
431 | ! initialize tridiagonal coefficient matrix |
---|
432 | |
---|
433 | b(:,:)=0. |
---|
434 | a(:,:,:)=0. |
---|
435 | |
---|
436 | ! initialize fpart(i) (see below) |
---|
437 | ! Bunsen solubility coefficient (for air-water interface) |
---|
438 | ! (transformation coefficient for concentration<-->partial pressure |
---|
439 | |
---|
440 | fpart(:,:)=1. |
---|
441 | |
---|
442 | !partie du code initial qui initialisait le profil de concentration en CH4 |
---|
443 | !-> cf. stomate_io |
---|
444 | !dans stomate_io, je ne connais pas rcmax => remplace par scmax |
---|
445 | !!$ DO i=1,nvert |
---|
446 | !!$ DO ipoint = 1,npts |
---|
447 | !!$ mwater=ns+NINT(pwater) |
---|
448 | !!$ IF (i .LE. mwater) THEN |
---|
449 | !!$ uo(ipoint,i)=rcmax(ipoint) |
---|
450 | !!$ ELSE |
---|
451 | !!$ uo(ipoint,i)=catm |
---|
452 | !!$ ENDIF |
---|
453 | !!$ uold2(ipoint,i)=uo(ipoint,i) |
---|
454 | !!$ ENDDO |
---|
455 | !!$ ENDDO |
---|
456 | |
---|
457 | !mwater=ns+NINT(pwater) |
---|
458 | |
---|
459 | firstcall = .FALSE. |
---|
460 | ENDIF |
---|
461 | |
---|
462 | b(:,:)=0. |
---|
463 | a(:,:,:)=0. |
---|
464 | |
---|
465 | |
---|
466 | j_point(:)=0 |
---|
467 | DO i=1,ns |
---|
468 | DO ipoint = 1,npts |
---|
469 | nroot=ns-iroot(ipoint) |
---|
470 | IF (i .LE. nroot-1) THEN |
---|
471 | !(a): no plant transport below rooting depth nroot |
---|
472 | root(ipoint,i)=0. |
---|
473 | ELSEIF ((i .GE. nroot) .and. (iroot(ipoint) .EQ. 0)) THEN |
---|
474 | root(ipoint,i)=0 |
---|
475 | ELSEIF ((i .GE. nroot) .and. (iroot(ipoint) .NE. 0)) THEN |
---|
476 | ! (b): efficiency of plant-mediated transportlinearly increasing with |
---|
477 | ! decreasing soil depth (gets larger when going up) |
---|
478 | j_point(ipoint)=j_point(ipoint)+1 |
---|
479 | root(ipoint,i)=j_point(ipoint)*2./iroot(ipoint) |
---|
480 | ENDIF |
---|
481 | ENDDO |
---|
482 | ENDDO |
---|
483 | |
---|
484 | ! (c): no plant transport above soil surface |
---|
485 | DO i=ns+1,nvert |
---|
486 | DO ipoint = 1,npts |
---|
487 | root(ipoint,i)=0. |
---|
488 | ENDDO |
---|
489 | ENDDO |
---|
490 | |
---|
491 | !! define vertical root distribution root(ns) for plant-mediated transport |
---|
492 | !! (a): no plant transport below rooting depth nroot |
---|
493 | ! nroot=ns-iroot |
---|
494 | ! DO i=1,nroot-1 |
---|
495 | ! root(i)=0. |
---|
496 | ! ENDDO |
---|
497 | !! (b): efficiency of plant-mediated transportlinearly increasing with |
---|
498 | !! decreasing soil depth (gets larger when going up) |
---|
499 | ! j=0 |
---|
500 | ! DO i=nroot,ns |
---|
501 | ! IF (iroot .NE. 0) THEN |
---|
502 | ! j=j+1 |
---|
503 | ! root(i)=j*2./iroot |
---|
504 | ! ELSE |
---|
505 | ! root(i)=0. |
---|
506 | ! ENDIF |
---|
507 | ! ENDDO |
---|
508 | !! (c): no plant transport above soil surface |
---|
509 | ! DO i=ns+1,nvert |
---|
510 | ! root(i)=0. |
---|
511 | ! ENDDO |
---|
512 | |
---|
513 | !ancienne subroutine Sinter********************************************** |
---|
514 | ! interpolates Tsoil ==> value for each cm |
---|
515 | ! interpolate soil temperature from input soil layers to 1cm-thick |
---|
516 | ! model soil layers ==> sou(i) |
---|
517 | |
---|
518 | |
---|
519 | !plus d interpolation necessaire entre les temperatures du sol de differents niveaux mises en entree du modele de Walter |
---|
520 | !met directement les temperatures des differentes couches de sol d'ORCHIDEE |
---|
521 | !ATTENTION si dpu_cste varie!? => a changer |
---|
522 | DO i=1,75 |
---|
523 | DO ipoint = 1,npts |
---|
524 | sou(ipoint,i)=stempdiag(ipoint,10)-273.16 |
---|
525 | ENDDO |
---|
526 | ENDDO |
---|
527 | DO i=76,113 |
---|
528 | DO ipoint = 1,npts |
---|
529 | sou(ipoint,i)=stempdiag(ipoint,9)-273.16 |
---|
530 | ENDDO |
---|
531 | ENDDO |
---|
532 | DO i=114,131 |
---|
533 | DO ipoint = 1,npts |
---|
534 | sou(ipoint,i)=stempdiag(ipoint,8)-273.16 |
---|
535 | ENDDO |
---|
536 | ENDDO |
---|
537 | DO i=132,141 |
---|
538 | DO ipoint = 1,npts |
---|
539 | sou(ipoint,i)=stempdiag(ipoint,7)-273.16 |
---|
540 | ENDDO |
---|
541 | ENDDO |
---|
542 | DO i=142,146 |
---|
543 | DO ipoint = 1,npts |
---|
544 | sou(ipoint,i)=stempdiag(ipoint,6)-273.16 |
---|
545 | ENDDO |
---|
546 | ENDDO |
---|
547 | DO i=147,148 |
---|
548 | DO ipoint = 1,npts |
---|
549 | sou(ipoint,i)=stempdiag(ipoint,5)-273.16 |
---|
550 | ENDDO |
---|
551 | ENDDO |
---|
552 | |
---|
553 | sou(:,149)=stempdiag(:,4)-273.16 |
---|
554 | |
---|
555 | DO i=150,ns |
---|
556 | DO ipoint = 1,npts |
---|
557 | sou(ipoint,i)=stempdiag(ipoint,3)-273.16 |
---|
558 | ENDDO |
---|
559 | ENDDO |
---|
560 | |
---|
561 | |
---|
562 | |
---|
563 | !!$ DO i=1,71 |
---|
564 | !!$ DO ipoint = 1,npts |
---|
565 | !!$ sou(ipoint,i)=std(ipoint)+((193.+i)/265.)*(std5(ipoint)-std(ipoint)) |
---|
566 | !!$ ENDDO |
---|
567 | !!$ ENDDO |
---|
568 | !!$ |
---|
569 | !!$ j=0 |
---|
570 | !!$ DO i=72,137 |
---|
571 | !!$ j=j+1 |
---|
572 | !!$ DO ipoint = 1,npts |
---|
573 | !!$ sou(ipoint,i)=std5(ipoint)+((j-1.)/66.)*(std4(ipoint)-std5(ipoint)) |
---|
574 | !!$ ENDDO |
---|
575 | !!$ ENDDO |
---|
576 | !!$ |
---|
577 | !!$ j=0 |
---|
578 | !!$ DO i=138,149 |
---|
579 | !!$ j=j+1 |
---|
580 | !!$ DO ipoint = 1,npts |
---|
581 | !!$ sou(ipoint,i)=std4(ipoint)+((j-1.)/12.)*(std3(ipoint)-std4(ipoint)) |
---|
582 | !!$ ENDDO |
---|
583 | !!$ ENDDO |
---|
584 | !!$ |
---|
585 | !!$ DO i=150,ns |
---|
586 | !!$ DO ipoint = 1,npts |
---|
587 | !!$ sou(ipoint,i)=std3(ipoint) |
---|
588 | !!$ ENDDO |
---|
589 | !!$ ENDDO |
---|
590 | |
---|
591 | |
---|
592 | !!ancienne subroutine Sinter********************************************** |
---|
593 | !! interpolates Tsoil ==> value for each cm |
---|
594 | !! interpolate soil temperature from input soil layers to 1cm-thick |
---|
595 | !! model soil layers ==> sou(i) |
---|
596 | ! do i=1,71 |
---|
597 | ! sou(i)=std+((193.+i)/265.)*(std5-std) |
---|
598 | ! enddo |
---|
599 | ! j=0 |
---|
600 | ! do i=72,137 |
---|
601 | ! j=j+1 |
---|
602 | ! sou(i)=std5+((j-1.)/66.)*(std4-std5) |
---|
603 | ! enddo |
---|
604 | ! j=0 |
---|
605 | ! do i=138,149 |
---|
606 | ! j=j+1 |
---|
607 | ! sou(i)=std4+((j-1.)/12.)*(std3-std4) |
---|
608 | ! enddo |
---|
609 | ! do i=150,ns |
---|
610 | ! sou(i)=std3 |
---|
611 | ! enddo |
---|
612 | !!************************************************************************ |
---|
613 | |
---|
614 | !utilise le LAI calcule par ORCHIDEE plutot que de calculer son propre LAI comme dans le modele |
---|
615 | !de Walter initial |
---|
616 | |
---|
617 | ! define fgrow: growing stage of plants (LAI, after (Dickinson, 1993)) |
---|
618 | ! different thresholds for regions with different annual mean temperature |
---|
619 | ! tmean |
---|
620 | ! if (tmean .le. 5.) then |
---|
621 | ! xgrow=2. |
---|
622 | ! else |
---|
623 | ! xgrow=7. |
---|
624 | ! endif |
---|
625 | ! if (sou(101) .le. xgrow) then |
---|
626 | ! fgrow=0.0001 |
---|
627 | ! else |
---|
628 | !! maximum for regions with large tmean |
---|
629 | ! if (sou(101) .ge. (xgrow+10.)) then |
---|
630 | ! fgrow=4. |
---|
631 | ! else |
---|
632 | ! fgrow=0.0001+3.999*(1-((((xgrow+10.)-sou(101))/10.)**2.)) |
---|
633 | ! endif |
---|
634 | ! endif |
---|
635 | |
---|
636 | !WRITE(*,*) 'veget', veget |
---|
637 | !WRITE(*,*) 'veget_max', veget_max |
---|
638 | |
---|
639 | fgrow(:) = 0.0 |
---|
640 | DO i=1,11 |
---|
641 | DO ipoint = 1,npts |
---|
642 | fgrow(ipoint) = fgrow(ipoint) + veget_max(ipoint,i)*lai(ipoint,i) |
---|
643 | ENDDO |
---|
644 | ENDDO |
---|
645 | DO ipoint = 1,npts |
---|
646 | IF (sum_veget_nat(ipoint) .NE. 0.) THEN |
---|
647 | fgrow(ipoint) = fgrow(ipoint)/sum_veget_nat(ipoint) |
---|
648 | ELSE |
---|
649 | fgrow(ipoint) = 0.0 |
---|
650 | ENDIF |
---|
651 | ENDDO |
---|
652 | |
---|
653 | |
---|
654 | ! calculation of n0,n1,n2 / determine coefficients of diffusion |
---|
655 | ! n0: lower boundary (soil depth or thaw depth) |
---|
656 | ! n1: water table (if water<soil surface), else soil surface |
---|
657 | ! n2: soil surface ( " " " ), " water table |
---|
658 | ! n3: upper boundary (fixed: n3=n2+4) |
---|
659 | ! global model: NO thaw depth (set nthaw=max soil depth [=150]) |
---|
660 | ! for 1-dim model: read thaw depth in Sread.f & comment line 'nthaw=150' out |
---|
661 | ! instead write: 'nthat=pthaw(itime)' |
---|
662 | !nthaw=150 |
---|
663 | nthaw=150 !modifie nthaw car j ai modiife ns dans stomate_cste_wetlands.f90 |
---|
664 | |
---|
665 | DO ipoint = 1,npts |
---|
666 | n0(ipoint)=ns-MIN(isoil(ipoint),nthaw) |
---|
667 | nw(ipoint)=MAX(n0(ipoint),ns+nwater(ipoint)) |
---|
668 | nw(ipoint)=MIN(nw(ipoint),(nvert-4)) |
---|
669 | ENDDO |
---|
670 | |
---|
671 | |
---|
672 | DO ipoint = 1,npts ! boucle la plus externe |
---|
673 | |
---|
674 | !le cas suivant n arrive jamais!!!: la WTD est toujours en-dessous ou au niveau du sol |
---|
675 | !au niveau de la vectorisation=>n2 vaut toujours ns |
---|
676 | ! if water table is above soil surface |
---|
677 | ! IF (nw(ipoint) .GT. ns) THEN |
---|
678 | ! n1(ipoint)=ns |
---|
679 | ! n2(ipoint)=nw(ipoint) |
---|
680 | ! diffdo(ipoint)=diffsw |
---|
681 | ! diffup(ipoint)=diffwa |
---|
682 | ! rvmax(ipoint)=0. |
---|
683 | ! ENDIF |
---|
684 | |
---|
685 | !!!psstry+ |
---|
686 | ! define diffusion coefficients (after Penman) |
---|
687 | ! soil air |
---|
688 | diffsa=diffair*0.66*rpv |
---|
689 | ! soil water |
---|
690 | diffsw=0.0001*diffsa |
---|
691 | ! standing water |
---|
692 | diffwa=0.0001*diffair |
---|
693 | !!!psstry- |
---|
694 | |
---|
695 | ! if water table is below soil surface |
---|
696 | IF (nw(ipoint) .LT. ns) THEN |
---|
697 | n1(ipoint)=nw(ipoint) |
---|
698 | n2=ns |
---|
699 | diffdo(ipoint)=diffsw |
---|
700 | diffup(ipoint)=diffsa |
---|
701 | rvmax(ipoint)=xvmax |
---|
702 | ENDIF |
---|
703 | |
---|
704 | !WRITE(*,*) 'n2','ipoint',ipoint, n2, n0(ipoint), ns, n1(ipoint), n3 |
---|
705 | |
---|
706 | !!$! if water table is at soil surface |
---|
707 | !!$ IF (nw(ipoint) .EQ. ns) THEN |
---|
708 | !!$ n1(ipoint)=nw(ipoint) |
---|
709 | !!$ n2=nw(ipoint) |
---|
710 | !!$ diffdo(ipoint)=diffsw |
---|
711 | !!$ diffup(ipoint)=diffsw |
---|
712 | !!$ rvmax(ipoint)=0. |
---|
713 | !!$ endif |
---|
714 | ENDDO |
---|
715 | |
---|
716 | !n3 est constant finalement car pour tous les pixels n2=ns et ns=171(ou151) |
---|
717 | n3=n2+4 |
---|
718 | |
---|
719 | !ijump(:)=1 |
---|
720 | DO ipoint = 1,npts |
---|
721 | ! make sure that n0<n1<n2 |
---|
722 | IF ((n0(ipoint) .EQ. n2) .OR. (n0(ipoint)+1 .EQ. n2)) THEN |
---|
723 | ijump(ipoint)=0 |
---|
724 | ELSE |
---|
725 | ijump(ipoint)=1 |
---|
726 | ENDIF |
---|
727 | IF ((n1(ipoint) .EQ. n0(ipoint)) .AND. (n2 .GT. n0(ipoint)+1)) THEN |
---|
728 | n1(ipoint)=n0(ipoint)+1 |
---|
729 | diffdo(ipoint)=diffup(ipoint) |
---|
730 | sou(ipoint,n0(ipoint))=0. |
---|
731 | sou(ipoint,n1(ipoint))=0. |
---|
732 | ENDIF |
---|
733 | ENDDO |
---|
734 | |
---|
735 | |
---|
736 | ! Bunsen solubility coefficient (air-water interface) |
---|
737 | ! (transformation coefficient for concentration<-->partial pressure |
---|
738 | ! transformtion in water (do i=1,nw) and air (do i=nw+1,n)) |
---|
739 | |
---|
740 | !inverse l ordre des boucles pour ameliorer le vectorissation |
---|
741 | DO i=1,nvert |
---|
742 | DO ipoint = 1,npts |
---|
743 | IF (i .LE. nw(ipoint))THEN |
---|
744 | fpart(ipoint,i)=1. |
---|
745 | ELSE |
---|
746 | fpart(ipoint,i)=1./23. |
---|
747 | ENDIF |
---|
748 | ENDDO |
---|
749 | ENDDO |
---|
750 | |
---|
751 | |
---|
752 | ! DO i=1,nw |
---|
753 | ! fpart(i)=1. |
---|
754 | ! ENDDO |
---|
755 | ! DO i=nw+1,nvert |
---|
756 | ! fpart(i)=1./23. |
---|
757 | ! ENDDO |
---|
758 | |
---|
759 | |
---|
760 | |
---|
761 | ! calculate total daily production [wpro] / oxidation rate [goxid], |
---|
762 | ! difference in conc(day)-conc(day-1) [tout] / |
---|
763 | ! daily amount of CH4 transported through plants [dplant] |
---|
764 | ! necessary for calculation of diffusive flux from mass balance |
---|
765 | DO ipoint = 1,npts |
---|
766 | wpro(ipoint)=0. |
---|
767 | goxid(ipoint)=0. |
---|
768 | tout(ipoint)=0. |
---|
769 | dplant(ipoint)=0. |
---|
770 | negconc(ipoint)=0. |
---|
771 | ENDDO |
---|
772 | |
---|
773 | ! if water table falls (i.e. nw (new water tabel) < nwtm1 (old water table)): |
---|
774 | wtdiff=0. |
---|
775 | !supprime boucle if sur ijump: dans l approche choisie WTD=cste au cours du temps |
---|
776 | !perform calculation only if n0<n1<n2 |
---|
777 | !!! IF (ijump .EQ. 1) THEN |
---|
778 | |
---|
779 | ! update old u-values after change of water table/thaw depth |
---|
780 | ! water-air interface: conc ==> partial pressure |
---|
781 | ! transform from concentration into partial pressure unit |
---|
782 | |
---|
783 | DO i=1,n3 |
---|
784 | DO ipoint = 1,npts |
---|
785 | IF (i .GE. n0(ipoint)) THEN |
---|
786 | uo(ipoint,i)=uo(ipoint,i)*fpart(ipoint,i) |
---|
787 | ENDIF |
---|
788 | ENDDO |
---|
789 | ENDDO |
---|
790 | |
---|
791 | |
---|
792 | |
---|
793 | ! |
---|
794 | ! solve diffusion equation (Crank-Nicholson) |
---|
795 | ! Schwarz, Numerische Mathematik, Teubner, Stuttgart, p.476) |
---|
796 | ! Press et al., Numerical Recepies, FORTRAN, p. 838) |
---|
797 | !************************************************************* |
---|
798 | |
---|
799 | ! left hand side coefficient matrix: |
---|
800 | !************************************ |
---|
801 | |
---|
802 | ! left boundary |
---|
803 | |
---|
804 | |
---|
805 | p=0. |
---|
806 | DO ipoint = 1,npts |
---|
807 | amhalf=diffdo(ipoint) |
---|
808 | aphalf=diffdo(ipoint) |
---|
809 | a(ipoint,n0(ipoint),n0(ipoint))=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
810 | a(ipoint,n0(ipoint),n0(ipoint)+1)=-rkh*(aphalf+amhalf) |
---|
811 | ENDDO |
---|
812 | |
---|
813 | |
---|
814 | ! inner points |
---|
815 | DO ipoint = 1,npts |
---|
816 | j1(ipoint)=n0(ipoint)-1 |
---|
817 | j2(ipoint)=j1(ipoint)+1 |
---|
818 | j3(ipoint)=j2(ipoint)+1 |
---|
819 | ENDDO |
---|
820 | |
---|
821 | DO i=1,n2 !idem que precedement: sur-calcul |
---|
822 | DO ipoint = 1,npts |
---|
823 | IF ((i .GE. n0(ipoint)+1) .AND. (i .LE. n1(ipoint)-1)) THEN |
---|
824 | j1(ipoint)=j1(ipoint)+1 |
---|
825 | j2(ipoint)=j2(ipoint)+1 |
---|
826 | j3(ipoint)=j3(ipoint)+1 |
---|
827 | amhalf=diffdo(ipoint) |
---|
828 | aphalf=diffdo(ipoint) |
---|
829 | p=0. |
---|
830 | a(ipoint, i,j1(ipoint))=-rkh*amhalf |
---|
831 | a(ipoint,i,j2(ipoint))=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
832 | a(ipoint,i,j3(ipoint))=-rkh*aphalf |
---|
833 | ENDIF |
---|
834 | ENDDO |
---|
835 | ENDDO |
---|
836 | |
---|
837 | ! j1=n0-1 |
---|
838 | ! j2=j1+1 |
---|
839 | ! j3=j2+1 |
---|
840 | ! do i=n0+1,n1-1 |
---|
841 | ! j1=j1+1 |
---|
842 | ! j2=j2+1 |
---|
843 | ! j3=j3+1 |
---|
844 | ! amhalf=diffdo |
---|
845 | ! aphalf=diffdo |
---|
846 | ! p=0. |
---|
847 | ! a(i,j1)=-rkh*amhalf |
---|
848 | ! a(i,j2)=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
849 | ! a(i,j3)=-rkh*aphalf |
---|
850 | ! enddo |
---|
851 | |
---|
852 | |
---|
853 | diffgrdo=0.57*0.66*rpv |
---|
854 | diffgrup=0.0248*0.66*rpv |
---|
855 | |
---|
856 | !les parties suivantes sont differentes entre stomate_wet_ch4_pt_ter_0 et stomate_wet_ch4_pt_ter_m10 |
---|
857 | |
---|
858 | DO i = 1,n3-1 |
---|
859 | !sur-calcul: au depart j avais = n0,n3-1 |
---|
860 | !mais les nouvelles bornes choisies permettent d ameliorer le vectorisation |
---|
861 | DO ipoint = 1,npts |
---|
862 | |
---|
863 | !1er CAS: (nw(ipoint) .LT. ns-1) |
---|
864 | !diffusion coefficients water-air INTERFACE (soil water --> soil air) |
---|
865 | !par defaut ici on a toujours (nw(ipoint) .LT. ns-1) |
---|
866 | |
---|
867 | IF (i .EQ. n1(ipoint)) THEN |
---|
868 | !i=n1 |
---|
869 | j1(ipoint)=n1(ipoint)-1 |
---|
870 | j2(ipoint)=n1(ipoint) |
---|
871 | j3(ipoint)=n1(ipoint)+1 |
---|
872 | amhalf=diffdo(ipoint) |
---|
873 | aphalf=diffgrdo |
---|
874 | p=0. |
---|
875 | a(ipoint,i,j1(ipoint))=-rkh*amhalf |
---|
876 | a(ipoint,i,j2(ipoint))=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
877 | a(ipoint,i,j3(ipoint))=-rkh*aphalf |
---|
878 | ELSEIF (i .EQ. n1(ipoint)+1) THEN |
---|
879 | !i=n1+1 |
---|
880 | j1(ipoint)=n1(ipoint) |
---|
881 | j2(ipoint)=n1(ipoint)+1 |
---|
882 | j3(ipoint)=n1(ipoint)+2 |
---|
883 | amhalf=diffgrup |
---|
884 | aphalf=diffup(ipoint) |
---|
885 | p=0. |
---|
886 | a(ipoint,i,j1(ipoint))=-rkh*amhalf |
---|
887 | a(ipoint,i,j2(ipoint))=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
888 | a(ipoint,i,j3(ipoint))=-rkh*aphalf |
---|
889 | ELSEIF ((i .GE. n1(ipoint)+2) .AND. (i .LE. n2-1)) THEN |
---|
890 | !DO i=n1+2,n2-1 |
---|
891 | j1(ipoint)=j1(ipoint)+1 |
---|
892 | j2(ipoint)=j2(ipoint)+1 |
---|
893 | j3(ipoint)=j3(ipoint)+1 |
---|
894 | amhalf=diffup(ipoint) |
---|
895 | aphalf=diffup(ipoint) |
---|
896 | p=0. |
---|
897 | a(ipoint,i,j1(ipoint))=-rkh*amhalf |
---|
898 | a(ipoint,i,j2(ipoint))=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
899 | a(ipoint,i,j3(ipoint))=-rkh*aphalf |
---|
900 | !ENDDO |
---|
901 | ELSEIF (i .EQ. n2) THEN |
---|
902 | !i=n2 |
---|
903 | j1(ipoint)=n2-1 |
---|
904 | j2(ipoint)=n2 |
---|
905 | j3(ipoint)=n2+1 |
---|
906 | amhalf=diffup(ipoint) |
---|
907 | aphalf=diffair |
---|
908 | p=0. |
---|
909 | a(ipoint,i,j1(ipoint))=-rkh*amhalf |
---|
910 | a(ipoint,i,j2(ipoint))=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
911 | a(ipoint,i,j3(ipoint))=-rkh*aphalf |
---|
912 | ELSEIF ((i .GE. n2+1) .AND. (i .LE. n3-1)) THEN |
---|
913 | !DO i=n2+1,n3-1 |
---|
914 | j1(ipoint)=j1(ipoint)+1 |
---|
915 | j2(ipoint)=j2(ipoint)+1 |
---|
916 | j3(ipoint)=j3(ipoint)+1 |
---|
917 | amhalf=diffair |
---|
918 | aphalf=diffair |
---|
919 | p=0. |
---|
920 | a(ipoint,i,j1(ipoint))=-rkh*amhalf |
---|
921 | a(ipoint,i,j2(ipoint))=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
922 | a(ipoint,i,j3(ipoint))=-rkh*aphalf |
---|
923 | !ENDDO |
---|
924 | |
---|
925 | !!$ !2eme CAS: (nw .EQ. ns-1) |
---|
926 | !!$ ! diffusion coefficients water-air interface (soil water --> soil air) |
---|
927 | !!$ ELSEIF ((nw(ipoint) .EQ. ns-1) .AND. (i .EQ. n1(ipoint))) THEN |
---|
928 | !!$ !i=n1 |
---|
929 | !!$ j1(ipoint)=n1(ipoint)-1 |
---|
930 | !!$ j2(ipoint)=n1(ipoint) |
---|
931 | !!$ j3(ipoint)=n1(ipoint)+1 |
---|
932 | !!$ amhalf=diffdo(ipoint) |
---|
933 | !!$ aphalf=diffgrdo |
---|
934 | !!$ p=0. |
---|
935 | !!$ a(ipoint,i,j1(ipoint))=-rkh*amhalf |
---|
936 | !!$ a(ipoint,i,j2(ipoint))=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
937 | !!$ a(ipoint,i,j3(ipoint))=-rkh*aphalf |
---|
938 | !!$ ELSEIF ((nw(ipoint) .EQ. ns-1) .AND. (i .EQ. n2)) THEN |
---|
939 | !!$ !i=n2 |
---|
940 | !!$ j1(ipoint)=n2-1 |
---|
941 | !!$ j2(ipoint)=n2 |
---|
942 | !!$ j3(ipoint)=n2+1 |
---|
943 | !!$ amhalf=diffgrup |
---|
944 | !!$ aphalf=diffup(ipoint) |
---|
945 | !!$ p=0. |
---|
946 | !!$ a(ipoint,i,j1(ipoint))=-rkh*amhalf |
---|
947 | !!$ a(ipoint,i,j2(ipoint))=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
948 | !!$ a(ipoint,i,j3(ipoint))=-rkh*aphalf |
---|
949 | !!$ ELSEIF ((nw(ipoint) .EQ. ns-1) .AND. (i .EQ. n2+1)) THEN |
---|
950 | !!$ !i=n2+1 |
---|
951 | !!$ j1(ipoint)=n2 |
---|
952 | !!$ j2(ipoint)=n2+1 |
---|
953 | !!$ j3(ipoint)=n2+2 |
---|
954 | !!$ amhalf=diffup(ipoint) |
---|
955 | !!$ aphalf=diffair |
---|
956 | !!$ p=0. |
---|
957 | !!$ a(ipoint,i,j1(ipoint))=-rkh*amhalf |
---|
958 | !!$ a(ipoint,i,j2(ipoint))=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
959 | !!$ a(ipoint,i,j3(ipoint))=-rkh*aphalf |
---|
960 | !!$ ELSEIF ((nw(ipoint) .EQ. ns-1) .AND. (i .GE. n2+2) .AND. (i .LE. n3-1)) THEN |
---|
961 | !!$ !DO i=n2+2,n3-1 |
---|
962 | !!$ j1(ipoint)=j1(ipoint)+1 |
---|
963 | !!$ j2(ipoint)=j2(ipoint)+1 |
---|
964 | !!$ j3(ipoint)=j3(ipoint)+1 |
---|
965 | !!$ amhalf=diffair |
---|
966 | !!$ aphalf=diffair |
---|
967 | !!$ p=0. |
---|
968 | !!$ a(ipoint,i,j1(ipoint))=-rkh*amhalf |
---|
969 | !!$ a(ipoint,i,j2(ipoint))=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
970 | !!$ a(ipoint,i,j3(ipoint))=-rkh*aphalf |
---|
971 | !!$ !enddo |
---|
972 | !!$ |
---|
973 | !!$ !3eme CAS: (nw .EQ. ns) |
---|
974 | !!$ ! diffusion coefficients water-air interface (soil water --> soil air) |
---|
975 | !!$ ELSEIF ((nw(ipoint) .EQ. ns) .AND. (i .EQ. n1(ipoint))) THEN |
---|
976 | !!$ !i=n1 |
---|
977 | !!$ j1(ipoint)=n1(ipoint)-1 |
---|
978 | !!$ j2(ipoint)=n1(ipoint) |
---|
979 | !!$ j3(ipoint)=n1(ipoint)+1 |
---|
980 | !!$ amhalf=diffdo(ipoint) |
---|
981 | !!$ aphalf=diffgrdo |
---|
982 | !!$ p=0. |
---|
983 | !!$ a(ipoint,i,j1(ipoint))=-rkh*amhalf |
---|
984 | !!$ a(ipoint,i,j2(ipoint))=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
985 | !!$ a(ipoint,i,j3(ipoint))=-rkh*aphalf |
---|
986 | !!$ ELSEIF ((nw(ipoint) .EQ. ns) .AND. (i .EQ. n1(ipoint)+1)) THEN |
---|
987 | !!$ !i=n1+1 |
---|
988 | !!$ j1(ipoint)=n1(ipoint) |
---|
989 | !!$ j2(ipoint)=n1(ipoint)+1 |
---|
990 | !!$ j3(ipoint)=n1(ipoint)+2 |
---|
991 | !!$ amhalf=diffgrup |
---|
992 | !!$ aphalf=diffair |
---|
993 | !!$ p=0. |
---|
994 | !!$ a(ipoint,i,j1(ipoint))=-rkh*amhalf |
---|
995 | !!$ a(ipoint,i,j2(ipoint))=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
996 | !!$ a(ipoint,i,j3(ipoint))=-rkh*aphalf |
---|
997 | !!$ ELSEIF ((nw(ipoint) .EQ. ns) .AND. (i .GE. n1(ipoint)+2) .AND. (i .LE. n3-1)) THEN |
---|
998 | !!$ ! do i=n1+2,n3-1 |
---|
999 | !!$ j1(ipoint)=j1(ipoint)+1 |
---|
1000 | !!$ j2(ipoint)=j2(ipoint)+1 |
---|
1001 | !!$ j3(ipoint)=j3(ipoint)+1 |
---|
1002 | !!$ amhalf=diffair |
---|
1003 | !!$ aphalf=diffair |
---|
1004 | !!$ p=0. |
---|
1005 | !!$ a(ipoint,i,j1(ipoint))=-rkh*amhalf |
---|
1006 | !!$ a(ipoint,i,j2(ipoint))=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
1007 | !!$ a(ipoint,i,j3(ipoint))=-rkh*aphalf |
---|
1008 | !!$ !ENDDO |
---|
1009 | |
---|
1010 | !4eme CAS: (nw .GE. ns): non traite dans ORCHIDEE! la WTD n est jamais au-dessus de la surface |
---|
1011 | !de toute facon pas d oxydation dans ce cas dans le modele de Wlater (juste modif du transport) |
---|
1012 | ENDIF |
---|
1013 | ENDDO |
---|
1014 | ENDDO |
---|
1015 | |
---|
1016 | |
---|
1017 | ! right boundary |
---|
1018 | |
---|
1019 | DO ipoint = 1,npts |
---|
1020 | amhalf=diffair |
---|
1021 | aphalf=diffair |
---|
1022 | p=0. |
---|
1023 | a(ipoint,n3,n3-1)=-rkh*amhalf |
---|
1024 | a(ipoint,n3,n3)=2.+rkh*(amhalf+aphalf-h**2*p) |
---|
1025 | ENDDO |
---|
1026 | |
---|
1027 | |
---|
1028 | ! |
---|
1029 | ! begin of time loop (time steps per day) |
---|
1030 | !***************************************** |
---|
1031 | |
---|
1032 | DO iday=1,nday |
---|
1033 | |
---|
1034 | |
---|
1035 | |
---|
1036 | !***************************************** |
---|
1037 | |
---|
1038 | ! right hand side coefficient matrix: |
---|
1039 | !************************************ |
---|
1040 | |
---|
1041 | ! left boundary |
---|
1042 | |
---|
1043 | DO ipoint = 1,npts |
---|
1044 | i=n0(ipoint) |
---|
1045 | amhalf=diffdo(ipoint) |
---|
1046 | aphalf=diffdo(ipoint) |
---|
1047 | p=0. |
---|
1048 | ! q10 dependent production rate (no production if T<0oC, i.e.sou(i)<=0oC) |
---|
1049 | ! forg(i): vertical distribution of substrate in soil (calculated in Scalc.f) |
---|
1050 | ! fin(itime): describes seasonality of NPP |
---|
1051 | ! rq10: Q10 value of production rate (from para.h) |
---|
1052 | ! r0pl: tuning parameter sr0pl (scaled in Scalc.f) |
---|
1053 | IF (sou(ipoint,n0(ipoint)) .GT. 0.) THEN |
---|
1054 | rexp=(sou(ipoint,n0(ipoint))-MAX(0.,tmean(ipoint)))/10. |
---|
1055 | ! source=forg(n0)*r0pl*(rq10**rexp)*fin(itime) |
---|
1056 | ! source=forg(n0)*(rq10**rexp)*carbon(itime)*(0.001239567) |
---|
1057 | ! source=forg(n0)*(rq10**rexp)*carbon(itime)*r0pl*(0.0012) |
---|
1058 | ! source=5.0*(rq10**rexp) |
---|
1059 | source=forg(ipoint,n0(ipoint))*(rq10**rexp)*substrat(ipoint)*alpha_PFT(ipoint) |
---|
1060 | ELSE |
---|
1061 | source=0. |
---|
1062 | ENDIF |
---|
1063 | q=source |
---|
1064 | wpro(ipoint)=wpro(ipoint)+q |
---|
1065 | ! transform from concentration into partial pressure unit |
---|
1066 | q=q*fpart(ipoint,i) |
---|
1067 | aold1=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1068 | aold2=rkh*(aphalf+amhalf) |
---|
1069 | aold3=2.*rk*q |
---|
1070 | b(ipoint,n0(ipoint))=aold1*uo(ipoint,n0(ipoint))+aold2*uo(ipoint,n0(ipoint)+1)+aold3 |
---|
1071 | ENDDO |
---|
1072 | |
---|
1073 | !DO ipoint =1,npts |
---|
1074 | ! IF (iday .EQ. 2) THEN |
---|
1075 | ! WRITE(*,*) 'matrice b apres left boundary',b(ipoint,:) |
---|
1076 | ! ENDIF |
---|
1077 | !ENDDO |
---|
1078 | |
---|
1079 | |
---|
1080 | ! inner points (from left to right) |
---|
1081 | |
---|
1082 | !do i=n0+1,n1-1 |
---|
1083 | DO i = 1,n2 !sur calcul |
---|
1084 | DO ipoint = 1,npts |
---|
1085 | IF ((i .GE. n0(ipoint)+1) .AND. (i .LE. n1(ipoint)-1)) THEN |
---|
1086 | amhalf=diffdo(ipoint) |
---|
1087 | aphalf=diffdo(ipoint) |
---|
1088 | p=0. |
---|
1089 | ! q10 dependent production rate (no production if T<0oC, i.e.sou(i)<=0oC) |
---|
1090 | ! forg(i): vertical distribution of substrate in soil (calculated in Scalc.f) |
---|
1091 | ! fin(itime): describes seasonality of NPP |
---|
1092 | ! rq10: Q10 value of production rate (from para.h) |
---|
1093 | ! r0pl: tuning parameter sr0pl (scaled in Scalc.f) |
---|
1094 | IF (sou(ipoint,i) .GT. 0.) THEN |
---|
1095 | rexp=(sou(ipoint,i)-MAX(0.,tmean(ipoint)))/10. |
---|
1096 | ! source=forg(i)*r0pl*(rq10**rexp)*fin(itime) |
---|
1097 | ! source=forg(i)*(rq10**rexp)*carbon(itime)*(0.001239567) |
---|
1098 | !modifie le calcul de la production --> utilise le carbone comme proxi du substrat |
---|
1099 | source=forg(ipoint,i)*(rq10**rexp)*substrat(ipoint)*alpha_PFT(ipoint) |
---|
1100 | ELSE |
---|
1101 | source=0. |
---|
1102 | ENDIF |
---|
1103 | q=source |
---|
1104 | wpro(ipoint)=wpro(ipoint)+q |
---|
1105 | ! transform from concentration into partial pressure unit |
---|
1106 | q=q*fpart(ipoint,i) |
---|
1107 | aold1=rkh*amhalf |
---|
1108 | aold2=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1109 | aold3=rkh*aphalf |
---|
1110 | aold4=2.*rk*q |
---|
1111 | b(ipoint,i)=aold1*uo(ipoint,i-1)+aold2*uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4 |
---|
1112 | ENDIF |
---|
1113 | ENDDO |
---|
1114 | ENDDO |
---|
1115 | |
---|
1116 | !DO ipoint =1,npts |
---|
1117 | ! IF (iday .EQ. 2) THEN |
---|
1118 | ! WRITE(*,*) 'matrice b apres inner point avant cas',b(ipoint,:) |
---|
1119 | ! ENDIF |
---|
1120 | !ENDDO |
---|
1121 | |
---|
1122 | diffgrdo=0.57*0.66*rpv |
---|
1123 | diffgrup=0.0248*0.66*rpv |
---|
1124 | |
---|
1125 | !toujours vrai ici: nw(ipoint) .LT. ns-1) |
---|
1126 | |
---|
1127 | DO i = 1,n3-1 !sur-calcul: au depart j avais i = n0,n3-1 |
---|
1128 | DO ipoint = 1,npts |
---|
1129 | |
---|
1130 | !1er CAS: (nw(ipoint) .LT. ns-1) |
---|
1131 | !diffusion coefficients water-air INTERFACE (soil water --> soil air) |
---|
1132 | !!$ IF ((nw(ipoint) .LT. ns-1) .AND. (i .EQ. n1(ipoint))) THEN |
---|
1133 | |
---|
1134 | IF (i .EQ. n1(ipoint)) THEN |
---|
1135 | ! diffusion coefficients water-air interface (soil water --> soil air) |
---|
1136 | !i=n1 |
---|
1137 | amhalf=diffdo(ipoint) |
---|
1138 | aphalf=diffgrdo |
---|
1139 | p=0. |
---|
1140 | ! q10 dependent production rate (no production if T<0oC, i.e.sou(i)<=0oC) |
---|
1141 | ! forg(i): vertical distribution of substrate in soil (calculated in Scalc.f) |
---|
1142 | ! fin(itime): describes seasonality of NPP |
---|
1143 | ! rq10: Q10 value of production rate (from para.h) |
---|
1144 | ! r0pl: tuning parameter sr0pl (scaled in Scalc.f) |
---|
1145 | IF (sou(ipoint,i) .GT. 0.) THEN |
---|
1146 | rexp=(sou(ipoint,i)-MAX(0.,tmean(ipoint)))/10. |
---|
1147 | !source=forg(i)*r0pl*(rq10**rexp)*fin(itime) |
---|
1148 | source=forg(ipoint,i)*(rq10**rexp)*substrat(ipoint)*alpha_PFT(ipoint) |
---|
1149 | ELSE |
---|
1150 | source=0. |
---|
1151 | ENDIF |
---|
1152 | q=source |
---|
1153 | wpro(ipoint)=wpro(ipoint)+q |
---|
1154 | ! transform from concentration into partial pressure unit |
---|
1155 | q=q*fpart(ipoint,i) |
---|
1156 | aold1=rkh*amhalf |
---|
1157 | aold2=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1158 | aold3=rkh*aphalf |
---|
1159 | aold4=2.*rk*q |
---|
1160 | b(ipoint,i)=aold1*uo(ipoint,i-1)+aold2*uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4 |
---|
1161 | |
---|
1162 | |
---|
1163 | ELSEIF (i .EQ. n1(ipoint)+1) THEN |
---|
1164 | !i=n1+1 |
---|
1165 | amhalf=diffgrup |
---|
1166 | aphalf=diffup(ipoint) |
---|
1167 | p=0. |
---|
1168 | ! oxidation following Michaelis-Menten kinetics |
---|
1169 | IF (uo(ipoint,i) .GT. 0.) THEN |
---|
1170 | ! transform from partial pressure into concentration unit |
---|
1171 | uo(ipoint,i)=uo(ipoint,i)/fpart(ipoint,i) |
---|
1172 | ! Michaelis Menten Equation |
---|
1173 | q=-rvmax(ipoint)*uo(ipoint,i)/(uo(ipoint,i)+rkm) |
---|
1174 | ! temperature dependence (oxq10 from para.h) |
---|
1175 | q=(q/oxq10)*(oxq10**((sou(ipoint,i)-MAX(0.,tmean(ipoint)))/10.)) |
---|
1176 | goxid(ipoint)=goxid(ipoint)+q |
---|
1177 | ! transform from concentration into partial pressure unit |
---|
1178 | uo(ipoint,i)=uo(ipoint,i)*fpart(ipoint,i) |
---|
1179 | ! in case of negative concentrations: no oxidation |
---|
1180 | ! (can occur if number of time steps per day (nday) is small and water |
---|
1181 | ! table varies strongly) (help: increase nday (or change numerical scheme!)) |
---|
1182 | ELSE |
---|
1183 | q=0. |
---|
1184 | goxid(ipoint)=goxid(ipoint)+q |
---|
1185 | ENDIF |
---|
1186 | ! transform from concentration into partial pressure unit |
---|
1187 | q=q*fpart(ipoint,i) |
---|
1188 | aold1=rkh*amhalf |
---|
1189 | aold2=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1190 | aold3=rkh*aphalf |
---|
1191 | aold4=2.*rk*q |
---|
1192 | !!!!pss debug+ |
---|
1193 | ! write(2, *) 'pssdebug, i, ipoint, q, fpart(ipoint, i), rkh, amhalf, aphalf, h, p, rk, aold1, aold2, aold3, aold4, uo(ipoint, i), uo(ipoint, i-1), uo(ipoint, i+1)', i, ipoint, q, fpart(ipoint, i), rkh, amhalf, aphalf, h, p, rk, aold1, aold2, aold3, aold4, uo(ipoint, i), uo(ipoint, i-1), uo(ipoint, i+1) |
---|
1194 | b(ipoint,i)=aold1*uo(ipoint,i-1)+aold2*uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4 |
---|
1195 | |
---|
1196 | |
---|
1197 | ELSEIF ((i .GE. n1(ipoint)+2) .AND. (i .LE. n2-1)) THEN |
---|
1198 | !do i=n1+2,n2-1 |
---|
1199 | amhalf=diffup(ipoint) |
---|
1200 | aphalf=diffup(ipoint) |
---|
1201 | p=0. |
---|
1202 | ! oxidation following Michaelis-Menten kinetics |
---|
1203 | IF (uo(ipoint,i) .GT. 0.) THEN |
---|
1204 | ! transform from partial pressure into concentration unit |
---|
1205 | uo(ipoint,i)=uo(ipoint,i)/fpart(ipoint,i) |
---|
1206 | ! Michaelis Menten Equation |
---|
1207 | q=-rvmax(ipoint)*uo(ipoint,i)/(uo(ipoint,i)+rkm) |
---|
1208 | ! temperature dependence (oxq10 from para.h) |
---|
1209 | q=(q/oxq10)*(oxq10**((sou(ipoint,i)-MAX(0.,tmean(ipoint)))/10.)) |
---|
1210 | goxid(ipoint)=goxid(ipoint)+q |
---|
1211 | ! transform from concentration into partial pressure unit |
---|
1212 | uo(ipoint,i)=uo(ipoint,i)*fpart(ipoint,i) |
---|
1213 | ! in case of negative concentrations: no oxidation |
---|
1214 | ! (can occur if number of time steps per day (nday) is small and water |
---|
1215 | ! table varies strongly) (help: increase nday (or change numerical scheme!)) |
---|
1216 | ELSE |
---|
1217 | q=0. |
---|
1218 | goxid(ipoint)=goxid(ipoint)+q |
---|
1219 | ENDIF |
---|
1220 | ! transform from concentration into partial pressure unit |
---|
1221 | q=q*fpart(ipoint,i) |
---|
1222 | aold1=rkh*amhalf |
---|
1223 | aold2=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1224 | aold3=rkh*aphalf |
---|
1225 | aold4=2.*rk*q |
---|
1226 | b(ipoint,i)=aold1*uo(ipoint,i-1)+aold2*uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4 |
---|
1227 | |
---|
1228 | ELSEIF (i .EQ. n2) THEN |
---|
1229 | !i=n2 |
---|
1230 | amhalf=diffup(ipoint) |
---|
1231 | aphalf=diffair |
---|
1232 | p=0. |
---|
1233 | ! oxidation following Michaelis-Menten kinetics |
---|
1234 | IF (uo(ipoint,i) .GT. 0.) THEN |
---|
1235 | ! transform from partial pressure into concentration unit |
---|
1236 | uo(ipoint,i)=uo(ipoint,i)/fpart(ipoint,i) |
---|
1237 | ! Michaelis Menten Equation |
---|
1238 | q=-rvmax(ipoint)*uo(ipoint,i)/(uo(ipoint,i)+rkm) |
---|
1239 | ! temperature dependence (oxq10 from para.h) |
---|
1240 | q=(q/oxq10)*(oxq10**((sou(ipoint,i)-MAX(0.,tmean(ipoint)))/10.)) |
---|
1241 | goxid(ipoint)=goxid(ipoint)+q |
---|
1242 | ! transform from concentration into partial pressure unit |
---|
1243 | uo(ipoint,i)=uo(ipoint,i)*fpart(ipoint,i) |
---|
1244 | ! in case of negative concentrations: no oxidation |
---|
1245 | ! (can occur if number of time steps per day (nday) is small and water |
---|
1246 | ! table varies strongly) (help: increase nday (or change numerical scheme!)) |
---|
1247 | ELSE |
---|
1248 | q=0. |
---|
1249 | goxid(ipoint)=goxid(ipoint)+q |
---|
1250 | ENDIF |
---|
1251 | ! transform from concentration into partial pressure unit |
---|
1252 | q=q*fpart(ipoint,i) |
---|
1253 | aold1=rkh*amhalf |
---|
1254 | aold2=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1255 | aold3=rkh*aphalf |
---|
1256 | aold4=2.*rk*q |
---|
1257 | b(ipoint,i)=aold1*uo(ipoint,i-1)+aold2*uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4 |
---|
1258 | |
---|
1259 | ELSEIF ((i .GE. n2+1) .AND. (i .LE. n3-1)) THEN |
---|
1260 | !DO i=n2+1,n3-1 |
---|
1261 | amhalf=diffair |
---|
1262 | aphalf=diffair |
---|
1263 | p=0. |
---|
1264 | q=0. |
---|
1265 | aold1=rkh*amhalf |
---|
1266 | aold2=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1267 | aold3=rkh*aphalf |
---|
1268 | aold4=2.*rk*q |
---|
1269 | b(ipoint,i)=aold1*uo(ipoint,i-1)+aold2*uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4 |
---|
1270 | |
---|
1271 | ENDIF |
---|
1272 | ENDDO |
---|
1273 | ENDDO |
---|
1274 | |
---|
1275 | !!$ |
---|
1276 | !!$ |
---|
1277 | !!$ DO i = 1,n3-1 !sur-calcul: au depart j avais mis i = n0,n3-1 |
---|
1278 | !!$ DO ipoint = 1,npts |
---|
1279 | !!$ |
---|
1280 | !!$ !2nd CAS: (nw(ipoint) .EQ. ns) |
---|
1281 | !!$ !diffusion coefficients water-air INTERFACE (soil water --> soil air) |
---|
1282 | !!$ IF ((nw(ipoint) .EQ. ns) .AND. (i .EQ. n1(ipoint))) THEN |
---|
1283 | !!$ !i=n1 |
---|
1284 | !!$ amhalf=diffdo(ipoint) |
---|
1285 | !!$ aphalf=diffgrdo |
---|
1286 | !!$ p=0. |
---|
1287 | !!$! q10 dependent production rate (no production if T<0oC, i.e.sou(i)<=0oC) |
---|
1288 | !!$! forg(i): vertical distribution of substrate in soil (calculated in Scalc.f) |
---|
1289 | !!$! fin(itime): describes seasonality of NPP |
---|
1290 | !!$! rq10: Q10 value of production rate (from para.h) |
---|
1291 | !!$! r0pl: tuning parameter sr0pl (scaled in Scalc.f) |
---|
1292 | !!$ IF (sou(ipoint,i) .GT. 0.) THEN |
---|
1293 | !!$ rexp=(sou(ipoint,i)-MAX(0.,tmean(ipoint)))/10. |
---|
1294 | !!$! source=forg(i)*r0pl*(rq10**rexp)*fin(itime) |
---|
1295 | !!$! source=forg(i)*(rq10**rexp)*carbon(itime)*(0.001239567) |
---|
1296 | !!$! source=5.0*(rq10**rexp) |
---|
1297 | !!$ source=forg(ipoint,i)*(rq10**rexp)*substrat(ipoint)*alpha_CH4 |
---|
1298 | !!$ ELSE |
---|
1299 | !!$ source=0. |
---|
1300 | !!$ ENDIF |
---|
1301 | !!$ q=source |
---|
1302 | !!$ wpro(ipoint)=wpro(ipoint)+q |
---|
1303 | !!$! transform from concentration into partial pressure unit |
---|
1304 | !!$ q=q*fpart(ipoint,i) |
---|
1305 | !!$ aold1=rkh*amhalf |
---|
1306 | !!$ aold2=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1307 | !!$ aold3=rkh*aphalf |
---|
1308 | !!$ aold4=2.*rk*q |
---|
1309 | !!$ b(ipoint,i)=aold1*uo(ipoint,i-1)+aold2*uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4 |
---|
1310 | !!$ |
---|
1311 | !!$ ELSEIF ((nw(ipoint) .EQ. ns) .AND. (i .EQ. n1(ipoint)+1)) THEN |
---|
1312 | !!$ !i=n1+1 |
---|
1313 | !!$ amhalf=diffgrup |
---|
1314 | !!$ aphalf=diffair |
---|
1315 | !!$ p=0. |
---|
1316 | !!$ q=0. |
---|
1317 | !!$ aold1=rkh*amhalf |
---|
1318 | !!$ aold2=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1319 | !!$ aold3=rkh*aphalf |
---|
1320 | !!$ aold4=2.*rk*q |
---|
1321 | !!$ b(ipoint,i)=aold1*uo(ipoint,i-1)+aold2*uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4 |
---|
1322 | !!$ |
---|
1323 | !!$ ELSEIF ((nw(ipoint) .EQ. ns) .AND. (i .GE. n1(ipoint)+2) .AND. (i .LE. n3-1)) THEN |
---|
1324 | !!$ !DO i=n1+2,n3-1 |
---|
1325 | !!$ amhalf=diffair |
---|
1326 | !!$ aphalf=diffair |
---|
1327 | !!$ p=0. |
---|
1328 | !!$ q=0. |
---|
1329 | !!$ aold1=rkh*amhalf |
---|
1330 | !!$ aold2=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1331 | !!$ aold3=rkh*aphalf |
---|
1332 | !!$ aold4=2.*rk*q |
---|
1333 | !!$ b(ipoint,i)=aold1*uo(ipoint,i-1)+aold2*uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4 |
---|
1334 | !!$ !ENDDO |
---|
1335 | !!$ |
---|
1336 | !!$ ENDIF |
---|
1337 | !!$ ENDDO |
---|
1338 | !!$ ENDDO |
---|
1339 | !!$ |
---|
1340 | !!$ |
---|
1341 | |
---|
1342 | |
---|
1343 | |
---|
1344 | !!$ DO i = 1,n3-1 !sur-calcul: au depart j avais mis i = n0,n3-1 |
---|
1345 | !!$ DO ipoint = 1,npts |
---|
1346 | !!$ |
---|
1347 | !!$ !3eme CAS: (nw(ipoint) .EQ. ns-1) |
---|
1348 | !!$ !diffusion coefficients water-air INTERFACE (soil water --> soil air) |
---|
1349 | !!$ IF ((nw(ipoint) .EQ. ns-1) .AND. (i .EQ. n1(ipoint))) THEN |
---|
1350 | !!$ !i=n1 |
---|
1351 | !!$ amhalf=diffdo(ipoint) |
---|
1352 | !!$ aphalf=diffgrdo |
---|
1353 | !!$ p=0. |
---|
1354 | !!$! q10 dependent production rate (no production if T<0oC, i.e.sou(i)<=0oC) |
---|
1355 | !!$! forg(i): vertical distribution of substrate in soil (calculated in Scalc.f) |
---|
1356 | !!$! fin(itime): describes seasonality of NPP |
---|
1357 | !!$! rq10: Q10 value of production rate (from para.h) |
---|
1358 | !!$! r0pl: tuning parameter sr0pl (scaled in Scalc.f) |
---|
1359 | !!$ IF (sou(ipoint,i) .GT. 0.) THEN |
---|
1360 | !!$ rexp=(sou(ipoint,i)-MAX(0.,tmean(ipoint)))/10. |
---|
1361 | !!$ !source=forg(i)*r0pl*(rq10**rexp)*fin(itime) |
---|
1362 | !!$ source=forg(ipoint,i)*(rq10**rexp)*substrat(ipoint)*alpha_CH4 |
---|
1363 | !!$ ELSE |
---|
1364 | !!$ source=0. |
---|
1365 | !!$ ENDIF |
---|
1366 | !!$ q=source |
---|
1367 | !!$ wpro(ipoint)=wpro(ipoint)+q |
---|
1368 | !!$! transform from concentration into partial pressure unit |
---|
1369 | !!$ q=q*fpart(ipoint,i) |
---|
1370 | !!$ aold1=rkh*amhalf |
---|
1371 | !!$ aold2=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1372 | !!$ aold3=rkh*aphalf |
---|
1373 | !!$ aold4=2.*rk*q |
---|
1374 | !!$ b(ipoint,i)=aold1*uo(ipoint,i-1)+aold2*uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4 |
---|
1375 | !!$ |
---|
1376 | !!$ |
---|
1377 | !!$ ELSEIF ((nw(ipoint) .EQ. ns-1) .AND. (i .EQ. n2)) THEN |
---|
1378 | !!$ !i=n2 |
---|
1379 | !!$ amhalf=diffgrup |
---|
1380 | !!$ aphalf=diffup(ipoint) |
---|
1381 | !!$ p=0. |
---|
1382 | !!$! oxidation following Michaelis-Menten kinetics |
---|
1383 | !!$ IF (uo(ipoint,i) .GT. 0.) THEN |
---|
1384 | !!$! transform from partial pressure into concentration unit |
---|
1385 | !!$ uo(ipoint,i)=uo(ipoint,i)/fpart(ipoint,i) |
---|
1386 | !!$! Michaelis Menten Equation |
---|
1387 | !!$ q=-rvmax(ipoint)*uo(ipoint,i)/(uo(ipoint,i)+rkm) |
---|
1388 | !!$! temperature dependence (oxq10 from para.h) |
---|
1389 | !!$ q=(q/oxq10)*(oxq10**((sou(ipoint,i)-MAX(0.,tmean(ipoint)))/10.)) |
---|
1390 | !!$ goxid(ipoint)=goxid(ipoint)+q |
---|
1391 | !!$! transform from concentration into partial pressure unit |
---|
1392 | !!$ uo(ipoint,i)=uo(ipoint,i)*fpart(ipoint,i) |
---|
1393 | !!$! in case of negative concentrations: no oxidation |
---|
1394 | !!$! (can occur if number of time steps per day (nday) is small and water |
---|
1395 | !!$! table varies strongly) (help: increase nday (or change numerical scheme!)) |
---|
1396 | !!$ ELSE |
---|
1397 | !!$ q=0. |
---|
1398 | !!$ goxid(ipoint)=goxid(ipoint)+q |
---|
1399 | !!$ ENDIF |
---|
1400 | !!$ ! transform from concentration into partial pressure unit |
---|
1401 | !!$ q=q*fpart(ipoint,i) |
---|
1402 | !!$ aold1=rkh*amhalf |
---|
1403 | !!$ aold2=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1404 | !!$ aold3=rkh*aphalf |
---|
1405 | !!$ aold4=2.*rk*q |
---|
1406 | !!$ b(ipoint,i)=aold1*uo(ipoint,i-1)+aold2*uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4 |
---|
1407 | !!$ |
---|
1408 | !!$ ELSEIF ((nw(ipoint) .EQ. ns-1) .AND. (i .EQ. n2+1)) THEN |
---|
1409 | !!$ !i=n2+1 |
---|
1410 | !!$ amhalf=diffup(ipoint) |
---|
1411 | !!$ aphalf=diffair |
---|
1412 | !!$ p=0. |
---|
1413 | !!$ q=0. |
---|
1414 | !!$ goxid(ipoint)=goxid(ipoint)+q |
---|
1415 | !!$ aold1=rkh*amhalf |
---|
1416 | !!$ aold2=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1417 | !!$ aold3=rkh*aphalf |
---|
1418 | !!$ aold4=2.*rk*q |
---|
1419 | !!$ b(ipoint,i)=aold1*uo(ipoint,i-1)+aold2*uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4 |
---|
1420 | !!$ |
---|
1421 | !!$ ELSEIF ((nw(ipoint) .EQ. ns-1) .AND. (i .GE. n2+2) .AND. (i .LE. n3-1)) THEN |
---|
1422 | !!$ |
---|
1423 | !!$ !DO i=n2+2,n3-1 |
---|
1424 | !!$ amhalf=diffair |
---|
1425 | !!$ aphalf=diffair |
---|
1426 | !!$ p=0. |
---|
1427 | !!$ q=0. |
---|
1428 | !!$ aold1=rkh*amhalf |
---|
1429 | !!$ aold2=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1430 | !!$ aold3=rkh*aphalf |
---|
1431 | !!$ aold4=2.*rk*q |
---|
1432 | !!$ b(ipoint,i)=aold1*uo(ipoint,i-1)+aold2*uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4 |
---|
1433 | !!$ !ENDDO |
---|
1434 | !!$ |
---|
1435 | !!$ ENDIF |
---|
1436 | !!$ ENDDO |
---|
1437 | !!$ ENDDO |
---|
1438 | |
---|
1439 | !4eme CAS: (nw .GE. ns): non traite! ma water table ne va jamais au-dessus de la |
---|
1440 | !surface en ce qui me concerne. de toute facon, jamais d oxydation dans Walter dans |
---|
1441 | !la colonne d eau au-dessus de la surface: cala modifie juste le transport |
---|
1442 | |
---|
1443 | |
---|
1444 | |
---|
1445 | !DO ipoint =1,npts |
---|
1446 | ! IF (iday .EQ. 2) THEN |
---|
1447 | ! WRITE(*,*) 'matrice b apres inner point apres cas',b(ipoint,:) |
---|
1448 | ! ENDIF |
---|
1449 | !ENDDO |
---|
1450 | |
---|
1451 | ! right boundary |
---|
1452 | DO ipoint = 1,npts |
---|
1453 | amhalf=diffair |
---|
1454 | aphalf=diffair |
---|
1455 | p=0. |
---|
1456 | q=0. |
---|
1457 | aoldm2=rkh*amhalf |
---|
1458 | aoldm1=2.-rkh*(amhalf+aphalf-h**2*p) |
---|
1459 | aoldn=2.*(rkh*aphalf*catm+rk*q) |
---|
1460 | b(ipoint,n3)=aoldm2*uo(ipoint,n3-1)+aoldm1*uo(ipoint,n3)+aoldn |
---|
1461 | ENDDO |
---|
1462 | |
---|
1463 | |
---|
1464 | |
---|
1465 | |
---|
1466 | !********************************************************* |
---|
1467 | ! solution of tridiagonal system using subroutine tridag |
---|
1468 | ! calculate NEW concentration profiles uo(i) |
---|
1469 | |
---|
1470 | |
---|
1471 | !DO i=n0,n3 |
---|
1472 | DO i=1,n3 |
---|
1473 | DO ipoint =1,npts |
---|
1474 | IF (i .GE. n0(ipoint)) THEN |
---|
1475 | tria(ipoint,i-n0(ipoint)+1)=0. |
---|
1476 | trib(ipoint,i-n0(ipoint)+1)=0. |
---|
1477 | tric(ipoint,i-n0(ipoint)+1)=0. |
---|
1478 | trir(ipoint,i-n0(ipoint)+1)=0. |
---|
1479 | triu(ipoint,i-n0(ipoint)+1)=0. |
---|
1480 | ENDIF |
---|
1481 | ENDDO |
---|
1482 | enddo |
---|
1483 | |
---|
1484 | |
---|
1485 | !do i=n0,n3 |
---|
1486 | DO i=1,n3 |
---|
1487 | DO ipoint =1,npts |
---|
1488 | IF (i .GE. n0(ipoint)) THEN |
---|
1489 | trib(ipoint,i-n0(ipoint)+1)=a(ipoint,i,i) |
---|
1490 | ENDIF |
---|
1491 | ENDDO |
---|
1492 | ENDDO |
---|
1493 | |
---|
1494 | !do i=n0+1,n3 |
---|
1495 | DO i=1,n3 |
---|
1496 | DO ipoint =1,npts |
---|
1497 | IF (i .GE. n0(ipoint)+1) THEN |
---|
1498 | tria(ipoint,i-n0(ipoint)+1)=a(ipoint,i,i-1) |
---|
1499 | ENDIF |
---|
1500 | ENDDO |
---|
1501 | ENDDO |
---|
1502 | |
---|
1503 | !do i=n0,n3-1 |
---|
1504 | DO i=1,n3-1 |
---|
1505 | DO ipoint =1,npts |
---|
1506 | IF (i .GE. n0(ipoint)) THEN |
---|
1507 | tric(ipoint,i-n0(ipoint)+1)=a(ipoint,i,i+1) |
---|
1508 | ENDIF |
---|
1509 | ENDDO |
---|
1510 | ENDDO |
---|
1511 | |
---|
1512 | ! do i=n0,n3 |
---|
1513 | DO i=1,n3 |
---|
1514 | DO ipoint =1,npts |
---|
1515 | IF (i .GE. n0(ipoint)) THEN |
---|
1516 | trir(ipoint,i-n0(ipoint)+1)=b(ipoint,i) |
---|
1517 | ENDIF |
---|
1518 | ENDDO |
---|
1519 | ENDDO |
---|
1520 | |
---|
1521 | DO ipoint =1,npts |
---|
1522 | ndim(ipoint)=n3-n0(ipoint)+1 |
---|
1523 | ENDDO |
---|
1524 | |
---|
1525 | !WRITE(*,*) ' ' |
---|
1526 | !CALL tridiag(tria,trib,tric,trir,triu,ndim) |
---|
1527 | !REMPLACEMENT DE LA ROUTINE TRIDAG |
---|
1528 | !WRITE(*,*) ' ' |
---|
1529 | |
---|
1530 | DO ipoint =1,npts |
---|
1531 | bet(ipoint)=trib(ipoint,1) |
---|
1532 | triu(ipoint,1)=trir(ipoint,1)/bet(ipoint) |
---|
1533 | ENDDO |
---|
1534 | |
---|
1535 | !DO j=2,ndim |
---|
1536 | !enleve j=2,500 et met a la place j=2,n3 |
---|
1537 | DO j=2,n3 |
---|
1538 | DO ipoint =1,npts |
---|
1539 | IF (j .LE. ndim(ipoint)) THEN |
---|
1540 | gam(ipoint,j)=tric(ipoint,j-1)/bet(ipoint) |
---|
1541 | ! bet(ipoint)=trib(ipoint,j)-tria(ipoint,j)*gam(ipoint,j) |
---|
1542 | IF ((trib(ipoint,j)-tria(ipoint,j)*gam(ipoint,j)) .EQ. 0) THEN |
---|
1543 | bet(ipoint)=trib(ipoint,j)-tria(ipoint,j)*gam(ipoint,j) + 10 |
---|
1544 | ELSE |
---|
1545 | bet(ipoint)=trib(ipoint,j)-tria(ipoint,j)*gam(ipoint,j) |
---|
1546 | ENDIF |
---|
1547 | triu(ipoint,j)=(trir(ipoint,j)-tria(ipoint,j)*triu(ipoint,j-1))/bet(ipoint) |
---|
1548 | ENDIF |
---|
1549 | ENDDO |
---|
1550 | ENDDO |
---|
1551 | |
---|
1552 | !DO j=ndim-1,1,-1 |
---|
1553 | DO j=n3-1,1,-1 |
---|
1554 | DO ipoint =1,npts |
---|
1555 | IF (j .LE. ndim(ipoint)-1) THEN |
---|
1556 | triu(ipoint,j)=triu(ipoint,j)-gam(ipoint,j+1)*triu(ipoint,j+1) |
---|
1557 | ENDIF |
---|
1558 | ENDDO |
---|
1559 | ENDDO |
---|
1560 | !FIN DE REMPLACEMENT DE LA ROUTINE TRIDAG |
---|
1561 | |
---|
1562 | |
---|
1563 | !DO i=n0,n3 |
---|
1564 | DO i=1,n3 |
---|
1565 | DO ipoint =1,npts |
---|
1566 | IF (i .GE. n0(ipoint)) THEN |
---|
1567 | uo(ipoint,i)=triu(ipoint,i-n0(ipoint)+1) |
---|
1568 | ! if (due to small number of time steps per day and strong variation in |
---|
1569 | ! water table (i.e. large change in vertical profile of diffusion coefficient) |
---|
1570 | ! negative methane concentrations occur: set to zero |
---|
1571 | ! should usually not happen - can be avoided by increasing nday (from 24 to |
---|
1572 | ! 240 or 2400 or ...) or possibly using a different numerical scheme |
---|
1573 | IF (uo(ipoint,i) .LT. 0.) THEN |
---|
1574 | negconc(ipoint)=negconc(ipoint)-(uo(ipoint,i)/fpart(ipoint,i)) |
---|
1575 | ENDIF |
---|
1576 | uo(ipoint,i)=MAX(0.,triu(ipoint,i-n0(ipoint)+1)) |
---|
1577 | ! transform from partial pressure into concentration unit |
---|
1578 | uo(ipoint,i)=uo(ipoint,i)/fpart(ipoint,i) |
---|
1579 | ENDIF |
---|
1580 | ENDDO |
---|
1581 | ENDDO |
---|
1582 | |
---|
1583 | |
---|
1584 | |
---|
1585 | !la nouvelle valeur de uo va etre utilisée comme valeur initiale au pas de temps suivant.... |
---|
1586 | !************************************ |
---|
1587 | !************************************ |
---|
1588 | |
---|
1589 | |
---|
1590 | !************************************************************** |
---|
1591 | |
---|
1592 | ! plant-mediated transport: |
---|
1593 | !************************** |
---|
1594 | |
---|
1595 | DO ipoint =1,npts |
---|
1596 | tplant(ipoint)=0. |
---|
1597 | ENDDO |
---|
1598 | |
---|
1599 | !DO i=n0,n3 |
---|
1600 | DO i=1,n3 |
---|
1601 | DO ipoint =1,npts |
---|
1602 | IF (i .GE. n0(ipoint)) THEN |
---|
1603 | ! calculated fraction of methane that enters plants: |
---|
1604 | ! tveg: vegetation type from Sread.f |
---|
1605 | ! dveg: scaling factor from para.h |
---|
1606 | ! fgrow: growing stage of plants as calculated above (LAI) |
---|
1607 | ! root(i): vertical root distribution as calculated above |
---|
1608 | cplant=MIN(tveg(ipoint)*dveg*fgrow(ipoint)*root(ipoint,i)*10.*rk,1.) |
---|
1609 | ! amount (concentration) of methane that enters plants |
---|
1610 | cplant=cplant*uo(ipoint,i) |
---|
1611 | dplant(ipoint)=dplant(ipoint)+cplant |
---|
1612 | ! amount (concentration) of methane that remains in soil |
---|
1613 | uo(ipoint,i)=uo(ipoint,i)-cplant |
---|
1614 | ! only fraction (1.-pox) of methane entering plants is emitted |
---|
1615 | ! to the atmosphere - the rest (pox) is oxidized in the rhizosphere |
---|
1616 | tplant(ipoint)=tplant(ipoint)+cplant*(1.-pox) |
---|
1617 | ENDIF |
---|
1618 | ENDDO |
---|
1619 | ENDDO |
---|
1620 | ! transform amount (concentration) into flux unit (mg/m^2*d) |
---|
1621 | flplant(:,iday)=tplant(:)*funit |
---|
1622 | |
---|
1623 | !DO ipoint =1,npts |
---|
1624 | ! IF (iday .EQ. 1) THEN |
---|
1625 | ! WRITE(*,*) 'tveg',tveg,'fgrow',fgrow |
---|
1626 | ! WRITE(*,*) 'root',root |
---|
1627 | ! ENDIF |
---|
1628 | !ENDDO |
---|
1629 | |
---|
1630 | !DO ipoint =1,npts |
---|
1631 | ! IF (iday .EQ. 1) THEN |
---|
1632 | ! WRITE(*,*) 'matrice uo apres plante',uo |
---|
1633 | ! ENDIF |
---|
1634 | ! ENDDO |
---|
1635 | |
---|
1636 | |
---|
1637 | ! bubble transport: |
---|
1638 | !****************** |
---|
1639 | |
---|
1640 | DO ipoint =1,npts |
---|
1641 | tmax(ipoint)=0. |
---|
1642 | ENDDO |
---|
1643 | |
---|
1644 | !DO i=n0,nw |
---|
1645 | DO i=1,n2 |
---|
1646 | DO ipoint =1,npts |
---|
1647 | IF ((i .GE. n0(ipoint)) .AND. (i .LE. nw(ipoint))) THEN |
---|
1648 | ! methane in water: if concentration above threshold concentration rcmax |
---|
1649 | ! (calculated in Scalc.f): bubble formation |
---|
1650 | cdiff=uo(ipoint,i)-rcmax(ipoint) |
---|
1651 | cbub=MAX(0.,cdiff) |
---|
1652 | tmax(ipoint)=tmax(ipoint)+cbub |
---|
1653 | uo(ipoint,i)=MIN(uo(ipoint,i),rcmax(ipoint)) |
---|
1654 | ENDIF |
---|
1655 | ENDDO |
---|
1656 | ENDDO |
---|
1657 | |
---|
1658 | |
---|
1659 | !test Estelle |
---|
1660 | ! IF (tmax .GT. 0.) then |
---|
1661 | ! WRITE(*,*) 'check', tmax, n0, nw, uo(nw), cdiff |
---|
1662 | ! ENDIF |
---|
1663 | ! if water table is above soil surface: bubble flux is emitted directly |
---|
1664 | ! to atmosphere (mg/m^2*d) |
---|
1665 | |
---|
1666 | DO ipoint =1,npts |
---|
1667 | IF (nw(ipoint) .GE. ns) THEN |
---|
1668 | flbub(ipoint,iday)=tmax(ipoint)*funit |
---|
1669 | flbubdiff(ipoint,iday)=0. |
---|
1670 | tmax(ipoint)=0. |
---|
1671 | ! if water table is below soil surface: |
---|
1672 | ELSE |
---|
1673 | flbub(ipoint,iday)=0. |
---|
1674 | flbubdiff(ipoint,iday)=0. |
---|
1675 | ! 1. if water table is 5 or less cm below soil surface: flbubdiff is |
---|
1676 | ! calculated in the same way as flbub and added to the total flux (mg/m^2*d) |
---|
1677 | IF (nw(ipoint) .GE. ns-5) THEN |
---|
1678 | flbubdiff(ipoint,iday)=tmax(ipoint)*funit |
---|
1679 | tmax(ipoint)=0. |
---|
1680 | ENDIF |
---|
1681 | ENDIF |
---|
1682 | ENDDO |
---|
1683 | |
---|
1684 | !test Estelle |
---|
1685 | ! IF (flbub(iday) .gt. 0.) then |
---|
1686 | ! WRITE(*,*) 'buble', iday,flbub(iday),tmax, cdiff, cbub, nw |
---|
1687 | ! ENDIF |
---|
1688 | |
---|
1689 | |
---|
1690 | ! 2. if water table is more than 5 cm below soil surface: no bubble flux |
---|
1691 | ! is directly emitted to the atmosphere, BUT methane amount from bubble is |
---|
1692 | ! added to the first layer above the soil surface |
---|
1693 | ! (the distinction more/less than 5 cm below soil has been made, because |
---|
1694 | ! adding tmax to uo(nw+1) if the water table is less than 5cm below the soil |
---|
1695 | ! surface 'disturbed' the system too much and more time steps per day (larger |
---|
1696 | ! nday) is needed) (this solution avoids this) |
---|
1697 | DO ipoint =1,npts |
---|
1698 | uo(ipoint,nw(ipoint)+1)=uo(ipoint,nw(ipoint)+1)+tmax(ipoint) |
---|
1699 | ENDDO |
---|
1700 | |
---|
1701 | !DO ipoint =1,npts |
---|
1702 | ! IF (iday .EQ. 1) THEN |
---|
1703 | ! WRITE(*,*) 'matrice uo apres bubble',uo |
---|
1704 | ! ENDIF |
---|
1705 | !ENDDO |
---|
1706 | |
---|
1707 | ! calculate diffusive flux from the concentration gradient (mg/m^2*d): |
---|
1708 | !********************************************************************** |
---|
1709 | DO ipoint =1,npts |
---|
1710 | fdifftime(ipoint,iday)=((uo(ipoint,n3-2)-uo(ipoint,n3-1))/0.1)*38.4*diffair |
---|
1711 | ENDDO |
---|
1712 | |
---|
1713 | ! transform from concentration into partial pressure unit |
---|
1714 | !do i=n0,n3 |
---|
1715 | DO i=1,n3 |
---|
1716 | DO ipoint =1,npts |
---|
1717 | IF (i .GE. n0(ipoint)) THEN |
---|
1718 | uo(ipoint,i)=uo(ipoint,i)*fpart(ipoint,i) |
---|
1719 | ENDIF |
---|
1720 | ENDDO |
---|
1721 | enddo |
---|
1722 | |
---|
1723 | |
---|
1724 | |
---|
1725 | ! DO ipoint =1,npts |
---|
1726 | ! IF (iday .EQ. 1) THEN |
---|
1727 | ! WRITE(*,*) 'matrice uo FIN TIMELOOP',uo |
---|
1728 | ! ENDIF |
---|
1729 | ! ENDDO |
---|
1730 | |
---|
1731 | ! |
---|
1732 | ! end of (iday=1,nday) loop, i.e. time step per day |
---|
1733 | !*************************************************** |
---|
1734 | |
---|
1735 | ENDDO |
---|
1736 | |
---|
1737 | !*************************************************** |
---|
1738 | |
---|
1739 | |
---|
1740 | ! water-air interface |
---|
1741 | ! transform from partial pressure into concentration unit |
---|
1742 | !do i=n0,n3 |
---|
1743 | DO i=1,n3 |
---|
1744 | DO ipoint =1,npts |
---|
1745 | IF (i .GE. n0(ipoint)) THEN |
---|
1746 | uo(ipoint,i)=uo(ipoint,i)/fpart(ipoint,i) |
---|
1747 | ENDIF |
---|
1748 | ENDDO |
---|
1749 | ENDDO |
---|
1750 | |
---|
1751 | |
---|
1752 | ! calculate conc kept in soil - |
---|
1753 | ! necessary to calculate fluxdiff (diffusive flux from mass balance) |
---|
1754 | DO i=1,nvert |
---|
1755 | DO ipoint =1,npts |
---|
1756 | uout(ipoint,i)=uo(ipoint,i) |
---|
1757 | ENDDO |
---|
1758 | ENDDO |
---|
1759 | |
---|
1760 | !do i=n0,n3 |
---|
1761 | DO i=1,n3 |
---|
1762 | DO ipoint =1,npts |
---|
1763 | IF (i .GE. n0(ipoint)) THEN |
---|
1764 | tout(ipoint)=tout(ipoint)+(uout(ipoint,i)-uold2(ipoint,i)) |
---|
1765 | ENDIF |
---|
1766 | ENDDO |
---|
1767 | ENDDO |
---|
1768 | |
---|
1769 | !DO ipoint =1,npts |
---|
1770 | ! tout(ipoint)=0. |
---|
1771 | !ENDDO |
---|
1772 | |
---|
1773 | tout(:)=0. |
---|
1774 | |
---|
1775 | !do i=n0,n3 |
---|
1776 | DO i=1,n3 |
---|
1777 | DO ipoint =1,npts |
---|
1778 | IF (i .GE. n0(ipoint)) THEN |
---|
1779 | tout(ipoint)=tout(ipoint)+(uout(ipoint,i)-uold2(ipoint,i)) |
---|
1780 | ENDIF |
---|
1781 | ENDDO |
---|
1782 | ENDDO |
---|
1783 | |
---|
1784 | |
---|
1785 | DO i=1,nvert |
---|
1786 | DO ipoint =1,npts |
---|
1787 | uold2(ipoint,i)=uout(ipoint,i) |
---|
1788 | ENDDO |
---|
1789 | ENDDO |
---|
1790 | |
---|
1791 | |
---|
1792 | !***************************** |
---|
1793 | !vire concout |
---|
1794 | |
---|
1795 | ! write CH4conc to varible concout after each time interval (at iday=nday) |
---|
1796 | ! concout(i,itime) is needed to write CH4conc profiles in Soutput.f |
---|
1797 | ! do i=1,n3 |
---|
1798 | ! concout(i,itime)=uout(i) |
---|
1799 | ! enddo |
---|
1800 | ! do i=n3+1,n |
---|
1801 | ! transform from partial pressure into concentration unit |
---|
1802 | ! concout(i,itime)=catm/fpart(i) |
---|
1803 | ! enddo |
---|
1804 | |
---|
1805 | ! calculate daily means of bubble flux, fluxbubdiff, plant flux, and |
---|
1806 | ! diffusive flux (from concentration gradient) (mg/m^2*d) |
---|
1807 | |
---|
1808 | !******************************* |
---|
1809 | !la j enleve la dimension temps des sorties... |
---|
1810 | !******************************* |
---|
1811 | |
---|
1812 | |
---|
1813 | |
---|
1814 | DO ipoint =1,npts |
---|
1815 | fluxbub(ipoint)=0. |
---|
1816 | fluxbubdiff(ipoint)=0. |
---|
1817 | fluxplant(ipoint)=0. |
---|
1818 | fdiffday(ipoint)=0. |
---|
1819 | ENDDO |
---|
1820 | |
---|
1821 | |
---|
1822 | DO i=1,nday |
---|
1823 | DO ipoint =1,npts |
---|
1824 | fluxbub(ipoint)=fluxbub(ipoint)+flbub(ipoint,i) |
---|
1825 | fluxbubdiff(ipoint)=fluxbubdiff(ipoint)+flbubdiff(ipoint,i) |
---|
1826 | fluxplant(ipoint)=fluxplant(ipoint)+flplant(ipoint,i) |
---|
1827 | fdiffday(ipoint)=fdiffday(ipoint)+fdifftime(ipoint,i) |
---|
1828 | ENDDO |
---|
1829 | ENDDO |
---|
1830 | |
---|
1831 | fluxbub(:)=fluxbub(:)/nday |
---|
1832 | fluxbubdiff(:)=fluxbubdiff(:)/nday |
---|
1833 | fluxplant(:)=fluxplant(:)/nday |
---|
1834 | fdiffday(:)=fdiffday(:)/nday |
---|
1835 | |
---|
1836 | |
---|
1837 | ! |
---|
1838 | ! calculate diffusive flux ( from mass balance): |
---|
1839 | !************************************************ |
---|
1840 | ! use amounts in concentration unit and tranform in flux unit (mg/m^2*d) |
---|
1841 | ! wpro: amount of methane produced during entire day |
---|
1842 | ! goxid: amount of methane oxidized during entire day |
---|
1843 | ! mutiply by rk (depending on number of time steps per day) |
---|
1844 | ! dplant: amount of methane that entered plants during entire day |
---|
1845 | ! (includes both: methane emitted through plants AND oxidized in rhizosphere) |
---|
1846 | ! tout: additional methane stored in soil during entire day |
---|
1847 | ! wtdiff: see above |
---|
1848 | ! fluxbub: bubble flux (if water table > soil surface) |
---|
1849 | ! bubble flux (if water table less than 5 cm below soil surface (see above)) |
---|
1850 | ! negconc: see above |
---|
1851 | fluxdiff(:)= ((rk*(wpro(:)+goxid(:)) - dplant(:)-tout(:)-wtdiff)/nday)*funit & |
---|
1852 | & - fluxbub(:) - fluxbubdiff(:) + negconc(:)*(funit/nday) |
---|
1853 | |
---|
1854 | ! fluxbubdiff is added to diffusive flux (mg/m^2*d) |
---|
1855 | |
---|
1856 | |
---|
1857 | fluxdiff(:) = fluxdiff(:)+fluxbubdiff(:) |
---|
1858 | fdiffday(:) = fdiffday(:)+fluxbubdiff(:) |
---|
1859 | |
---|
1860 | |
---|
1861 | !enleve boucle if sur ihump |
---|
1862 | |
---|
1863 | ! ELSE |
---|
1864 | ! fluxdiff=0. |
---|
1865 | ! fdiffday=0. |
---|
1866 | ! fluxbub=0. |
---|
1867 | ! fluxplant=0. |
---|
1868 | ! fluxtot=0. |
---|
1869 | ! !do i=1,n |
---|
1870 | ! ! concout(i,itime)=0. |
---|
1871 | ! !enddo |
---|
1872 | ! |
---|
1873 | !! end of (if (ijump .eq. 1)) 'loop' |
---|
1874 | !!*********************************** |
---|
1875 | ! ENDIF |
---|
1876 | |
---|
1877 | |
---|
1878 | |
---|
1879 | |
---|
1880 | |
---|
1881 | ! to calculate wtdiff (see above) |
---|
1882 | ! define uwp1 |
---|
1883 | |
---|
1884 | !non necessaire: la WTD ne se modifie pas d un pas de tps a l autre |
---|
1885 | |
---|
1886 | ! uwp0=uo(ipoint,nw) |
---|
1887 | ! uwp1=uo(ipoint,nw+1) |
---|
1888 | ! uwp2=uo(ipoint,nw+2) |
---|
1889 | |
---|
1890 | |
---|
1891 | |
---|
1892 | ! calculate total flux fluxtot (using diffusive flux from mass balance) |
---|
1893 | ! and sum of diffusive flux from mass balance and bubble flux (for plots) |
---|
1894 | ! (mg/m^2*d) |
---|
1895 | ! do i=1,ntime |
---|
1896 | |
---|
1897 | fluxtot(:)=fluxdiff(:)+fluxbub(:)+fluxplant(:) |
---|
1898 | fbubdiff(:)=fluxdiff(:)+fluxbub(:) |
---|
1899 | |
---|
1900 | |
---|
1901 | DO ipoint =1,npts |
---|
1902 | ch4_flux_density_dif(ipoint) = fluxdiff(ipoint) |
---|
1903 | ch4_flux_density_bub(ipoint) = fluxbub(ipoint) |
---|
1904 | ch4_flux_density_pla(ipoint) = fluxplant(ipoint) |
---|
1905 | ENDDO |
---|
1906 | |
---|
1907 | |
---|
1908 | WHERE (fluxtot(:) .LE. 0.0) |
---|
1909 | ch4_flux_density_tot(:) = 0.0 |
---|
1910 | ELSEWHERE |
---|
1911 | ch4_flux_density_tot(:) = fluxtot(:) |
---|
1912 | endwhere |
---|
1913 | |
---|
1914 | !pour contrebalancer valeur par defaut mise pour les grid-cells ou veget_max=0 |
---|
1915 | WHERE ((veget_max(:,1) .LE. 0.95) .AND. (sum_veget_nat_ss_nu(:) .GT. 0.1)) |
---|
1916 | ch4_flux_density_tot(:) = ch4_flux_density_tot(:) |
---|
1917 | ELSEWHERE |
---|
1918 | ch4_flux_density_tot(:) = 0.0 |
---|
1919 | ENDWHERE |
---|
1920 | |
---|
1921 | |
---|
1922 | |
---|
1923 | END SUBROUTINE ch4_wet_flux_density_wet2 |
---|
1924 | !******************************************************************** |
---|
1925 | !!$ |
---|
1926 | !!$SUBROUTINE tridag(a,b,c,r,u,n) |
---|
1927 | !!$! |
---|
1928 | !!$! Press et al., Numerical Recipes, FORTRAN, p.43 |
---|
1929 | !!$!************************************************ |
---|
1930 | !!$ |
---|
1931 | !!$ ! Domain size |
---|
1932 | !!$ INTEGER(i_std), INTENT(in) :: n |
---|
1933 | !!$ REAL(r_std),DIMENSION(n), INTENT(in) :: a |
---|
1934 | !!$ REAL(r_std),DIMENSION(n), INTENT(in) :: b |
---|
1935 | !!$ REAL(r_std),DIMENSION(n), INTENT(in) :: c |
---|
1936 | !!$ REAL(r_std),DIMENSION(n), INTENT(in) :: r |
---|
1937 | !!$ REAL(r_std),DIMENSION(n), INTENT(out) :: u |
---|
1938 | !!$ INTEGER(stnd),PARAMETER :: nmax=500 |
---|
1939 | !!$ INTEGER(i_std) :: j |
---|
1940 | !!$ REAL(r_std),DIMENSION(nmax) :: gam |
---|
1941 | !!$ REAL(r_std) :: bet |
---|
1942 | !!$ |
---|
1943 | !!$ bet=b(1) |
---|
1944 | !!$ |
---|
1945 | !!$ u(1)=r(1)/bet |
---|
1946 | !!$ |
---|
1947 | !!$ do j=2,n |
---|
1948 | !!$ gam(j)=c(j-1)/bet |
---|
1949 | !!$ bet=b(j)-a(j)*gam(j) |
---|
1950 | !!$ u(j)=(r(j)-a(j)*u(j-1))/bet |
---|
1951 | !!$ enddo |
---|
1952 | !!$ |
---|
1953 | !!$ do j=n-1,1,-1 |
---|
1954 | !!$ u(j)=u(j)-gam(j+1)*u(j+1) |
---|
1955 | !!$ enddo |
---|
1956 | !!$ |
---|
1957 | !!$END SUBROUTINE tridag |
---|
1958 | |
---|
1959 | |
---|
1960 | END MODULE stomate_wet_ch4_pt_ter_wet2 |
---|