Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

stomate_wet_ch4_pt_ter_0.f90 @ 5491

Last change on this file since 5491 was 5080, checked in by chunjing.qiu, 6 years ago
soil freezing, soil moisture, fwet bugs fixed
File size: 66.4 KB

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

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

Context Navigation

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

Download in other formats: