Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

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

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_wet1.f90 @ 7346

Download in other formats: