Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

stomate_wet_ch4_pt_ter_wet3.f90 @ 7346

Last change on this file since 7346 was 5080, checked in by chunjing.qiu, 6 years ago
soil freezing, soil moisture, fwet bugs fixed
File size: 66.5 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_wt3 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_wt3 cm)
10
11	MODULE stomate_wet_ch4_pt_ter_wet3
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_wet3,ch4_wet_flux_density_clear_wet3
28
29	! first call
30	LOGICAL, SAVE :: firstcall = .TRUE.
31
32	CONTAINS
33
34	SUBROUTINE ch4_wet_flux_density_clear_wet3
35	firstcall=.TRUE.
36	END SUBROUTINE ch4_wet_flux_density_clear_wet3
37
38
39	SUBROUTINE ch4_wet_flux_density_wet3 ( 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=-15
369	rwater(:) = pwater_wet3
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	b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4
1193
1194
1195	ELSEIF ((i .GE. n1(ipoint)+2) .AND. (i .LE. n2-1)) THEN
1196	!do i=n1+2,n2-1
1197	amhalf=diffup(ipoint)
1198	aphalf=diffup(ipoint)
1199	p=0.
1200	! oxidation following Michaelis-Menten kinetics
1201	IF (uo(ipoint,i) .GT. 0.) THEN
1202	! transform from partial pressure into concentration unit
1203	uo(ipoint,i)=uo(ipoint,i)/fpart(ipoint,i)
1204	! Michaelis Menten Equation
1205	q=-rvmax(ipoint)*uo(ipoint,i)/(uo(ipoint,i)+rkm)
1206	! temperature dependence (oxq10 from para.h)
1207	q=(q/oxq10)(oxq10*((sou(ipoint,i)-MAX(0.,tmean(ipoint)))/10.))
1208	goxid(ipoint)=goxid(ipoint)+q
1209	! transform from concentration into partial pressure unit
1210	uo(ipoint,i)=uo(ipoint,i)*fpart(ipoint,i)
1211	! in case of negative concentrations: no oxidation
1212	! (can occur if number of time steps per day (nday) is small and water
1213	! table varies strongly) (help: increase nday (or change numerical scheme!))
1214	ELSE
1215	q=0.
1216	goxid(ipoint)=goxid(ipoint)+q
1217	ENDIF
1218	! transform from concentration into partial pressure unit
1219	q=q*fpart(ipoint,i)
1220	aold1=rkh*amhalf
1221	aold2=2.-rkh(amhalf+aphalf-h2p)
1222	aold3=rkh*aphalf
1223	aold4=2.rkq
1224	b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4
1225
1226	ELSEIF (i .EQ. n2) THEN
1227	!i=n2
1228	amhalf=diffup(ipoint)
1229	aphalf=diffair
1230	p=0.
1231	! oxidation following Michaelis-Menten kinetics
1232	IF (uo(ipoint,i) .GT. 0.) THEN
1233	! transform from partial pressure into concentration unit
1234	uo(ipoint,i)=uo(ipoint,i)/fpart(ipoint,i)
1235	! Michaelis Menten Equation
1236	q=-rvmax(ipoint)*uo(ipoint,i)/(uo(ipoint,i)+rkm)
1237	! temperature dependence (oxq10 from para.h)
1238	q=(q/oxq10)(oxq10*((sou(ipoint,i)-MAX(0.,tmean(ipoint)))/10.))
1239	goxid(ipoint)=goxid(ipoint)+q
1240	! transform from concentration into partial pressure unit
1241	uo(ipoint,i)=uo(ipoint,i)*fpart(ipoint,i)
1242	! in case of negative concentrations: no oxidation
1243	! (can occur if number of time steps per day (nday) is small and water
1244	! table varies strongly) (help: increase nday (or change numerical scheme!))
1245	ELSE
1246	q=0.
1247	goxid(ipoint)=goxid(ipoint)+q
1248	ENDIF
1249	! transform from concentration into partial pressure unit
1250	q=q*fpart(ipoint,i)
1251	aold1=rkh*amhalf
1252	aold2=2.-rkh(amhalf+aphalf-h2p)
1253	aold3=rkh*aphalf
1254	aold4=2.rkq
1255	b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4
1256
1257	ELSEIF ((i .GE. n2+1) .AND. (i .LE. n3-1)) THEN
1258	!DO i=n2+1,n3-1
1259	amhalf=diffair
1260	aphalf=diffair
1261	p=0.
1262	q=0.
1263	aold1=rkh*amhalf
1264	aold2=2.-rkh(amhalf+aphalf-h2p)
1265	aold3=rkh*aphalf
1266	aold4=2.rkq
1267	b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4
1268
1269	ENDIF
1270	ENDDO
1271	ENDDO
1272
1273	!!$
1274	!!$
1275	!!$ DO i = 1,n3-1 !sur-calcul: au depart j avais mis i = n0,n3-1
1276	!!$ DO ipoint = 1,npts
1277	!!$
1278	!!$ !2nd CAS: (nw(ipoint) .EQ. ns)
1279	!!$ !diffusion coefficients water-air INTERFACE (soil water --> soil air)
1280	!!$ IF ((nw(ipoint) .EQ. ns) .AND. (i .EQ. n1(ipoint))) THEN
1281	!!$ !i=n1
1282	!!$ amhalf=diffdo(ipoint)
1283	!!$ aphalf=diffgrdo
1284	!!$ p=0.
1285	!!$! q10 dependent production rate (no production if T<0oC, i.e.sou(i)<=0oC)
1286	!!$! forg(i): vertical distribution of substrate in soil (calculated in Scalc.f)
1287	!!$! fin(itime): describes seasonality of NPP
1288	!!$! rq10: Q10 value of production rate (from para.h)
1289	!!$! r0pl: tuning parameter sr0pl (scaled in Scalc.f)
1290	!!$ IF (sou(ipoint,i) .GT. 0.) THEN
1291	!!$ rexp=(sou(ipoint,i)-MAX(0.,tmean(ipoint)))/10.
1292	!!$! source=forg(i)r0pl(rq10*rexp)fin(itime)
1293	!!$! source=forg(i)(rq10rexp)carbon(itime)*(0.001239567)
1294	!!$! source=5.0(rq10*rexp)
1295	!!$ source=forg(ipoint,i)(rq10rexp)substrat(ipoint)*alpha_CH4
1296	!!$ ELSE
1297	!!$ source=0.
1298	!!$ ENDIF
1299	!!$ q=source
1300	!!$ wpro(ipoint)=wpro(ipoint)+q
1301	!!$! transform from concentration into partial pressure unit
1302	!!$ q=q*fpart(ipoint,i)
1303	!!$ aold1=rkh*amhalf
1304	!!$ aold2=2.-rkh(amhalf+aphalf-h2p)
1305	!!$ aold3=rkh*aphalf
1306	!!$ aold4=2.rkq
1307	!!$ b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4
1308	!!$
1309	!!$ ELSEIF ((nw(ipoint) .EQ. ns) .AND. (i .EQ. n1(ipoint)+1)) THEN
1310	!!$ !i=n1+1
1311	!!$ amhalf=diffgrup
1312	!!$ aphalf=diffair
1313	!!$ p=0.
1314	!!$ q=0.
1315	!!$ aold1=rkh*amhalf
1316	!!$ aold2=2.-rkh(amhalf+aphalf-h2p)
1317	!!$ aold3=rkh*aphalf
1318	!!$ aold4=2.rkq
1319	!!$ b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4
1320	!!$
1321	!!$ ELSEIF ((nw(ipoint) .EQ. ns) .AND. (i .GE. n1(ipoint)+2) .AND. (i .LE. n3-1)) THEN
1322	!!$ !DO i=n1+2,n3-1
1323	!!$ amhalf=diffair
1324	!!$ aphalf=diffair
1325	!!$ p=0.
1326	!!$ q=0.
1327	!!$ aold1=rkh*amhalf
1328	!!$ aold2=2.-rkh(amhalf+aphalf-h2p)
1329	!!$ aold3=rkh*aphalf
1330	!!$ aold4=2.rkq
1331	!!$ b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4
1332	!!$ !ENDDO
1333	!!$
1334	!!$ ENDIF
1335	!!$ ENDDO
1336	!!$ ENDDO
1337	!!$
1338	!!$
1339
1340
1341
1342	!!$ DO i = 1,n3-1 !sur-calcul: au depart j avais mis i = n0,n3-1
1343	!!$ DO ipoint = 1,npts
1344	!!$
1345	!!$ !3eme CAS: (nw(ipoint) .EQ. ns-1)
1346	!!$ !diffusion coefficients water-air INTERFACE (soil water --> soil air)
1347	!!$ IF ((nw(ipoint) .EQ. ns-1) .AND. (i .EQ. n1(ipoint))) THEN
1348	!!$ !i=n1
1349	!!$ amhalf=diffdo(ipoint)
1350	!!$ aphalf=diffgrdo
1351	!!$ p=0.
1352	!!$! q10 dependent production rate (no production if T<0oC, i.e.sou(i)<=0oC)
1353	!!$! forg(i): vertical distribution of substrate in soil (calculated in Scalc.f)
1354	!!$! fin(itime): describes seasonality of NPP
1355	!!$! rq10: Q10 value of production rate (from para.h)
1356	!!$! r0pl: tuning parameter sr0pl (scaled in Scalc.f)
1357	!!$ IF (sou(ipoint,i) .GT. 0.) THEN
1358	!!$ rexp=(sou(ipoint,i)-MAX(0.,tmean(ipoint)))/10.
1359	!!$ !source=forg(i)r0pl(rq10*rexp)fin(itime)
1360	!!$ source=forg(ipoint,i)(rq10rexp)substrat(ipoint)*alpha_CH4
1361	!!$ ELSE
1362	!!$ source=0.
1363	!!$ ENDIF
1364	!!$ q=source
1365	!!$ wpro(ipoint)=wpro(ipoint)+q
1366	!!$! transform from concentration into partial pressure unit
1367	!!$ q=q*fpart(ipoint,i)
1368	!!$ aold1=rkh*amhalf
1369	!!$ aold2=2.-rkh(amhalf+aphalf-h2p)
1370	!!$ aold3=rkh*aphalf
1371	!!$ aold4=2.rkq
1372	!!$ b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4
1373	!!$
1374	!!$
1375	!!$ ELSEIF ((nw(ipoint) .EQ. ns-1) .AND. (i .EQ. n2)) THEN
1376	!!$ !i=n2
1377	!!$ amhalf=diffgrup
1378	!!$ aphalf=diffup(ipoint)
1379	!!$ p=0.
1380	!!$! oxidation following Michaelis-Menten kinetics
1381	!!$ IF (uo(ipoint,i) .GT. 0.) THEN
1382	!!$! transform from partial pressure into concentration unit
1383	!!$ uo(ipoint,i)=uo(ipoint,i)/fpart(ipoint,i)
1384	!!$! Michaelis Menten Equation
1385	!!$ q=-rvmax(ipoint)*uo(ipoint,i)/(uo(ipoint,i)+rkm)
1386	!!$! temperature dependence (oxq10 from para.h)
1387	!!$ q=(q/oxq10)(oxq10*((sou(ipoint,i)-MAX(0.,tmean(ipoint)))/10.))
1388	!!$ goxid(ipoint)=goxid(ipoint)+q
1389	!!$! transform from concentration into partial pressure unit
1390	!!$ uo(ipoint,i)=uo(ipoint,i)*fpart(ipoint,i)
1391	!!$! in case of negative concentrations: no oxidation
1392	!!$! (can occur if number of time steps per day (nday) is small and water
1393	!!$! table varies strongly) (help: increase nday (or change numerical scheme!))
1394	!!$ ELSE
1395	!!$ q=0.
1396	!!$ goxid(ipoint)=goxid(ipoint)+q
1397	!!$ ENDIF
1398	!!$ ! transform from concentration into partial pressure unit
1399	!!$ q=q*fpart(ipoint,i)
1400	!!$ aold1=rkh*amhalf
1401	!!$ aold2=2.-rkh(amhalf+aphalf-h2p)
1402	!!$ aold3=rkh*aphalf
1403	!!$ aold4=2.rkq
1404	!!$ b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4
1405	!!$
1406	!!$ ELSEIF ((nw(ipoint) .EQ. ns-1) .AND. (i .EQ. n2+1)) THEN
1407	!!$ !i=n2+1
1408	!!$ amhalf=diffup(ipoint)
1409	!!$ aphalf=diffair
1410	!!$ p=0.
1411	!!$ q=0.
1412	!!$ goxid(ipoint)=goxid(ipoint)+q
1413	!!$ aold1=rkh*amhalf
1414	!!$ aold2=2.-rkh(amhalf+aphalf-h2p)
1415	!!$ aold3=rkh*aphalf
1416	!!$ aold4=2.rkq
1417	!!$ b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4
1418	!!$
1419	!!$ ELSEIF ((nw(ipoint) .EQ. ns-1) .AND. (i .GE. n2+2) .AND. (i .LE. n3-1)) THEN
1420	!!$
1421	!!$ !DO i=n2+2,n3-1
1422	!!$ amhalf=diffair
1423	!!$ aphalf=diffair
1424	!!$ p=0.
1425	!!$ q=0.
1426	!!$ aold1=rkh*amhalf
1427	!!$ aold2=2.-rkh(amhalf+aphalf-h2p)
1428	!!$ aold3=rkh*aphalf
1429	!!$ aold4=2.rkq
1430	!!$ b(ipoint,i)=aold1uo(ipoint,i-1)+aold2uo(ipoint,i)+aold3*uo(ipoint,i+1)+aold4
1431	!!$ !ENDDO
1432	!!$
1433	!!$ ENDIF
1434	!!$ ENDDO
1435	!!$ ENDDO
1436
1437	!4eme CAS: (nw .GE. ns): non traite! ma water table ne va jamais au-dessus de la
1438	!surface en ce qui me concerne. de toute facon, jamais d oxydation dans Walter dans
1439	!la colonne d eau au-dessus de la surface: cala modifie juste le transport
1440
1441
1442
1443	!DO ipoint =1,npts
1444	! IF (iday .EQ. 2) THEN
1445	! WRITE(,) 'matrice b apres inner point apres cas',b(ipoint,:)
1446	! ENDIF
1447	!ENDDO
1448
1449	! right boundary
1450	DO ipoint = 1,npts
1451	amhalf=diffair
1452	aphalf=diffair
1453	p=0.
1454	q=0.
1455	aoldm2=rkh*amhalf
1456	aoldm1=2.-rkh(amhalf+aphalf-h2p)
1457	aoldn=2.(rkhaphalfcatm+rkq)
1458	b(ipoint,n3)=aoldm2uo(ipoint,n3-1)+aoldm1uo(ipoint,n3)+aoldn
1459	ENDDO
1460
1461
1462
1463
1464	!*********************************************************
1465	! solution of tridiagonal system using subroutine tridag
1466	! calculate NEW concentration profiles uo(i)
1467
1468
1469	!DO i=n0,n3
1470	DO i=1,n3
1471	DO ipoint =1,npts
1472	IF (i .GE. n0(ipoint)) THEN
1473	tria(ipoint,i-n0(ipoint)+1)=0.
1474	trib(ipoint,i-n0(ipoint)+1)=0.
1475	tric(ipoint,i-n0(ipoint)+1)=0.
1476	trir(ipoint,i-n0(ipoint)+1)=0.
1477	triu(ipoint,i-n0(ipoint)+1)=0.
1478	ENDIF
1479	ENDDO
1480	enddo
1481
1482
1483	!do i=n0,n3
1484	DO i=1,n3
1485	DO ipoint =1,npts
1486	IF (i .GE. n0(ipoint)) THEN
1487	trib(ipoint,i-n0(ipoint)+1)=a(ipoint,i,i)
1488	ENDIF
1489	ENDDO
1490	ENDDO
1491
1492	!do i=n0+1,n3
1493	DO i=1,n3
1494	DO ipoint =1,npts
1495	IF (i .GE. n0(ipoint)+1) THEN
1496	tria(ipoint,i-n0(ipoint)+1)=a(ipoint,i,i-1)
1497	ENDIF
1498	ENDDO
1499	ENDDO
1500
1501	!do i=n0,n3-1
1502	DO i=1,n3-1
1503	DO ipoint =1,npts
1504	IF (i .GE. n0(ipoint)) THEN
1505	tric(ipoint,i-n0(ipoint)+1)=a(ipoint,i,i+1)
1506	ENDIF
1507	ENDDO
1508	ENDDO
1509
1510	! do i=n0,n3
1511	DO i=1,n3
1512	DO ipoint =1,npts
1513	IF (i .GE. n0(ipoint)) THEN
1514	trir(ipoint,i-n0(ipoint)+1)=b(ipoint,i)
1515	ENDIF
1516	ENDDO
1517	ENDDO
1518
1519	DO ipoint =1,npts
1520	ndim(ipoint)=n3-n0(ipoint)+1
1521	ENDDO
1522
1523	!WRITE(,) ' '
1524	!CALL tridiag(tria,trib,tric,trir,triu,ndim)
1525	!REMPLACEMENT DE LA ROUTINE TRIDAG
1526	!WRITE(,) ' '
1527
1528	DO ipoint =1,npts
1529	bet(ipoint)=trib(ipoint,1)
1530	triu(ipoint,1)=trir(ipoint,1)/bet(ipoint)
1531	ENDDO
1532
1533	!DO j=2,ndim
1534	!enleve j=2,500 et met a la place j=2,n3
1535	DO j=2,n3
1536	DO ipoint =1,npts
1537	IF (j .LE. ndim(ipoint)) THEN
1538	gam(ipoint,j)=tric(ipoint,j-1)/bet(ipoint)
1539	! bet(ipoint)=trib(ipoint,j)-tria(ipoint,j)*gam(ipoint,j)
1540	IF ((trib(ipoint,j)-tria(ipoint,j)*gam(ipoint,j)) .EQ. 0) THEN
1541	bet(ipoint)=trib(ipoint,j)-tria(ipoint,j)*gam(ipoint,j) + 10
1542	ELSE
1543	bet(ipoint)=trib(ipoint,j)-tria(ipoint,j)*gam(ipoint,j)
1544	ENDIF
1545	triu(ipoint,j)=(trir(ipoint,j)-tria(ipoint,j)*triu(ipoint,j-1))/bet(ipoint)
1546	ENDIF
1547	ENDDO
1548	ENDDO
1549
1550	!DO j=ndim-1,1,-1
1551	DO j=n3-1,1,-1
1552	DO ipoint =1,npts
1553	IF (j .LE. ndim(ipoint)-1) THEN
1554	triu(ipoint,j)=triu(ipoint,j)-gam(ipoint,j+1)*triu(ipoint,j+1)
1555	ENDIF
1556	ENDDO
1557	ENDDO
1558	!FIN DE REMPLACEMENT DE LA ROUTINE TRIDAG
1559
1560
1561	!DO i=n0,n3
1562	DO i=1,n3
1563	DO ipoint =1,npts
1564	IF (i .GE. n0(ipoint)) THEN
1565	uo(ipoint,i)=triu(ipoint,i-n0(ipoint)+1)
1566	! if (due to small number of time steps per day and strong variation in
1567	! water table (i.e. large change in vertical profile of diffusion coefficient)
1568	! negative methane concentrations occur: set to zero
1569	! should usually not happen - can be avoided by increasing nday (from 24 to
1570	! 240 or 2400 or ...) or possibly using a different numerical scheme
1571	IF (uo(ipoint,i) .LT. 0.) THEN
1572	negconc(ipoint)=negconc(ipoint)-(uo(ipoint,i)/fpart(ipoint,i))
1573	ENDIF
1574	uo(ipoint,i)=MAX(0.,triu(ipoint,i-n0(ipoint)+1))
1575	! transform from partial pressure into concentration unit
1576	uo(ipoint,i)=uo(ipoint,i)/fpart(ipoint,i)
1577	ENDIF
1578	ENDDO
1579	ENDDO
1580
1581
1582
1583	!la nouvelle valeur de uo va etre utilisée comme valeur initiale au pas de temps suivant....
1584	!************************************
1585	!************************************
1586
1587
1588	!**************************************************************
1589
1590	! plant-mediated transport:
1591	!**************************
1592
1593	DO ipoint =1,npts
1594	tplant(ipoint)=0.
1595	ENDDO
1596
1597	!DO i=n0,n3
1598	DO i=1,n3
1599	DO ipoint =1,npts
1600	IF (i .GE. n0(ipoint)) THEN
1601	! calculated fraction of methane that enters plants:
1602	! tveg: vegetation type from Sread.f
1603	! dveg: scaling factor from para.h
1604	! fgrow: growing stage of plants as calculated above (LAI)
1605	! root(i): vertical root distribution as calculated above
1606	cplant=MIN(tveg(ipoint)dvegfgrow(ipoint)root(ipoint,i)10.*rk,1.)
1607	! amount (concentration) of methane that enters plants
1608	cplant=cplant*uo(ipoint,i)
1609	dplant(ipoint)=dplant(ipoint)+cplant
1610	! amount (concentration) of methane that remains in soil
1611	uo(ipoint,i)=uo(ipoint,i)-cplant
1612	! only fraction (1.-pox) of methane entering plants is emitted
1613	! to the atmosphere - the rest (pox) is oxidized in the rhizosphere
1614	tplant(ipoint)=tplant(ipoint)+cplant*(1.-pox)
1615	ENDIF
1616	ENDDO
1617	ENDDO
1618	! transform amount (concentration) into flux unit (mg/m^2*d)
1619	flplant(:,iday)=tplant(:)*funit
1620
1621	!DO ipoint =1,npts
1622	! IF (iday .EQ. 1) THEN
1623	! WRITE(,) 'tveg',tveg,'fgrow',fgrow
1624	! WRITE(,) 'root',root
1625	! ENDIF
1626	!ENDDO
1627
1628	!DO ipoint =1,npts
1629	! IF (iday .EQ. 1) THEN
1630	! WRITE(,) 'matrice uo apres plante',uo
1631	! ENDIF
1632	! ENDDO
1633
1634
1635	! bubble transport:
1636	!******************
1637
1638	DO ipoint =1,npts
1639	tmax(ipoint)=0.
1640	ENDDO
1641
1642	!DO i=n0,nw
1643	DO i=1,n2
1644	DO ipoint =1,npts
1645	IF ((i .GE. n0(ipoint)) .AND. (i .LE. nw(ipoint))) THEN
1646	! methane in water: if concentration above threshold concentration rcmax
1647	! (calculated in Scalc.f): bubble formation
1648	cdiff=uo(ipoint,i)-rcmax(ipoint)
1649	cbub=MAX(0.,cdiff)
1650	tmax(ipoint)=tmax(ipoint)+cbub
1651	uo(ipoint,i)=MIN(uo(ipoint,i),rcmax(ipoint))
1652	ENDIF
1653	ENDDO
1654	ENDDO
1655
1656
1657	!test Estelle
1658	! IF (tmax .GT. 0.) then
1659	! WRITE(,) 'check', tmax, n0, nw, uo(nw), cdiff
1660	! ENDIF
1661	! if water table is above soil surface: bubble flux is emitted directly
1662	! to atmosphere (mg/m^2*d)
1663
1664	DO ipoint =1,npts
1665	IF (nw(ipoint) .GE. ns) THEN
1666	flbub(ipoint,iday)=tmax(ipoint)*funit
1667	flbubdiff(ipoint,iday)=0.
1668	tmax(ipoint)=0.
1669	! if water table is below soil surface:
1670	ELSE
1671	flbub(ipoint,iday)=0.
1672	flbubdiff(ipoint,iday)=0.
1673	! 1. if water table is 5 or less cm below soil surface: flbubdiff is
1674	! calculated in the same way as flbub and added to the total flux (mg/m^2*d)
1675	IF (nw(ipoint) .GE. ns-5) THEN
1676	flbubdiff(ipoint,iday)=tmax(ipoint)*funit
1677	tmax(ipoint)=0.
1678	ENDIF
1679	ENDIF
1680	ENDDO
1681
1682	!test Estelle
1683	! IF (flbub(iday) .gt. 0.) then
1684	! WRITE(,) 'buble', iday,flbub(iday),tmax, cdiff, cbub, nw
1685	! ENDIF
1686
1687
1688	! 2. if water table is more than 5 cm below soil surface: no bubble flux
1689	! is directly emitted to the atmosphere, BUT methane amount from bubble is
1690	! added to the first layer above the soil surface
1691	! (the distinction more/less than 5 cm below soil has been made, because
1692	! adding tmax to uo(nw+1) if the water table is less than 5cm below the soil
1693	! surface 'disturbed' the system too much and more time steps per day (larger
1694	! nday) is needed) (this solution avoids this)
1695	DO ipoint =1,npts
1696	uo(ipoint,nw(ipoint)+1)=uo(ipoint,nw(ipoint)+1)+tmax(ipoint)
1697	ENDDO
1698
1699	!DO ipoint =1,npts
1700	! IF (iday .EQ. 1) THEN
1701	! WRITE(,) 'matrice uo apres bubble',uo
1702	! ENDIF
1703	!ENDDO
1704
1705	! calculate diffusive flux from the concentration gradient (mg/m^2*d):
1706	!**********************************************************************
1707	DO ipoint =1,npts
1708	fdifftime(ipoint,iday)=((uo(ipoint,n3-2)-uo(ipoint,n3-1))/0.1)38.4diffair
1709	ENDDO
1710
1711	! transform from concentration into partial pressure unit
1712	!do i=n0,n3
1713	DO i=1,n3
1714	DO ipoint =1,npts
1715	IF (i .GE. n0(ipoint)) THEN
1716	uo(ipoint,i)=uo(ipoint,i)*fpart(ipoint,i)
1717	ENDIF
1718	ENDDO
1719	enddo
1720
1721
1722
1723	! DO ipoint =1,npts
1724	! IF (iday .EQ. 1) THEN
1725	! WRITE(,) 'matrice uo FIN TIMELOOP',uo
1726	! ENDIF
1727	! ENDDO
1728
1729	!
1730	! end of (iday=1,nday) loop, i.e. time step per day
1731	!***************************************************
1732
1733	ENDDO
1734
1735	!***************************************************
1736
1737
1738	! water-air interface
1739	! transform from partial pressure into concentration unit
1740	!do i=n0,n3
1741	DO i=1,n3
1742	DO ipoint =1,npts
1743	IF (i .GE. n0(ipoint)) THEN
1744	uo(ipoint,i)=uo(ipoint,i)/fpart(ipoint,i)
1745	ENDIF
1746	ENDDO
1747	ENDDO
1748
1749
1750	! calculate conc kept in soil -
1751	! necessary to calculate fluxdiff (diffusive flux from mass balance)
1752	DO i=1,nvert
1753	DO ipoint =1,npts
1754	uout(ipoint,i)=uo(ipoint,i)
1755	ENDDO
1756	ENDDO
1757
1758	!do i=n0,n3
1759	DO i=1,n3
1760	DO ipoint =1,npts
1761	IF (i .GE. n0(ipoint)) THEN
1762	tout(ipoint)=tout(ipoint)+(uout(ipoint,i)-uold2(ipoint,i))
1763	ENDIF
1764	ENDDO
1765	ENDDO
1766
1767	!DO ipoint =1,npts
1768	! tout(ipoint)=0.
1769	!ENDDO
1770
1771	tout(:)=0.
1772
1773	!do i=n0,n3
1774	DO i=1,n3
1775	DO ipoint =1,npts
1776	IF (i .GE. n0(ipoint)) THEN
1777	tout(ipoint)=tout(ipoint)+(uout(ipoint,i)-uold2(ipoint,i))
1778	ENDIF
1779	ENDDO
1780	ENDDO
1781
1782
1783	DO i=1,nvert
1784	DO ipoint =1,npts
1785	uold2(ipoint,i)=uout(ipoint,i)
1786	ENDDO
1787	ENDDO
1788
1789
1790	!*****************************
1791	!vire concout
1792
1793	! write CH4conc to varible concout after each time interval (at iday=nday)
1794	! concout(i,itime) is needed to write CH4conc profiles in Soutput.f
1795	! do i=1,n3
1796	! concout(i,itime)=uout(i)
1797	! enddo
1798	! do i=n3+1,n
1799	! transform from partial pressure into concentration unit
1800	! concout(i,itime)=catm/fpart(i)
1801	! enddo
1802
1803	! calculate daily means of bubble flux, fluxbubdiff, plant flux, and
1804	! diffusive flux (from concentration gradient) (mg/m^2*d)
1805
1806	!*******************************
1807	!la j enleve la dimension temps des sorties...
1808	!*******************************
1809
1810
1811
1812	DO ipoint =1,npts
1813	fluxbub(ipoint)=0.
1814	fluxbubdiff(ipoint)=0.
1815	fluxplant(ipoint)=0.
1816	fdiffday(ipoint)=0.
1817	ENDDO
1818
1819
1820	DO i=1,nday
1821	DO ipoint =1,npts
1822	fluxbub(ipoint)=fluxbub(ipoint)+flbub(ipoint,i)
1823	fluxbubdiff(ipoint)=fluxbubdiff(ipoint)+flbubdiff(ipoint,i)
1824	fluxplant(ipoint)=fluxplant(ipoint)+flplant(ipoint,i)
1825	fdiffday(ipoint)=fdiffday(ipoint)+fdifftime(ipoint,i)
1826	ENDDO
1827	ENDDO
1828
1829	fluxbub(:)=fluxbub(:)/nday
1830	fluxbubdiff(:)=fluxbubdiff(:)/nday
1831	fluxplant(:)=fluxplant(:)/nday
1832	fdiffday(:)=fdiffday(:)/nday
1833
1834
1835	!
1836	! calculate diffusive flux ( from mass balance):
1837	!************************************************
1838	! use amounts in concentration unit and tranform in flux unit (mg/m^2*d)
1839	! wpro: amount of methane produced during entire day
1840	! goxid: amount of methane oxidized during entire day
1841	! mutiply by rk (depending on number of time steps per day)
1842	! dplant: amount of methane that entered plants during entire day
1843	! (includes both: methane emitted through plants AND oxidized in rhizosphere)
1844	! tout: additional methane stored in soil during entire day
1845	! wtdiff: see above
1846	! fluxbub: bubble flux (if water table > soil surface)
1847	! bubble flux (if water table less than 5 cm below soil surface (see above))
1848	! negconc: see above
1849	fluxdiff(:)= ((rk(wpro(:)+goxid(:)) - dplant(:)-tout(:)-wtdiff)/nday)funit &
1850	& - fluxbub(:) - fluxbubdiff(:) + negconc(:)*(funit/nday)
1851
1852	! fluxbubdiff is added to diffusive flux (mg/m^2*d)
1853
1854
1855	fluxdiff(:) = fluxdiff(:)+fluxbubdiff(:)
1856	fdiffday(:) = fdiffday(:)+fluxbubdiff(:)
1857
1858
1859	!enleve boucle if sur ihump
1860
1861	! ELSE
1862	! fluxdiff=0.
1863	! fdiffday=0.
1864	! fluxbub=0.
1865	! fluxplant=0.
1866	! fluxtot=0.
1867	! !do i=1,n
1868	! ! concout(i,itime)=0.
1869	! !enddo
1870	!
1871	!! end of (if (ijump .eq. 1)) 'loop'
1872	!!***********************************
1873	! ENDIF
1874
1875
1876
1877
1878
1879	! to calculate wtdiff (see above)
1880	! define uwp1
1881
1882	!non necessaire: la WTD ne se modifie pas d un pas de tps a l autre
1883
1884	! uwp0=uo(ipoint,nw)
1885	! uwp1=uo(ipoint,nw+1)
1886	! uwp2=uo(ipoint,nw+2)
1887
1888
1889
1890	! calculate total flux fluxtot (using diffusive flux from mass balance)
1891	! and sum of diffusive flux from mass balance and bubble flux (for plots)
1892	! (mg/m^2*d)
1893	! do i=1,ntime
1894
1895	fluxtot(:)=fluxdiff(:)+fluxbub(:)+fluxplant(:)
1896	fbubdiff(:)=fluxdiff(:)+fluxbub(:)
1897
1898
1899	DO ipoint =1,npts
1900	ch4_flux_density_dif(ipoint) = fluxdiff(ipoint)
1901	ch4_flux_density_bub(ipoint) = fluxbub(ipoint)
1902	ch4_flux_density_pla(ipoint) = fluxplant(ipoint)
1903	ENDDO
1904
1905
1906	WHERE (fluxtot(:) .LE. 0.0)
1907	ch4_flux_density_tot(:) = 0.0
1908	ELSEWHERE
1909	ch4_flux_density_tot(:) = fluxtot(:)
1910	endwhere
1911
1912	!pour contrebalancer valeur par defaut mise pour les grid-cells ou veget_max=0
1913	WHERE ((veget_max(:,1) .LE. 0.95) .AND. (sum_veget_nat_ss_nu(:) .GT. 0.1))
1914	ch4_flux_density_tot(:) = ch4_flux_density_tot(:)
1915	ELSEWHERE
1916	ch4_flux_density_tot(:) = 0.0
1917	ENDWHERE
1918
1919
1920
1921	END SUBROUTINE ch4_wet_flux_density_wet3
1922	!********************************************************************
1923	!!$
1924	!!$SUBROUTINE tridag(a,b,c,r,u,n)
1925	!!$!
1926	!!$! Press et al., Numerical Recipes, FORTRAN, p.43
1927	!!$!************************************************
1928	!!$
1929	!!$ ! Domain size
1930	!!$ INTEGER(i_std), INTENT(in) :: n
1931	!!$ REAL(r_std),DIMENSION(n), INTENT(in) :: a
1932	!!$ REAL(r_std),DIMENSION(n), INTENT(in) :: b
1933	!!$ REAL(r_std),DIMENSION(n), INTENT(in) :: c
1934	!!$ REAL(r_std),DIMENSION(n), INTENT(in) :: r
1935	!!$ REAL(r_std),DIMENSION(n), INTENT(out) :: u
1936	!!$ INTEGER(stnd),PARAMETER :: nmax=500
1937	!!$ INTEGER(i_std) :: j
1938	!!$ REAL(r_std),DIMENSION(nmax) :: gam
1939	!!$ REAL(r_std) :: bet
1940	!!$
1941	!!$ bet=b(1)
1942	!!$
1943	!!$ u(1)=r(1)/bet
1944	!!$
1945	!!$ do j=2,n
1946	!!$ gam(j)=c(j-1)/bet
1947	!!$ bet=b(j)-a(j)*gam(j)
1948	!!$ u(j)=(r(j)-a(j)*u(j-1))/bet
1949	!!$ enddo
1950	!!$
1951	!!$ do j=n-1,1,-1
1952	!!$ u(j)=u(j)-gam(j+1)*u(j+1)
1953	!!$ enddo
1954	!!$
1955	!!$END SUBROUTINE tridag
1956
1957
1958	END MODULE stomate_wet_ch4_pt_ter_wet3

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd_mict_peat_ch4/src_stomate/stomate_wet_ch4_pt_ter_wet3.f90 @ 7346

Download in other formats: