Context Navigation

stomate_wet_ch4_pt_ter_wet.f90 @ 6737

Last change on this file since 6737 was 4719, checked in by albert.jornet, 7 years ago

Merge: from revisions [4491:4695/trunk/ORCHIDEE]

Merge done in [4671:4718/perso/albert.jornet/MICT_MERGE]

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

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

Context Navigation

source: branches/publications/ORCHIDEE_GLUC_r6545/src_stomate/stomate_wet_ch4_pt_ter_wet.f90 @ 6737

Download in other formats: