Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

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

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_wet4.f90 @ 5491

Download in other formats: