Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

stomate_wet_ch4_pt_ter_wet2.f90 @ 5491

Last change on this file since 5491 was 5080, checked in by chunjing.qiu, 6 years ago
soil freezing, soil moisture, fwet bugs fixed
File size: 66.9 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
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(:)=1641/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=diffair0.66rpv
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-h2p)
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-h2p)
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-h2p)
849	! a(i,j3)=-rkh*aphalf
850	! enddo
851
852
853	diffgrdo=0.570.66rpv
854	diffgrup=0.02480.66rpv
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-h2p)
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-h2p)
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-h2p)
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-h2p)
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-h2p)
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-h2p)
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-h2p)
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-h2p)
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-h2p)
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-h2p)
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-h2p)
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-h2p)
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-h2p)
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)(rq10rexp)carbon(itime)*(0.001239567)
1057	! source=forg(n0)(rq10rexp)carbon(itime)r0pl(0.0012)
1058	! source=5.0(rq10*rexp)
1059	source=forg(ipoint,n0(ipoint))(rq10rexp)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-h2p)
1068	aold2=rkh*(aphalf+amhalf)
1069	aold3=2.rkq
1070	b(ipoint,n0(ipoint))=aold1uo(ipoint,n0(ipoint))+aold2uo(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)(rq10rexp)carbon(itime)*(0.001239567)
1098	!modifie le calcul de la production --> utilise le carbone comme proxi du substrat
1099	source=forg(ipoint,i)(rq10rexp)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-h2p)
1109	aold3=rkh*aphalf
1110	aold4=2.rkq
1111	b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(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.570.66rpv
1123	diffgrup=0.02480.66rpv
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)(rq10rexp)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-h2p)
1158	aold3=rkh*aphalf
1159	aold4=2.rkq
1160	b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(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-h2p)
1190	aold3=rkh*aphalf
1191	aold4=2.rkq
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)=aold1uo(ipoint,i-1)+aold2uo(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-h2p)
1224	aold3=rkh*aphalf
1225	aold4=2.rkq
1226	b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(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-h2p)
1255	aold3=rkh*aphalf
1256	aold4=2.rkq
1257	b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(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-h2p)
1267	aold3=rkh*aphalf
1268	aold4=2.rkq
1269	b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(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)(rq10rexp)carbon(itime)*(0.001239567)
1296	!!$! source=5.0(rq10*rexp)
1297	!!$ source=forg(ipoint,i)(rq10rexp)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-h2p)
1307	!!$ aold3=rkh*aphalf
1308	!!$ aold4=2.rkq
1309	!!$ b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(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-h2p)
1319	!!$ aold3=rkh*aphalf
1320	!!$ aold4=2.rkq
1321	!!$ b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(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-h2p)
1331	!!$ aold3=rkh*aphalf
1332	!!$ aold4=2.rkq
1333	!!$ b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(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)(rq10rexp)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-h2p)
1372	!!$ aold3=rkh*aphalf
1373	!!$ aold4=2.rkq
1374	!!$ b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(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-h2p)
1404	!!$ aold3=rkh*aphalf
1405	!!$ aold4=2.rkq
1406	!!$ b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(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-h2p)
1417	!!$ aold3=rkh*aphalf
1418	!!$ aold4=2.rkq
1419	!!$ b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(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-h2p)
1430	!!$ aold3=rkh*aphalf
1431	!!$ aold4=2.rkq
1432	!!$ b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(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-h2p)
1459	aoldn=2.(rkhaphalfcatm+rkq)
1460	b(ipoint,n3)=aoldm2uo(ipoint,n3-1)+aoldm1uo(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)dvegfgrow(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.4diffair
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

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

Context Navigation

source: branches/publications/ORCHIDEE-PEAT_r5488/src_stomate/stomate_wet_ch4_pt_ter_wet2.f90 @ 5491

Download in other formats: