Context Navigation

stomate_wet_ch4_pt_ter_0.f90 @ 6940

Last change on this file since 6940 was 6940, checked in by jinfeng.chang, 4 years ago
add missing files for ORCHIDEE-GMv3.2
File size: 66.4 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE-GMv3.2/ORCHIDEE/src_stomate/stomate_wet_ch4_pt_ter_0.f90 @ 6940

Download in other formats: