Context Navigation

isotrac_routines_mod.F90 @ 4582

Last change on this file since 4582 was 4491, checked in by crisi, 19 months ago
Bug corrections in LMDZiso, especially for water tagging
Property svn:executable set to ``*
File size: 94.4 KB

Line
1	#ifdef ISO
2	#ifdef ISOTRAC
3	! $Id: $
4
5	MODULE isotrac_routines_mod
6	! on créé ce module pour éviter dépendances circulaires.
7	! isotopes_routines_mod a besoin de isotrac et isotopes_verif
8	! isotopes_verif a besoin de isotopes et isotrac
9	! isotrac n'a besoin que de isotopes
10	USE infotrac_phy, ONLY: ntraciso=>ntiso, niso, index_trac=>itZonIso, ntraceurs_zone=>nzone
11	IMPLICIT NONE
12
13	CONTAINS
14
15
16	subroutine uncompress_commun_zone(ncas, cas, &
17	& xtp_cas,xtp,xtwater_cas,xtwater,xtevap_cas,xtevap, &
18	& ncum,izone)
19
20	USE isotopes_mod, ONLY: ridicule,iso_eau
21
22	implicit none
23
24	! decompression des outputs communs à tous les cas dans
25	! appel_stewart
26	! cas des traceurs. Ici, aucun risque de revap franche.
27
28	! inputs
29	integer ncas,ncum
30	integer cas(ncum)
31	real xtevap_cas(niso,ncum)
32	real xtp_cas(niso,ncum)
33	real xtwater_cas(niso,ncum)
34	integer izone
35	! outputs
36	real xtwater(ntraciso,ncum)
37	real xtp(ntraciso,ncum)
38	real xtevap(ntraciso,ncum)
39
40	! locals
41	integer il,ixt,iiso,ixt_revap
42
43
44	do il=1,ncas
45	do iiso=1,niso
46	ixt=index_trac(izone,iiso)
47	xtevap(ixt,cas(il))=xtevap_cas(iiso,il)
48	xtp(ixt,cas(il))=xtp_cas(iiso,il)
49	xtwater(ixt,cas(il))=xtwater_cas(iiso,il)
50	enddo !do iiso=1,niso
51	enddo !do il=1,ncas
52
53	end subroutine uncompress_commun_zone
54
55
56	subroutine uncompress_commun_zone_revap(ncas, cas, &
57	& xtp_cas,xtp,xtwater_cas,xtwater, &
58	& xtevap_cas,xtevap, &
59	& ncum,izone,Eqi_stewart,fac_ftmr_cas, &
60	#ifdef ISOVERIF
61	& Exi_cas,Exi, &
62	#endif
63	& xtp_avantevap_cas,liq,hdiag)
64
65	USE isotopes_mod, ONLY: ridicule,iso_eau,iso_HDO,ridicule_evap
66	USE isotrac_mod, only: option_revap,evap_franche,izone_revap, &
67	& ridicule_trac
68	#ifdef ISOVERIF
69	USE isotopes_verif_mod
70	#endif
71
72	implicit none
73
74	! decompression des outputs communs à tous les cas dans
75	! appel_stewart
76	! cas des traceurs: mais ici, risque de révap franche -> on fat
77	! plus attention
78
79	! inputs
80	integer ncas,ncum
81	integer cas(ncum)
82	real xtevap_cas(niso,ncum)
83	real xtp_cas(niso,ncum)
84	real xtwater_cas(niso,ncum)
85	integer izone
86	real Eqi_stewart(ncum)
87	real xtp_avantevap_cas(niso,ncum)
88	real fac_ftmr_cas(ncum)
89	integer liq
90	real hdiag(ncas)
91
92	! outputs
93	real xtwater(ntraciso,ncum)
94	real xtp(ntraciso,ncum)
95	real xtevap(ntraciso,ncum)
96
97	! locals
98	integer il,ixt,iiso,ixt_revap
99	real xtaddp_tag(niso,ncum)
100	#ifdef ISOVERIF
101	integer ieau,iHDO
102	real Exi_cas(niso,ncum)
103	real Exi(ntraciso,ncum)
104	!USE isotopes_verif, ONLY:
105	#endif
106
107	! write(,) 'compress_stewart 315 tmp: ',
108	! : 'entrée dans uncompress_commun_zone_revap'
109
110	do il=1,ncas
111	if ((option_revap.eq.1).and. &
112	& ((((liq.eq.1).and. &
113	& (Eqi_stewart(il)*fac_ftmr_cas(il).gt.evap_franche) &
114	& .and.(hdiag(il).lt.0.99)).or. &
115	& ((liq.eq.0).and. &
116	& (Eqi_stewart(il)*fac_ftmr_cas(il).ge.0.0))).or. &
117	& (izone.eq.izone_revap))) then
118	! if ((option_revap.eq.1).and.
119	! : (((Eqi_stewart(il)*fac_ftmr_cas(il).gt.evap_franche)
120	! : .or.((liq.eq.0).and.
121	! : (Eqi_stewart(il)*fac_ftmr_cas(il).ge.0.0))).or.
122	! : (izone.eq.izone_revap))) then
123	! on met la revap dans izone_revap si option_revap=1 et si:
124	! * evap glace (non fractionnante)
125	! * ou evap liq suffisemment forte pour que pas de flux
126	! d'isotopes négatifs.
127	! si option_revap=1 et izone=izone_revap, on met aussi dans izone_revap
128	!#ifdef ISOVERIF
129	! if (il.eq.1) then
130	! write(,) 'compress tmp 341: revap dans izone_revap'
131	! write(,) 'Eqi_stewart(il),fac_ftmr_cas(il),evap_franche=',
132	! : Eqi_stewart(il),fac_ftmr_cas(il),evap_franche
133	! write(,) 'il,xtp_cas(iso_eau,il)=',il,xtp_cas(iso_eau,il)
134	! write(,) 'il,xtp_avantevap_cas(iso_eau,il)=',il,
135	! : xtp_avantevap_cas(iso_eau,il)
136	! write(,) 'xtp(ixt_revap,cas(il))=',
137	! : xtp(index_trac(izone_revap,iso_eau),cas(il))
138	! endif
139	!#endif
140
141	! toute la révap franche va dans izone_revap
142	do iiso=1,niso
143	ixt=index_trac(izone,iiso)
144	ixt_revap=index_trac(izone_revap,iiso)
145	! le terme d'évap qui était pour la zone izone devient
146	! nul, et à la place on le met dans izone_revap&
147	xtevap(ixt_revap,cas(il))=xtevap(ixt_revap,cas(il)) &
148	& +xtevap_cas(iiso,il)
149
150	! ce qui a été ajouté à xtp par rapport à xtp_avantevap
151	! est ajouté à izone_revap, au lieu de izone
152	xtaddp_tag(iiso,il)=xtp_cas(iiso,il) &
153	& -xtp_avantevap_cas(iiso,il)
154	xtp(ixt_revap,cas(il))= &
155	& xtp(ixt_revap,cas(il)) &
156	& +xtaddp_tag(iiso,il)
157	enddo !do iiso=1,niso
158	#ifdef ISOVERIF
159	if (iso_HDO.gt.0) then
160	if (xtevap_cas(iso_eau,il).gt.ridicule_evap) then
161	call iso_verif_aberrant_choix( &
162	& xtevap_cas(iso_HDO,il),xtevap_cas(iso_eau,il), &
163	& ridicule_trac,deltalimtrac*2,'compress 344a')
164	endif
165	ieau=index_trac(izone_revap,iso_eau)
166	iHDO=index_trac(izone_revap,iso_HDO)
167	call iso_verif_aberrant_choix(xtevap(iHDO,cas(il)), &
168	& xtevap(ieau,cas(il)),ridicule_trac,deltalimtrac*2, &
169	& 'compress 344b')
170	endif !if (iso_HDO.gt.0) then
171	#endif
172
173	! l'évap des autres zones devient nulle
174	if (izone.ne.izone_revap) then
175	do iiso=1,niso
176	ixt=index_trac(izone,iiso)
177	xtevap(ixt,cas(il))=0.0
178	xtp(ixt,cas(il))=xtp_avantevap_cas(iiso,il)
179	enddo
180	endif
181
182	!#ifdef ISOVERIF
183	! if (il.eq.1) then
184	! write(,) 'compress tmp 341: revap dans izone_revap'
185	! write(,) 'xtp(ixt_revap,cas(il))=',
186	! : xtp(index_trac(izone_revap,iso_eau),cas(il))
187	! write(,) 'xtp(ixt,cas(il))=',
188	! : xtp(index_trac(izone,iso_eau),cas(il))
189	! endif
190	!#endif
191
192	else !if ((Eqi_stewart(il).gt.ridicule_evap*100)
193
194	do iiso=1,niso
195	ixt=index_trac(izone,iiso)
196	xtevap(ixt,cas(il))=xtevap_cas(iiso,il)
197	xtp(ixt,cas(il))=xtp_cas(iiso,il)
198	enddo !do iiso=1,niso
199	endif !if ((Eqi_stewart(il).gt.ridicule_evap*100)
200	enddo !do il=1,ncas
201
202	do il=1,ncas
203	do iiso=1,niso
204	ixt=index_trac(izone,iiso)
205	xtwater(ixt,cas(il))=xtwater_cas(iiso,il)
206	#ifdef ISOVERIF
207	Exi(ixt,cas(il))=Exi_cas(iiso,il)
208	#endif
209	enddo !do iiso=1,niso
210	enddo !do il=1,ncas
211	! write(,) 'compress_stewart 361 tmp: ',
212	! : 'sortie de uncompress_commun_zone_revap'
213
214	end subroutine uncompress_commun_zone_revap
215
216
217
218	subroutine compress_cond_facftmr_zone( &
219	& ncas, cas, &
220	& Eqi_prime_cas,Eqi_prime, &
221	& Pqisup_cas,Pqisup, &
222	& Pxtisup_cas,Pxtisup, &
223	& qp_avantevap_cas,qp_avantevap, &
224	& xtp_avantevap_cas,xtp_avantevap, &
225	& xtevapsup_cas,xtevap,&
226	& water_cas,water,&
227	#ifdef ISOVERIF
228	& evap_cas,evap, &
229	#endif
230	& nloc,ncum,nd,i,izone)
231
232	USE isotopes_mod, ONLY: iso_eau
233	#ifdef ISOVERIF
234	USE isotopes_verif_mod
235	#endif
236
237	implicit none
238
239	! compression dans le cas condensation_facftmr
240	! inputs
241	integer nd,ncum,nloc
242	integer ncas
243	integer cas(ncum)
244	integer i
245	integer izone
246	real &
247	& xtevapsup_cas(niso,ncum),xtevap(ntraciso,ncum), &
248	& water_cas(ncum),water(ncum)
249	real Eqi_prime_cas(ncum),Eqi_prime(ncum), &
250	& Pqisup_cas(ncum),Pqisup(ncum), &
251	& Pxtisup_cas(niso,ncum),Pxtisup(ntraciso,ncum), &
252	& qp_avantevap_cas(ncum),qp_avantevap(ncum), &
253	& xtp_avantevap_cas(niso,ncum),xtp_avantevap(ntraciso,ncum)
254	#ifdef ISOVERIF
255	real evap_cas(ncum),evap(ncum)
256	#endif
257	! locals
258	integer il,ixt,ieau,iiso
259
260	ieau=index_trac(izone,iso_eau)
261	do il=1,ncas
262	if (qp_avantevap(cas(il)).gt.0.0) then
263	Eqi_prime_cas(il)=Eqi_prime(cas(il)) &
264	& *(xtp_avantevap(ieau,cas(il))/qp_avantevap(cas(il)))
265	else !if (qp_avantevap_cas(cas(il)).gt.0.0) then
266	#ifdef ISOVERIF
267	call iso_verif_egalite_choix( &
268	& (Eqi_prime(cas(il))),0.0, &
269	& 'compress_stewart 495',errmax,errmaxrel)
270	#endif
271	Eqi_prime_cas(il)=0.0
272	endif !if (qp_avantevap_cas(cas(il)).gt.0.0) thens
273
274	if (Pqisup(cas(il))-Eqi_prime(cas(il)).gt.0.0) then
275	water_cas(il)=water(cas(il)) &
276	& *((Pxtisup(ieau,cas(il))-Eqi_prime_cas(il)) &
277	& /(Pqisup(cas(il))-Eqi_prime(cas(il))))
278	else !if (Pqisup(cas(il)).gt.0.0) then
279	#ifdef ISOVERIF
280	call iso_verif_egalite_choix(water(cas(il)),0.0, &
281	& 'compress_stewart 507',errmax,errmaxrel)
282	#endif
283	water_cas(il)=0.0
284	endif !if (Pqisup(cas(il)).gt.0.0) then
285
286	Pqisup_cas(il)=Pxtisup(ieau,cas(il))
287	qp_avantevap_cas(il)=xtp_avantevap(ieau,cas(il))
288
289	!#ifdef ISOVERIF
290	! call iso_verif_noNaN(water_cas(il),'compress_stewart 518')
291	! evap_cas(il)=evap(cas(il))
292	! : *(xtp_avantevap(ieau,cas(il))/qp_avantevap(cas(il)))
293	!! qp_cas(il)=xtp(ieau,cas(il))
294	!#endif
295	do iiso=1,niso
296	ixt=index_trac(izone,iiso)
297	Pxtisup_cas(iiso,il)=Pxtisup(ixt,cas(il))
298	xtp_avantevap_cas(iiso,il)=xtp_avantevap(ixt,cas(il))
299	xtevapsup_cas(iiso,il)=xtevap(ixt,cas(il))
300	enddo
301	enddo
302
303	end subroutine compress_cond_facftmr_zone
304
305	subroutine compress_cond_nofftmr_zone( &
306	& ncas, cas, &
307	& Eqi_prime_cas,Eqi_prime, &
308	& Pqisup_cas,Pqisup, &
309	& Pxtisup_cas,Pxtisup, &
310	& water_cas,water, &
311	& qp_avantevap_cas,qp_avantevap, &
312	& xtp_avantevap_cas,xtp_avantevap, &
313	& xt_cas,xt,q_cas,q, &
314	& xtevapsup_cas,xtevap, &
315	#ifdef ISOVERIF
316	& evap_cas,evap, &
317	#endif
318	& nloc,ncum,nd,i,izone)
319
320	USE isotopes_mod, ONLY: iso_eau
321	#ifdef ISOVERIF
322	USE isotopes_verif_mod
323	#endif
324
325	implicit none
326
327	! compression dans le cas condensation_facftmr
328	! inputs
329	integer nloc,nd,ncum
330	integer ncas
331	integer cas(ncum)
332	integer i
333	integer izone
334	real &
335	& xt_cas(niso,ncum),q_cas(ncum),xt(ntraciso,ncum),q(ncum), &
336	& xtevapsup_cas(niso,ncum),xtevap(ntraciso,ncum), &
337	& water_cas(ncum),water(ncum)
338	real Eqi_prime_cas(ncum),Eqi_prime(ncum), &
339	& Pqisup_cas(ncum),Pqisup(ncum), &
340	& Pxtisup_cas(niso,ncum),Pxtisup(ntraciso,ncum), &
341	& qp_avantevap_cas(ncum),qp_avantevap(ncum), &
342	& xtp_avantevap_cas(niso,ncum),xtp_avantevap(ntraciso,ncum)
343	#ifdef ISOVERIF
344	real evap_cas(ncum),evap(ncum)
345	#endif
346	! locals
347	integer il,ixt,ieau,iiso
348
349	ieau=index_trac(izone,iso_eau)
350	do il=1,ncas
351	if (qp_avantevap(cas(il)).gt.0) then
352	Eqi_prime_cas(il)=Eqi_prime(cas(il)) &
353	& *(xtp_avantevap(ieau,cas(il))/qp_avantevap(cas(il)))
354	else
355	Eqi_prime_cas(il)=Eqi_prime(cas(il)) &
356	& *(xt(ieau,cas(il))/q(cas(il)))
357	endif
358
359	if (Pqisup(cas(il))-Eqi_prime(cas(il)).gt.0.0) then
360	water_cas(il)=water(cas(il)) &
361	& *((Pxtisup(ieau,cas(il))-Eqi_prime_cas(il)) &
362	& /(Pqisup(cas(il))-Eqi_prime(cas(il))))
363	else !if (Pqisup(cas(il)).gt.0.0) then
364	#ifdef ISOVERIF
365	call iso_verif_egalite_choix(water(cas(il)),0.0, &
366	& 'compress_stewart 654',errmax,errmaxrel)
367	#endif
368	water_cas(il)=0.0
369	endif !if (Pqisup(cas(il)).gt.0.0) then
370
371	Pqisup_cas(il)=Pxtisup(ieau,cas(il))
372	qp_avantevap_cas(il)=xtp_avantevap(ieau,cas(il))
373	q_cas(il)=xt(ieau,cas(il))
374	!#ifdef ISOVERIF
375	! if (qp_avantevap(cas(il)).gt.0.0) then
376	! evap_cas(il)=evap(cas(il))
377	! : *(xtp_avantevap(ieau,cas(il))/qp_avantevap(cas(il)))
378	! else
379	! evap_cas(il)=evap(cas(il))
380	! : *(xt(ieau,cas(il))/q(cas(il)))
381	! endif
382	! qp_cas(il)=xtp(ieau,cas(il))
383	!#endif
384	do iiso=1,niso
385	ixt=index_trac(izone,iiso)
386	Pxtisup_cas(iiso,il)=Pxtisup(ixt,cas(il))
387	xtp_avantevap_cas(iiso,il)=xtp_avantevap(ixt,cas(il))
388	xt_cas(iiso,il)=xt(ixt,cas(il))
389	xtevapsup_cas(iiso,il)=xtevap(ixt,cas(il))
390	enddo
391	enddo
392
393	end subroutine compress_cond_nofftmr_zone
394
395	subroutine compress_noevap_zone( &
396	& ncas, cas, &
397	& Pqisup_cas,Pqisup, &
398	& Pxtisup_cas,Pxtisup, &
399	& xtp_avantevap_cas,xtp_avantevap, &
400	& xtevapsup_cas,xtevap, &
401	& water_cas,water, &
402	#ifdef ISOVERIF
403	& evap_cas,evap, &
404	#endif
405	& nloc,ncum,nd,i,izone)
406
407	USE isotopes_mod, ONLY: ridicule,iso_eau
408	#ifdef ISOVERIF
409	USE isotopes_verif_mod
410	#endif
411	implicit none
412
413	! compression dans le cas condensation_facftmr
414	integer nloc,nd,ncum
415	integer ncas
416	integer cas(ncum)
417	integer i
418	integer izone
419	real xtevapsup_cas(niso,ncum),xtevap(ntraciso,ncum), &
420	& water_cas(ncum),water(ncum)
421	real Pqisup_cas(ncum),Pqisup(ncum), &
422	& Pxtisup_cas(niso,ncum),Pxtisup(ntraciso,ncum), &
423	& xtp_avantevap_cas(niso,ncum),xtp_avantevap(ntraciso,ncum)
424	#ifdef ISOVERIF
425	real evap_cas(ncum),evap(ncum)
426	#endif
427	integer il,ixt,ieau,iiso
428
429	ieau=index_trac(izone,iso_eau)
430	do il=1,ncas
431	Pqisup_cas(il)=Pxtisup(ieau,cas(il))
432	if (Pqisup(cas(il)).gt.0.0) then
433	water_cas(il)=water(cas(il)) &
434	& *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
435	else
436	water_cas(il)=0.0
437	#ifdef ISOVERIF
438	call iso_verif_egalite_choix(water(cas(il)), &
439	& 0.0,'compress_stewart 709',errmax,errmaxrel)
440	#endif
441	endif
442	#ifdef ISOVERIF
443	! evap_cas(il)=evap(cas(il))
444	! : *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
445	! qp_cas(il)=xtp(ieau,cas(il))
446	#endif
447	do iiso=1,niso
448	ixt=index_trac(izone,iiso)
449	Pxtisup_cas(iiso,il)=Pxtisup(ixt,cas(il))
450	xtp_avantevap_cas(iiso,il)=xtp_avantevap(ixt,cas(il))
451	xtevapsup_cas(iiso,il)=xtevap(ixt,cas(il))
452	enddo
453	enddo
454
455	end subroutine compress_noevap_zone
456
457	subroutine compress_evap_liq_zone(iflag_con,ncas, &
458	& cas, &
459	& Pqisup_cas,Pqisup, &
460	& Pxtisup_cas,Pxtisup, &
461	& xtp_avantevap_cas,xtp_avantevap, &
462	& xtp_avantevaptrac_cas,qp_avantevaptrac_cas, &
463	& xtevapsup_cas,xtevap, &
464	& water_cas,water, &
465	& Eqi_stewart,Pqiinf_stewart,Eqi_prime_cas,&
466	& Pqiinf,Eqi_par,Pqiinf_par,Eqi_prime,ptrac, &
467	& Eqi,Eqi_cas, &
468	! & qp_cas,
469	#ifdef ISOVERIF
470	& evap_cas,evap, &
471	#endif
472	& nloc,ncum,nd,izone)
473
474	USE isotopes_mod, ONLY: ridicule,iso_eau
475	#ifdef ISOVERIF
476	USE isotopes_verif_mod
477	#endif
478	implicit none
479
480	! compression dans le cas condensation_facftmr
481	! inputs et outputs
482	integer iflag_con
483	integer izone
484	integer nloc,nd,ncum
485	integer ncas
486	integer cas(ncum)
487	! integer i
488	real xtevapsup_cas(niso,ncum),xtevap(ntraciso,ncum)
489	real xtp_avantevap_cas(niso,ncum), &
490	& xtp_avantevap(ntraciso,ncum)
491	real water_cas(ncum),water(ncum)
492	real xtp_avantevaptrac_cas(niso,ncum), &
493	& qp_avantevaptrac_cas(ncum)
494	real ptrac(ncum)
495	! real qp_cas(ncum)
496	#ifdef ISOVERIF
497	real evap_cas(ncum),evap(ncum)
498	#endif
499	real &
500	& Eqi_stewart(ncum),Pqiinf_stewart(ncum),Eqi_prime_cas(ncum), &
501	& Pqiinf(ncum),Eqi_par(ncum),Pqiinf_par(ncum), &
502	& Eqi_prime(ncum),Pqisup(ncum),Pqisup_cas(ncum), &
503	& Pxtisup(ntraciso,ncum),Pxtisup_cas(niso,ncum), &
504	& Eqi(ncum),Eqi_cas(ncum)
505	! locals
506	integer il,ixt,iiso,ieau
507
508	! write(,) 'compress 910: xtp_avantevap(iso_eau,cas(1))=',
509	! : xtp_avantevap(iso_eau,cas(1))
510	! write(,) 'compress_evap_liq_zone 510: ncas,ncum=',ncas,ncum
511	ptrac(:)=0. ! CR 31 mars 2023: initialisation de ptrac
512
513	ieau=index_trac(izone,iso_eau)
514	do il=1,ncas
515	Pqisup_cas(il)=Pxtisup(ieau,cas(il))
516
517	if (Pqisup(cas(il)).gt.0.0) then
518	ptrac(il)=Pxtisup(ieau,cas(il))/Pqisup(cas(il))
519	Eqi_prime_cas(il)=Eqi_prime(cas(il)) &
520	& *ptrac(il)
521	else
522	#ifdef ISOVERIF
523	call iso_verif_egalite(( &
524	& Eqi_prime(cas(il))),0.0, &
525	& 'compress_stewart 979')
526	#endif
527	Eqi_prime_cas(il)=0.0
528	endif
529
530	if (Pqisup(cas(il))-Eqi_prime(cas(il)).gt.0.0) then
531	water_cas(il)=water(cas(il)) &
532	& *((Pxtisup(ieau,cas(il))-Eqi_prime_cas(il)) &
533	& /(Pqisup(cas(il))-Eqi_prime(cas(il))))
534	else !if (Pqisup(cas(il)).gt.0.0) then
535	#ifdef ISOVERIF
536	call iso_verif_egalite_choix(water(cas(il)),0.0, &
537	& 'compress_stewart 507',errmax,errmaxrel)
538	#endif
539	water_cas(il)=0.0
540	endif !if (Pqisup(cas(il)).gt.0.0) then
541
542
543	! qp_cas(il)=qp(cas(il))
544	#ifdef ISOVERIF
545	! evap_cas(il)=evap(cas(il))
546	! & *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
547	! ! ce calcul est faux& normalement,
548	! evap_cas(il)=Eqi_prime_cas(il)//100.0/delP_cas(il)/SIGDg2.0
549	! & -(ieau,cas(il))
550	#endif
551	qp_avantevaptrac_cas(il)=xtp_avantevap(ieau,cas(il))
552	do iiso=1,niso
553	ixt=index_trac(izone,iiso)
554	Pxtisup_cas(iiso,il)=Pxtisup(ixt,cas(il))
555	xtp_avantevap_cas(iiso,il)=xtp_avantevap(iiso,cas(il))
556	xtp_avantevaptrac_cas(iiso,il)=xtp_avantevap(ixt,cas(il))
557	xtevapsup_cas(iiso,il)=xtevap(ixt,cas(il))
558	enddo
559	enddo !do il=1,ncas
560	! write(,) 'compress_stewart 931: ',
561	! & 'xtp_avantevap_cas(iso_eau,1)=',
562	! & xtp_avantevap_cas(iso_eau,1)
563	! write(,) 'xtp_avantevap(iso_eau,cas(1))=',
564	! & xtp_avantevap(iso_eau,cas(1))
565
566	! calculs des flux de masses à mettre en argument de stewart:
567	! comme l'eau n'est pas bien concervée dans les ddfts, on est
568	! obligé de bidouillé.
569	! 1) soit on considère Pqisup, Eqi, et Pqiinf_par=Pqisup-Eqi
570	! et on suppose que dans la réalité les compositions de
571	! Pqiinf sont les même que Pqiinf_par
572	! 2) soit on considère Pqisup, Eqi_par=Pqisup-Pqiinf, et Pqiinf,
573	! et on suppose que dans la réalité les compositions de
574	! Eqi_prime sont les même que Eqi_par
575	do il=1,ncas
576	if (Pqisup(cas(il)).gt.0.0) then
577	if ((water(cas(il)).gt.ridicule/100).and. &
578	& (Pqiinf_par(cas(il)).le.0.0)) then
579	! on ne peut pas utiliser la méthode 1, car KE prédit de l'eau
580	! alors que le bilan de masse n'enprédit pas.
581	! Peut-on utiliser la méthode 2?
582	Pqiinf_stewart(il)=Pqiinf(cas(il))*ptrac(il)
583	Eqi_stewart(il)=Eqi_par(cas(il)) *ptrac(il)
584	else !if ((water(il,i).gt.ridicule/100).and.(Pqiinf_par.le.0.0)) then
585	! il n'y a pas d'obstacles à l'utilisation de 1)
586	Pqiinf_stewart(il)=Pqiinf_par(cas(il))*ptrac(il)
587	if (iflag_con.eq.30) then
588	Eqi_stewart(il)=Eqi_prime(cas(il))*ptrac(il)
589	else
590	Eqi_stewart(il)=Eqi(cas(il))*ptrac(il)
591	endif
592	endif !if ((water(il,i).gt.ridicule/100).and.(Pqiinf_par.le.0.0)) then
593	else ! if (Pqisup(cas(il).gt.0.0) then
594	#ifdef ISOVERIF
595	call iso_verif_egalite((Pqiinf(cas(il))), &
596	& 0.0,'compress_stewart 1047a')
597	call iso_verif_egalite(( &
598	& Eqi_prime(cas(il))),0.0,'compress_stewart 1047b')
599	call iso_verif_egalite(( &
600	& Pqiinf_par(cas(il))),0.0,'compress_stewart 1047c')
601	call iso_verif_egalite((Eqi_par(cas(il))), &
602	& 0.0,'compress_stewart 1047d')
603	#endif
604	Pqiinf_stewart(il)=0.0
605	Eqi_stewart(il)=0.0
606	endif ! if (Pqisup(cas(il).gt.0.0) then
607	enddo !do il=1,ncas_evap_glace
608
609	! petite vérif
610	#ifdef ISOVERIF
611	do il=1,ncas
612	if ((abs(water_cas(il)).ge.ridicule/10.) &
613	& .and.(Pqiinf_stewart(il).le.0.0)) then
614	write(,) 'compress_stewart 498: evap liq:'
615	write(,) 'water(il,i)=', water_cas(il)
616	write(,) 'Pqiinf=',Pqiinf(cas(il))
617	write(,) 'Pqiinf_par=',Pqiinf_par(cas(il))
618	write(,) 'Pqiinf_stewart=',Pqiinf_stewart(il)
619	stop
620	endif
621	enddo !do il=1,ncas_evap_glace
622	#endif
623
624	end subroutine compress_evap_liq_zone
625
626	subroutine compress_evap_glace_zone(iflag_con, &
627	& ncas, cas, &
628	& water_cas,water, &
629	& Pqisup_cas,Pqisup, &
630	& Pxtisup_cas,Pxtisup, &
631	& xtp_avantevap_cas,xtp_avantevap, &
632	& xtp_avantevaptrac_cas,qp_avantevaptrac_cas, &
633	& xtevapsup_cas,xtevap, &
634	& Eqi_stewart,Pqiinf_stewart,Eqi_prime_cas,Eqi_cas, &
635	& Pqiinf,Eqi_par,Pqiinf_par,Eqi_prime,Eqi, &
636	! & qp_cas,
637	#ifdef ISOVERIF
638	& evap_cas,evap, &
639	#endif
640	& nloc,ncum,nd,i,frac_sublim,izone)
641
642	USE isotopes_mod, ONLY: ridicule,iso_eau
643	#ifdef ISOVERIF
644	USE isotopes_verif_mod
645	#endif
646	implicit none
647
648	! compression dans le cas condensation_facftmr
649	integer iflag_con
650	integer nloc,nd,ncum
651	integer ncas
652	integer cas(ncum)
653	integer i
654	integer izone
655	real &
656	& water_cas(ncum),water(ncum), &
657	& xtevapsup_cas(niso,ncum),xtevap(ntraciso,ncum)
658	! real qp_cas(ncum)
659	#ifdef ISOVERIF
660	real evap_cas(ncum),evap(ncum)
661	#endif
662	real &
663	& Pqisup_cas(ncum),Pqisup(ncum), &
664	& Pxtisup_cas(niso,ncum),Pxtisup(ntraciso,ncum), &
665	& xtp_avantevap_cas(niso,ncum),xtp_avantevap(ntraciso,ncum), &
666	& Eqi_stewart(ncum),Pqiinf_stewart(ncum),Eqi_prime_cas(ncum), &
667	& Pqiinf(ncum),Eqi_par(ncum),Pqiinf_par(ncum),Eqi_prime(ncum), &
668	& Eqi(ncum),Eqi_cas(ncum)
669
670	real xtp_avantevaptrac_cas(niso,ncum), &
671	& qp_avantevaptrac_cas(ncum)
672	integer frac_sublim
673	! locals
674	integer il,ixt,ieau,iiso
675
676	ieau=index_trac(izone,iso_eau)
677	do il=1,ncas
678	Pqisup_cas(il)=Pxtisup(ieau,cas(il))
679
680	if (Pqisup(cas(il)).gt.0.0) then
681	Eqi_prime_cas(il)=Eqi_prime(cas(il)) &
682	& *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
683	Eqi_cas(il)=Eqi(il) &
684	& *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
685	else
686	#ifdef ISOVERIF
687	call iso_verif_egalite(( &
688	& Eqi_prime(cas(il))),0.0, &
689	& 'compress_stewart 979b')
690	#endif
691	Eqi_prime_cas(il)=0.0
692	Eqi_cas(il)=0.0
693	endif
694
695	if (Pqisup(cas(il))-Eqi_prime(cas(il)).gt.0.0) then
696	water_cas(il)=water(cas(il)) &
697	& *((Pxtisup(ieau,cas(il))-Eqi_prime_cas(il)) &
698	& /(Pqisup(cas(il))-Eqi_prime(cas(il))))
699	else !if (Pqisup(cas(il)).gt.0.0) then
700	#ifdef ISOVERIF
701	call iso_verif_egalite_choix(water(cas(il)),0.0, &
702	& 'compress_stewart 507',errmax,errmaxrel)
703	#endif
704	water_cas(il)=0.0
705	endif !if (Pqisup(cas(il)).gt.0.0) then
706
707
708	qp_avantevaptrac_cas(il)=xtp_avantevap(ieau,cas(il))
709	! qp_cas(il)=xtp(ieau,cas(il))
710	#ifdef ISOVERIF
711	! evap_cas(il)=evap(cas(il))
712	! & *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
713	! ce calcul est faux& faire plutot:
714	! evap_cas(il)=Eqi_prime_cas(il)//100.0/delP_cas(il)/SIGDg2.0
715	! & -(ieau,cas(il))
716	#endif
717	do iiso=1,niso
718	ixt=index_trac(izone,iiso)
719	Pxtisup_cas(iiso,il)=Pxtisup(ixt,cas(il))
720	xtp_avantevap_cas(iiso,il)=xtp_avantevap(iiso,cas(il))
721	xtp_avantevaptrac_cas(iiso,il)=xtp_avantevap(ixt,cas(il))
722	xtevapsup_cas(iiso,il)=xtevap(ixt,cas(il))
723	enddo
724	enddo !do il=1,ncas
725
726	! calculs des flux de masses à mettre en argument de stewart:
727	! comme l'eau n'est pas bien concervée dans les ddfts, on est
728	! obligé de bidouillé.
729	! 1) soit on considère Pqisup, Eqi, et Pqiinf_par=Pqisup-Eqi
730	! et on suppose que dans la réalité les compositions de
731	! Pqiinf sont les même que Pqiinf_par
732	! 2) soit on considère Pqisup, Eqi_par=Pqisup-Pqiinf, et Pqiinf,
733	! et on suppose que dans la réalité les compositions de
734	! Eqi_prime sont les même que Eqi_par
735	do il=1,ncas
736	if (Pqisup(cas(il)).gt.0.0) then
737
738	if ((water(cas(il)).gt.ridicule/100).and. &
739	& (Pqiinf_par(cas(il)).le.0.0)) then
740	! on ne peut pas utiliser la méthode 1, car KE prédit de l'eau
741	! alors que le bilan de masse n'enprédit pas.
742	! Peut-on utiliser la méthode 2?
743	Pqiinf_stewart(il)=Pqiinf(cas(il)) &
744	& *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
745	Eqi_stewart(il)=Eqi_par(cas(il)) &
746	& *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
747	else !if ((water(il,i).gt.ridicule/100).and.(Pqiinf_par.le.0.0)) then
748	! il n'y a pas d'obstacles à l'utilisation de 1)
749	Pqiinf_stewart(il)=Pqiinf_par(cas(il)) &
750	& *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
751	if (iflag_con.eq.30) then
752	Eqi_stewart(il)=Eqi_prime(cas(il)) &
753	& *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
754	else
755	Eqi_stewart(il)=Eqi(cas(il)) &
756	& *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
757	endif
758	endif !if ((water(il,i).gt.ridicule/100).and.(Pqiinf_par.le.0.0)) then
759
760	else !if (Pqisup(cas(il).gt.0.0) then
761
762	#ifdef ISOVERIF
763	call iso_verif_egalite((Pqiinf(cas(il))), &
764	& 0.0,'compress_stewart 1347a')
765	call iso_verif_egalite(( &
766	& Eqi_prime(cas(il))),0.0,'compress_stewart 1347b')
767	call iso_verif_egalite(( &
768	& Pqiinf_par(cas(il))),0.0,'compress_stewart 1347c')
769	call iso_verif_egalite((Eqi_par(cas(il))), &
770	& 0.0,'compress_stewart 1347d')
771	#endif
772	Pqiinf_stewart(il)=0.0
773	Eqi_stewart(il)=0.0
774	endif !if (Pqisup(cas(il).gt.0.0) then
775	enddo !do il=1,ncas_evap_glace
776
777	! petite vérif
778	#ifdef ISOVERIF
779	do il=1,ncas
780	if ((abs(water_cas(il)).ge.ridicule/10.) &
781	& .and.(Pqiinf_stewart(il).le.0.0)) then
782	write(,) 'compress_stewart 498: evap glace:'
783	write(,) 'water(il,i)=', water_cas(il)
784	write(,) 'Pqiinf=',Pqiinf(cas(il))
785	write(,) 'Pqiinf_par=',Pqiinf_par(cas(il))
786	write(,) 'Pqiinf_stewart=',Pqiinf_stewart(il)
787	stop
788	endif
789	enddo !do il=1,ncas_evap_glace
790	#endif
791
792	end subroutine compress_evap_glace_zone
793
794	subroutine uncompress_ilp_zone( &
795	& ncas,cas, &
796	& zxtrfln_cas,zxt_cas,zxtrfl,zxtrfln,zxt,klon, &
797	& izone,Eqi,Exi,fac_ftmr, &
798	& xtrevap_tag,liq,hdiag)
799
800	USE isotopes_mod, ONLY: ridicule,iso_eau
801	USE isotrac_mod, only: option_revap,evap_franche
802	#ifdef ISOVERIF
803	use isotopes_verif_mod
804	#endif
805
806	implicit none
807
808	! inputs
809	integer ncas
810	integer cas(ncas)
811	integer klon
812	integer izone
813	real zxt_cas(niso,ncas),zxtrfln_cas(niso,ncas)
814	real Exi(niso,ncas)
815	real Eqi(ncas)
816	real fac_ftmr(ncas)
817	integer liq
818	real hdiag(ncas)
819
820	! output
821	real zxt(ntraciso,klon)
822	real zxtrfl(ntraciso,klon),zxtrfln(ntraciso,klon)
823	real xtrevap_tag(ntraciso,ncas)
824
825	! locals
826	integer il,ixt,iiso
827
828	! write(,) 'uncompress_stewart 1420 tmp: zxt=',
829	! : zxt(iso_eau:ntraciso:3,cas(1))
830	! write(,) 'Exi,fac_ftmr=',
831	! : Exi(iso_eau,1),fac_ftmr(1)
832	do il=1,ncas
833	do iiso=1,niso
834	ixt=index_trac(izone,iiso)
835	zxtrfln(ixt,cas(il))=zxtrfln_cas(iiso,il)
836	zxtrfl(ixt,cas(il))=zxtrfln_cas(iiso,il)
837	enddo !do iiso=1,niso
838	!#ifdef ISOVERIF
839	! if (il.eq.9) then
840	! write(,) 'uncompress 1521'
841	! write(,) 'il,Eqi,fac_ftmr,evap_franche,Exi(2,il)=',
842	! : il,Eqi(il),fac_ftmr(il),evap_franche,Exi(2,il)
843	! endif
844	!#endif
845
846	if ((option_revap.eq.1).and. &
847	& (((liq.eq.1).and.(Eqi(il)*fac_ftmr(il).gt.evap_franche) &
848	& .and.(hdiag(il).lt.0.99)).or. &
849	& ((liq.eq.0).and. &
850	& (Eqi(il)*fac_ftmr(il).ge.0.0)))) then
851	! le flux d'évap va dans un tag particulier
852	! -> zxt est inchangé mais xtrevap_tag(ixt,il) est incrémenté
853	do iiso=1,niso
854	ixt=index_trac(izone,iiso)
855	xtrevap_tag(ixt,il)=fac_ftmr(il)*Exi(iiso,il)
856	! zxt(ixt,cas(il))=zxt_cas(iiso,il)
857	! : -xtrevap_tag(ixt,il)
858	enddo !do iiso=1,niso
859
860	else !if ((Eqi(il)*fac_ftmr(il).gt.evap_franche).and.
861	! reequilibration standard
862	do iiso=1,niso
863	ixt=index_trac(izone,iiso)
864	zxt(ixt,cas(il))=zxt(ixt,cas(il)) &
865	& +fac_ftmr(il)*Exi(iiso,il)
866	zxt(ixt,cas(il))=max(0.0,zxt(ixt,cas(il)))
867	xtrevap_tag(ixt,il)=0.0
868	#ifdef ISOVERIF
869	call iso_verif_positif_choix(zxt(ixt,cas(il)), &
870	& 0.0,'compress 1508')
871	#endif
872	enddo ! do iiso=1,niso
873	endif !if ((Eqi(il)*fac_ftmr(il).gt.evap_franche).and.
874	enddo !do il=1,ncas
875	! write(,) 'compress_stewart 1453 tmp: zxt=',
876	! : zxt(iso_eau:ntraciso:3,cas(1))
877
878	end subroutine uncompress_ilp_zone
879
880	subroutine compress_ilp_evap_liq_zone( &
881	& ncas,cas, &
882	& zxt_cas,zxt, &
883	& zxtrfl_cas,zxtrfl_ancien, &
884	& zrfln_cas,zrfln, &
885	& zrfl_cas,zrfl_ancien, &
886	& zqev_diag_cas,zqev_diag, &
887	& klon,izone,ptrac)
888
889	USE isotopes_mod, ONLY: ridicule,iso_eau
890	#ifdef ISOVERIF
891	USE isotopes_verif_mod
892	#endif
893	implicit none
894
895	! inputs
896	integer ncas
897	integer cas(ncas)
898	integer klon
899	real zxt(ntraciso,klon)
900	real zxtrfl(ntraciso,klon)
901	real zrfl_ancien(klon)
902	real zqev_diag(klon)
903	real zrfln(klon)
904	real zxtrfl_ancien(ntraciso,klon)
905	integer izone
906
907	! outputs
908	real zxt_cas(niso,ncas)
909	real zxtrfl_cas(niso,ncas)
910	real zqev_diag_cas(ncas)
911	real zrfln_cas(ncas)
912	real zrfl_cas(ncas)
913	real ptrac(ncas)
914
915	! locals
916	integer il,ixt,ieau,iiso
917
918	ieau=index_trac(izone,iso_eau)
919	do il=1,ncas
920	do iiso=1,niso
921	ixt=index_trac(izone,iiso)
922	! la compo de la vap à l'extérieure reste la vapeur totale
923	zxt_cas(iiso,il)=zxt(iiso,cas(il))
924	! le flux de pluie est celui le flux de pluie lié à la zone
925	zxtrfl_cas(iiso,il)=zxtrfl_ancien(ixt,cas(il))
926	enddo
927	zrfl_cas(il)=zxtrfl_ancien(ieau,cas(il))
928
929	if (zrfl_ancien(cas(il)).gt.0.0) then
930	! proportion de izone dans l'évap = celle dans la goutte
931	ptrac(il)=zxtrfl_ancien(ieau,cas(il))/zrfl_ancien(cas(il))
932	zrfln_cas(il)=zrfln(cas(il))*ptrac(il)
933	zqev_diag_cas(il)=zqev_diag(cas(il))*ptrac(il)
934	!#ifdef ISOVERIF
935	! if (il.eq.9) then
936	! write(,) 'compress tmp: il, ptrac=',il,ptrac(il)
937	! write(,) 'ieau,zxtrfl_ancien(ieau,cas(il))=',
938	! : ieau,zxtrfl_ancien(ieau,cas(il))
939	! write(,) 'zrfl_ancien(cas(il))=',zrfl_ancien(cas(il))
940	! write(,) 'zrfl_cas(il)=',zrfl_cas(il)
941	! endif
942	!#endif
943
944	else !if (zrfl_ancien(cas(il)).gt.0.0) then
945	#ifdef ISOVERIF
946	! write(,) 'il,cas(il),zrfln,zrfl_ancien,zqev_diag=',
947	! : il,cas(il),zrfln(cas(il)),zrfl_ancien(cas(il)),
948	! : zqev_diag(cas(il))
949	call iso_verif_egalite(zqev_diag(cas(il)),0.0, &
950	& 'compress_stewart 1591a')
951	call iso_verif_egalite(zrfln(cas(il)),0.0, &
952	& 'compress_stewart 1591b')
953	#endif
954	zrfln_cas(il)=0.0
955	zqev_diag_cas(il)=0.0
956
957	endif !if (zrfl_ancien(cas(il)).gt.0.0) then
958	! les lignes suvantes ne sont pas à recalculer
959	! zt_cas(il)=zt(cas(il))
960	! delP(il)=paprs(cas(il),k)-paprs(cas(il),k+1)
961	enddo !do il=1,ncas
962
963	end subroutine compress_ilp_evap_liq_zone
964
965
966	subroutine compress_ilp_evap_glace_zone( &
967	& ncas,cas, &
968	& zxt_cas,zxt, &
969	& zxtrfl_cas,zxtrfl_ancien, &
970	& zrfln_cas,zrfln, &
971	& zrfl_cas, zrfl_ancien, &
972	& zqev_diag_cas,zqev_diag, &
973	& klon,izone)
974
975	USE isotopes_mod, ONLY: ridicule,iso_eau
976	#ifdef ISOVERIF
977	use isotopes_verif_mod
978	#endif
979
980	implicit none
981
982	! inputs
983	integer ncas
984	integer cas(ncas)
985	integer klon
986	real zxt(ntraciso,klon)
987	real zxtrfl_ancien(ntraciso,klon)
988	real zqev_diag(klon)
989	real zrfln(klon)
990	integer izone
991	real zrfl_ancien(klon)
992
993	! outputs
994	real zxt_cas(niso,ncas)
995	real zxtrfl_cas(niso,ncas)
996	real zqev_diag_cas(ncas)
997	real zrfln_cas(ncas)
998	real zrfl_cas(ncas)
999
1000
1001	! locals
1002	integer il,ixt,ieau,iiso
1003
1004	ieau=index_trac(izone,iso_eau)
1005	do il=1,ncas
1006	do iiso=1,niso
1007	ixt=index_trac(izone,iiso)
1008	zxt_cas(iiso,il)=zxt(iiso,cas(il))
1009	zxtrfl_cas(iiso,il)=zxtrfl_ancien(ixt,cas(il))
1010	enddo
1011	zrfl_cas(il)=zxtrfl_ancien(ieau,cas(il))
1012
1013	if (zrfl_ancien(cas(il)).gt.0.0) then
1014	zrfln_cas(il)=zrfln(cas(il)) &
1015	& *(zxtrfl_ancien(ieau,cas(il))/zrfl_ancien(cas(il)))
1016	! car la proportion de traceurs dans zqev_diag et la même
1017	! que dans zrfl_ancien. Comme zrfln=zrfl-zqev_diag*fac_ftmr
1018	! alors cette proportion de traceurs est aussi la même dans
1019	! zrfln
1020	zqev_diag_cas(il)=zqev_diag(cas(il)) &
1021	& *zxtrfl_ancien(ieau,cas(il))/zrfl_ancien(cas(il))
1022	else !if (zrfl_ancien(cas(il)).gt.0.0) then
1023	#ifdef ISOVERIF
1024	call iso_verif_egalite(zqev_diag(cas(il)),0.0, &
1025	& 'compress_stewart 1791a')
1026	call iso_verif_egalite(zrfln(cas(il)),0.0, &
1027	& 'compress_stewart 1791b')
1028	#endif
1029	zrfln_cas(il)=0.0
1030	zqev_diag_cas(il)=0.0
1031	endif !if (zrfl_ancien(cas(il)).gt.0.0) then
1032	enddo
1033
1034	end subroutine compress_ilp_evap_glace_zone
1035
1036
1037	subroutine ajoute_revap(ncas,cas, &
1038	& klon,izone,zxt,xtrevap_tag)
1039
1040	#ifdef ISOVERIF
1041	USE isotopes_verif_mod
1042	#endif
1043	USE isotrac_mod, only: izone_revap
1044	implicit none
1045
1046	! ajoute xtrevap_tag (evaps des différents traceurs d'isotopes)
1047	! dans la vapeur qui est taggée par la révap des gouttes.
1048
1049	! input
1050	integer ncas
1051	integer cas(ncas)
1052	integer klon
1053	integer izone
1054	real xtrevap_tag(ntraciso,ncas)
1055
1056	! inout
1057	real zxt(ntraciso,klon)
1058
1059	! local
1060	integer i,ixt,iiso,ixt_revap
1061	!#ifdef ISOVERIF
1062	! integer iso_verif_positif_choix_nostop
1063	!#endif
1064
1065	#ifdef ISOVERIF
1066	do i=1,ncas
1067	do ixt=1,ntraciso
1068	call iso_verif_positif_choix(zxt(ixt,cas(i)),0.0, &
1069	& 'ajoute_revap 29')
1070	enddo
1071	enddo !do i=1,ncas
1072	#endif
1073
1074	do i=1,ncas
1075	do iiso=1,niso
1076	ixt_revap=index_trac(izone_revap,iiso)
1077	do izone=1,ntraceurs_zone
1078	ixt=index_trac(izone,iiso)
1079	zxt(ixt_revap,cas(i))=zxt(ixt_revap,cas(i)) &
1080	& +xtrevap_tag(ixt,i)
1081	#ifdef ISOVERIF
1082	if (iso_verif_positif_choix_nostop(zxt(ixt_revap,cas(i)), &
1083	& 0.0,'ajoute_revap 46').eq.1) then
1084	write(,) 'i,iiso,izone,ixt=',i,iiso,izone,ixt
1085	write(,) 'xtrevap_tag(ixt,i)=',xtrevap_tag(ixt,i)
1086	! stop
1087	endif
1088	#endif
1089	zxt(ixt_revap,cas(i))=max(0.0,zxt(ixt_revap,cas(i)))
1090	enddo !do izone=1,ntraceurs_zone
1091	enddo !do iiso=1,niso
1092	enddo !do i=1,ncas_evap_liq
1093
1094	#ifdef ISOVERIF
1095	do i=1,ncas
1096	do ixt=1,ntraciso
1097	call iso_verif_positif_choix(zxt(ixt,cas(i)),0.0, &
1098	& 'ajoute_revap 40')
1099	enddo
1100	enddo !do i=1,ncas
1101	#endif
1102
1103	return
1104	end subroutine ajoute_revap
1105
1106
1107	function is_in_bassin(lat,lon,bassin)
1108	USE isotrac_mod, only: use_bassin_atlantic,use_bassin_medit, &
1109	& use_bassin_indian,use_bassin_austral,use_bassin_pacific, &
1110	& use_bassin_MerArabie,use_bassin_BengalGolf,use_bassin_SouthIndian, &
1111	& use_bassin_tropics,use_bassin_midlats,use_bassin_HighLats, &
1112	& bassin_atlantic,bassin_medit, &
1113	& bassin_indian,bassin_austral,bassin_pacific, &
1114	& bassin_MerArabie,bassin_BengalGolf,bassin_SouthIndian, &
1115	& bassin_tropics,bassin_midlats,bassin_HighLats
1116	implicit none
1117	! répond true si lat,lon se trouve dans le bassin numéroté bassin
1118
1119	! input
1120	integer bassin
1121	real lat,lon
1122
1123	! outut
1124	logical is_in_bassin
1125
1126	! locals
1127	!logical is_in_rectangle
1128	!logical is_in_triangle
1129
1130	is_in_bassin=.false.
1131	#ifdef ISOVERIF
1132	write(,) 'is_in_basin 84: entree,bassin=',bassin
1133	#endif
1134	if (use_bassin_atlantic .and. bassin==bassin_atlantic) then
1135	#ifdef ISOVERIF
1136	write(,) 'bassin Atlantique?'
1137	#endif
1138	if (is_in_rectangle(lon,lat,-67.0,28.0,20.0,-45.0)) then
1139	! boite sud
1140	is_in_bassin=.true.
1141	return
1142	endif
1143	if (is_in_rectangle(lon,lat,-100.0,40.0,-5.3,28.0)) then
1144	! ouest gibraltar
1145	is_in_bassin=.true.
1146	return
1147	endif
1148	if (is_in_rectangle(lon,lat,-100.0,48.0,0.0,40.0)) then
1149	! Ouest France
1150	is_in_bassin=.true.
1151	return
1152	endif
1153	if (is_in_rectangle(lon,lat,-90.0,80.0,10.0,46.0)) then
1154	! Atlantic Nord
1155	is_in_bassin=.true.
1156	return
1157	endif
1158	if (is_in_triangle(lon,lat, &
1159	& -62.0,0.0,-62.0,30.0,-112.0,30.0)) then
1160	! golfe du Mexique
1161	is_in_bassin=.true.
1162	return
1163	endif
1164
1165	else if (use_bassin_medit .and. bassin==bassin_medit) then
1166	#ifdef ISOVERIF
1167	write(,) 'bassin Medit?'
1168	#endif
1169	if (is_in_rectangle(lon,lat,0.0,48.0,45.0,29.0)) then
1170	is_in_bassin=.true.
1171	return
1172	endif
1173	if (is_in_rectangle(lon,lat,-5.3,42.0,45.0,29.0)) then
1174	is_in_bassin=.true.
1175	return
1176	endif
1177
1178	else if (use_bassin_indian .and. bassin==bassin_indian) then
1179	#ifdef ISOVERIF
1180	write(,) 'bassin indian?'
1181	#endif
1182	if (is_in_rectangle(lon,lat,20.0,30.0,110.0,-45.0)) then
1183	is_in_bassin=.true.
1184	return
1185	endif
1186	if (is_in_triangle(lon,lat, &
1187	& 90.0,30.0,90.0,-45.0,150.0,-45.0)) then
1188	! Ouest Australie
1189	is_in_bassin=.true.
1190	return
1191	endif
1192
1193	else if (use_bassin_SouthIndian .and. bassin==bassin_SouthIndian) then
1194	#ifdef ISOVERIF
1195	write(,) 'bassin indian hemisphere Sud?'
1196	#endif
1197	if (is_in_rectangle(lon,lat,20.0,0.0,120.0,-45.0)) then
1198	is_in_bassin=.true.
1199	return
1200	endif
1201
1202	else if (use_bassin_MerArabie .and. bassin==bassin_MerArabie) then
1203	#ifdef ISOVERIF
1204	write(,) 'bassin Mer d''Arabie?'
1205	#endif
1206	if (is_in_rectangle(lon,lat,20.0,30.0,76.0,0.0)) then
1207	is_in_bassin=.true.
1208	return
1209	endif
1210
1211	else if (use_bassin_BengalGolf .and. bassin==bassin_BengalGolf) then
1212	#ifdef ISOVERIF
1213	write(,) 'bassin Golfe du Bengale?'
1214	#endif
1215	if (is_in_rectangle(lon,lat,76.0,30.0,110.0,0.0)) then
1216	is_in_bassin=.true.
1217	return
1218	endif
1219
1220	else if (use_bassin_pacific .and. bassin==bassin_pacific) then
1221	#ifdef ISOVERIF
1222	write(,) 'bassin Pacific?'
1223	#endif
1224	if (is_in_rectangle(lon,lat,-180.0,80.0,-100.0,-45.0)) then
1225	! pacifique Est
1226	is_in_bassin=.true.
1227	return
1228	endif
1229	if (is_in_rectangle(lon,lat,110.0,80.0,180.0,28.0)) then
1230	! Pacifique Nord Ouest
1231	is_in_bassin=.true.
1232	return
1233	endif
1234	if (is_in_rectangle(lon,lat,120.0,80.0,180.0,-45.0)) then
1235	! Pacifique central Sud
1236	is_in_bassin=.true.
1237	return
1238	endif
1239	if (is_in_triangle(lon,lat, &
1240	& 90.0,28.0,150.0,-45.0,150.0,28.0)) then
1241	! Pacifque Sud Ouest
1242	is_in_bassin=.true.
1243	return
1244	endif
1245	if (is_in_triangle(lon,lat, &
1246	& -62.0,0.0,-112.0,30.0,-112.0,0.0)) then
1247	! Ouest Amérique centrale
1248	is_in_bassin=.true.
1249	return
1250	endif
1251	if (is_in_rectangle(lon,lat,-180.0,0.0,-67.0,-45.0)) then
1252	! Ouest Chili
1253	is_in_bassin=.true.
1254	return
1255	endif
1256
1257	else if (use_bassin_austral .and. bassin==bassin_austral) then
1258	#ifdef ISOVERIF
1259	write(,) 'bassin austral?'
1260	#endif
1261	if (lat.lt.-45.0+0.2) then
1262	is_in_bassin=.true.
1263	return
1264	endif
1265
1266	else if (use_bassin_HighLats .and. bassin==bassin_HighLats) then
1267	#ifdef ISOVERIF
1268	write(,) 'bassin hautes lats?'
1269	#endif
1270	if (abs(lat).gt.35.0) then
1271	is_in_bassin=.true.
1272	return
1273	endif
1274
1275	else if (use_bassin_tropics .and. bassin==bassin_tropics) then
1276	#ifdef ISOVERIF
1277	write(,) 'bassin tropics?'
1278	#endif
1279	if (abs(lat).lt.15.0) then
1280	is_in_bassin=.true.
1281	return
1282	endif
1283
1284	else if (use_bassin_midlats .and. bassin==bassin_midlats) then
1285	#ifdef ISOVERIF
1286	write(,) 'bassin mid lats?'
1287	#endif
1288	if ((abs(lat).ge.15.0).and. &
1289	& (abs(lat).le.35.0)) then
1290	is_in_bassin=.true.
1291	return
1292	endif
1293
1294	else
1295	write(,) 'iso_traceurs_routines 59: bassin inconnu'
1296	write(,) 'bassin_atlantic=' ,bassin_atlantic
1297	write(,) 'bassin_medit=' ,bassin_medit
1298	write(,) 'bassin_indian=' ,bassin_indian
1299	write(,) 'bassin_austral=' ,bassin_austral
1300	write(,) 'bassin_MerArabie=' ,bassin_MerArabie
1301	write(,) 'bassin_BengalGolf=' ,bassin_BengalGolf
1302	write(,) 'bassin_SouthIndian=' ,bassin_SouthIndian
1303	write(,) 'use_bassin_atlantic=' ,use_bassin_atlantic
1304	write(,) 'use_bassin_medit=' ,use_bassin_medit
1305	write(,) 'use_bassin_indian=' ,use_bassin_indian
1306	write(,) 'use_bassin_austral=' ,use_bassin_austral
1307	write(,) 'use_bassin_MerArabie=' ,use_bassin_MerArabie
1308	write(,) 'use_bassin_BengalGolf=' ,use_bassin_BengalGolf
1309	write(,) 'use_bassin_SouthIndian=' ,use_bassin_SouthIndian
1310	stop
1311	endif
1312
1313	return
1314	end function is_in_bassin
1315
1316	subroutine find_bassin(lat,lon,bassin)
1317	use isotrac_mod, only: izone_poubelle,ntraceurs_zone=>ntiso,option_traceurs, &
1318	& bassin_map
1319	#ifdef ISOVERIF
1320	use isotopes_verif_mod
1321	#endif
1322	implicit none
1323
1324	! inputs
1325	real lat,lon
1326	! output
1327	integer bassin
1328	!locals
1329	logical continu
1330	!logical is_in_bassin
1331
1332	continu=.true.
1333	bassin=1
1334
1335	#ifdef ISOVERIF
1336	write(,) ''
1337	write(,) 'find bassin lat,lon=',lat,lon
1338	#endif
1339	do while (continu)
1340	!#ifdef ISOVERIF
1341	!! write(,) 'find_bassin 169: lat,lon,bassin=',lat,lon,bassin
1342	!#endif
1343	if (is_in_bassin(lat,lon,bassin)) then
1344	continu=.false.
1345	#ifdef ISOVERIF
1346	write(,) 'find_bassin 173: trouve: bassin=',bassin
1347	#endif
1348	else
1349	!#ifdef ISOVERIF
1350	! write(,) 'find_bassin 175: pas trouve: bassin=',bassin
1351	!#endif
1352	bassin=bassin+1
1353	endif
1354	if (bassin.eq.izone_poubelle) then
1355	continu=.false.
1356	bassin=izone_poubelle
1357	!#ifdef ISOVERIF
1358	! write(,) 'find_bassin 179: poubelle: bassin=',bassin
1359	!#endif
1360	endif
1361	enddo
1362
1363	! normalement, le bassin est soit un bassin oce, soit un résidu
1364	! donc bassin<=ntraceurs_zone-1
1365	#ifdef ISOVERIF
1366	call iso_verif_positif(float(ntraceurs_zone-1-bassin), &
1367	& 'find_bassin 195')
1368	#endif
1369
1370	return
1371	end subroutine find_bassin
1372
1373	subroutine initialise_bassins_boites(presnivs)
1374	USE dimphy, only: klev
1375	USE geometry_mod, ONLY : longitude_deg, latitude_deg
1376	use isotrac_mod, only: bassin_map,option_traceurs,boite_map
1377	#ifdef ISOVERIF
1378	use isotopes_verif_mod
1379	#endif
1380	implicit none
1381	real presnivs(klev)
1382
1383	if (option_traceurs.eq.3) then
1384	! initialisation de bassin_map
1385	call bassin_map_init(latitude_deg,longitude_deg,bassin_map)
1386	else if (option_traceurs.eq.20) then
1387	! initialisation de bassin_map selon < ou > 35° lat
1388	write(,) 'physiq 1681: init de la map pour tag 20'
1389	call bassin_map_init_opt20(latitude_deg,bassin_map)
1390	else if (option_traceurs.eq.5) then
1391	call boite_AMMA_init(latitude_deg,longitude_deg,presnivs,boite_map)
1392	else if (option_traceurs.eq.21) then
1393	call boite_UT_extra_init(latitude_deg,longitude_deg,presnivs,boite_map)
1394	endif
1395
1396	return
1397	end subroutine initialise_bassins_boites
1398
1399	subroutine bassin_map_init(lat,lon,bassin_map)
1400	USE dimphy, only: klon
1401	#ifdef ISOVERIF
1402	use isotopes_verif_mod
1403	#endif
1404
1405	implicit none
1406
1407	! inputs
1408	real lat(klon),lon(klon)
1409	! output
1410	integer bassin_map(klon)
1411	! locals
1412	integer i
1413
1414	do i=1,klon
1415	call find_bassin(lat(i),lon(i),bassin_map(i))
1416	#ifdef ISOVERIF
1417	write(,) 'init 233: i,lat,lon,bassin=',i,lat(i),lon(i), &
1418	& bassin_map(i)
1419	#endif
1420	enddo
1421
1422	return
1423	end subroutine bassin_map_init
1424
1425	function is_in_rectangle(x,y,x1,y1,x2,y2)
1426	implicit none
1427	! inputs
1428	real x,y
1429	! point en haut à gauche
1430	real x1,y1
1431	! point en bas à droite
1432	real x2,y2
1433	! output
1434	logical is_in_rectangle
1435
1436	!#ifdef ISOVERIF
1437	! write(,) 'is_in_rectange 237: x,y=',x,y
1438	! write(,) 'x1,y1,x2,y2=',x1,y1,x2,y2
1439	!#endif
1440	if ((x-x2.lt.0.1).and.(x-x1.gt.-0.1).and. &
1441	& (y-y1.lt.0.1).and.(y-y2.gt.-0.1)) then
1442	is_in_rectangle=.true.
1443	else
1444	is_in_rectangle=.false.
1445	endif
1446	!#ifdef ISOVERIF
1447	! write(,) 'is_in_rectangle=',is_in_rectangle
1448	!#endif
1449	return
1450	end function is_in_rectangle
1451
1452	function is_in_triangle(x,y,x1,y1,x2,y2,x3,y3)
1453	implicit none
1454	! inputs
1455	real x,y
1456	! points dans le sens trigo
1457	! à gauche
1458	real x1,y1
1459	! en bas
1460	real x2,y2
1461	! à droite
1462	real x3,y3
1463	! output
1464	logical is_in_triangle
1465	! locals
1466	real det1
1467	real det2
1468	real det3
1469	!#ifdef ISOVERIF
1470	! write(,) 'is_in_triange 271: x,y=',x,y
1471	! write(,) 'x1,y1,x2,y2,x3,y3=',x1,y1,x2,y2,x3,y3
1472	!#endif
1473	det1=(x1-x)(y2-y)-(y1-y)(x2-x)
1474	det2=(x2-x)(y3-y)-(y2-y)(x3-x)
1475	det3=(x3-x)(y1-y)-(y3-y)(x1-x)
1476	is_in_triangle=.false.
1477	if ((det1det2.gt.0.0).and.(det2det3.gt.0.0)) then
1478	is_in_triangle=.true.
1479	else
1480	is_in_triangle=.false.
1481	endif
1482	!#ifdef ISOVERIF
1483	! write(,) 'det1,det2,det3,is_in_triangle',
1484	! : det1,det2,det3,is_in_triangle
1485	!#endif
1486	return
1487	end function is_in_triangle
1488
1489
1490	subroutine isotrac_recolorise_tmin(xt,t)
1491	USE dimphy, only: klon, klev
1492	USE isotrac_mod, only: zone_temp,nzone_temp
1493	#ifdef ISOVERIF
1494	USE isotopes_verif_mod
1495	#endif
1496	implicit none
1497
1498	! inout
1499	real xt(ntraciso,klon,klev)
1500	! input
1501	real t(klon,klev)
1502	! locals
1503	integer izone_temp
1504	integer ixt,ixt_recoit
1505	integer k,i,izone,iiso
1506	!integer find_index
1507
1508
1509	do k=1,klev
1510	do i=1,klon
1511	!#ifdef ISOVERIF
1512	! if (i.eq.1) then
1513	! write(,) 'recolorise 396: i,k,t=',i,k,t(i,k)
1514	! write(,) 'zone_temp=',zone_temp
1515	! endif
1516	!#endif
1517	! trouver la zone de cette température
1518	izone_temp=find_index(t(i,k),nzone_temp,zone_temp)
1519
1520	!#ifdef ISOVERIF
1521	! if (i.eq.1) then
1522	! write(,) 'recolorise 414: izone_temp=',izone_temp
1523	! endif
1524	!#endif
1525	do izone=1,nzone_temp-1
1526	! tous les tags de zone < nzone_temp se trouvant à des
1527	! températures plus basses sont convertis
1528	!#ifdef ISOVERIF
1529	! if (i.eq.1) then
1530	! write(,) 'recolorise 405: izone,xt_eau=',
1531	! : izone,xt(index_trac(izone,iso_eau),i,k)
1532	! endif
1533	!#endif
1534	if (izone.lt.izone_temp) then
1535	do iiso=1,niso
1536	ixt=index_trac(izone,iiso) ! emmetteur
1537	ixt_recoit=index_trac(izone_temp,iiso) ! recepteur
1538	xt(ixt_recoit,i,k)=xt(ixt_recoit,i,k)+xt(ixt,i,k)
1539	xt(ixt,i,k)=0.0
1540	enddo !do iiso=1,niso
1541	endif !if (izone.lt.izone_temp) then
1542	!#ifdef ISOVERIF
1543	! if (i.eq.1) then
1544	! write(,) 'recolorise 419: xt_eau,xt_recoit=',
1545	! : xt(index_trac(izone,iso_eau),i,k),
1546	! : xt(index_trac(izone_temp,iso_eau),i,k)
1547	! endif
1548	!#endif
1549	enddo !do izone=1,zone_pot(k)-1
1550
1551	! conversion de l'évap en surf et de la revap des gouttes
1552	do izone=nzone_temp+1,ntraceurs_zone
1553	do iiso=1,niso
1554	ixt=index_trac(izone,iiso) ! emmetteur
1555	ixt_recoit=index_trac(izone_temp,iiso) ! recepteur
1556	xt(ixt_recoit,i,k)=xt(ixt_recoit,i,k)+xt(ixt,i,k)
1557	xt(ixt,i,k)=0.0
1558	enddo !do iiso=1,niso
1559	enddo !do izone=nzone_temp+1,ntraceurs_zone
1560	enddo !do i=1,klon
1561	enddo !do k=1,klev
1562
1563	#ifdef ISOVERIF
1564	do k=1,klev
1565	do i=1,klon
1566	call iso_verif_traceur(xt(1,i,k),'recolorise 403')
1567	enddo !do i=1,klon
1568	enddo !do k=1,klev
1569	#endif
1570
1571	return
1572	end subroutine isotrac_recolorise_tmin
1573
1574	subroutine isotrac_recolorise_tmin_sfrev(xt,t)
1575	USE dimphy, only: klon,klev
1576	USE isotrac_mod, only: nzone_temp,zone_temp
1577	#ifdef ISOVERIF
1578	USE isotopes_verif_mod
1579	#endif
1580	implicit none
1581
1582	! recolorise selon la température minimum, mais les tags de
1583	! revap sont laissés en revap
1584
1585	! inout
1586	real xt(ntraciso,klon,klev)
1587	! input
1588	real t(klon,klev)
1589	! locals
1590	integer izone_temp
1591	integer ixt,ixt_recoit
1592	integer k,i,izone,iiso
1593	!integer find_index
1594
1595	do k=1,klev
1596	do i=1,klon
1597	! trouver la zone de cette température
1598	izone_temp=find_index(t(i,k),nzone_temp,zone_temp)
1599
1600	do izone=1,nzone_temp-1
1601	! tous les tags de zone < nzone_temp se trouvant à des
1602	! températures plus basses sont convertis
1603	! sauf la revap
1604	! le tag de la revap est nzone_temp+1=ntraceurs_zone
1605	if (izone.lt.izone_temp) then
1606	do iiso=1,niso
1607	ixt=index_trac(izone,iiso) ! emmetteur
1608	ixt_recoit=index_trac(izone_temp,iiso) ! recepteur
1609	xt(ixt_recoit,i,k)=xt(ixt_recoit,i,k)+xt(ixt,i,k)
1610	xt(ixt,i,k)=0.0
1611	enddo !do iiso=1,niso
1612	endif !if (izone.lt.izone_temp) then
1613	enddo !do izone=1,zone_pot(k)-1
1614
1615	enddo !do i=1,klon
1616	enddo !do k=1,klev
1617
1618	#ifdef ISOVERIF
1619	do k=1,klev
1620	do i=1,klon
1621	call iso_verif_traceur(xt(1,i,k),'recolorise 594')
1622	enddo !do i=1,klon
1623	enddo !do k=1,klev
1624	#endif
1625
1626	return
1627	end subroutine isotrac_recolorise_tmin_sfrev
1628
1629
1630	subroutine isotrac_recolorise_saturation(xt,rh,lat,pres)
1631	USE dimphy, only: klon,klev
1632	#ifdef ISOVERIF
1633	USE isotopes_verif_mod
1634	#endif
1635	implicit none
1636
1637	! recolorise selon la température minimum, mais les tags de
1638	! revap sont laissés en revap
1639
1640	! inout
1641	real xt(ntraciso,klon,klev)
1642	! input
1643	real rh(klon,klev)
1644	real lat(klon)
1645	real pres(klev)
1646	! locals
1647	integer izone_recoit
1648	integer ixt,ixt_recoit
1649	integer k,i,izone,iiso
1650	logical continu
1651	real rh_seuil
1652	parameter (rh_seuil=0.90)
1653	!integer index_zone_latpres
1654
1655	#ifdef ISOVERIF
1656	do k=1,klev
1657	do i=1,klon
1658	call iso_verif_traceur(xt(1,i,k),'recolorise 612')
1659	enddo !do i=1,klon
1660	enddo !do k=1,klev
1661	#endif
1662
1663	! on ne sature pas les 2 premières couches: on les laisse se
1664	! recharger en evap de surface
1665	do k=3,klev
1666	do i=1,klon
1667	if (rh(i,k).gt.rh_seuil) then
1668	izone_recoit=index_zone_latpres(lat(i),pres(k))
1669	do izone=1,ntraceurs_zone
1670	if (izone.ne.izone_recoit) then
1671	do iiso=1,niso
1672	ixt=index_trac(izone,iiso) ! emmetteur
1673	ixt_recoit=index_trac(izone_recoit,iiso) ! recepteur
1674	xt(ixt_recoit,i,k)=xt(ixt_recoit,i,k)+xt(ixt,i,k)
1675	xt(ixt,i,k)=0.0
1676	enddo !do iiso=1,niso
1677	endif
1678	enddo !do izone=1,ntraceurs_zone
1679	endif !if (rh(i,k).gt.rh_seuil) then
1680	enddo !do i=1,klon
1681	enddo !do k=1,klev
1682
1683	#ifdef ISOVERIF
1684	do k=1,klev
1685	do i=1,klon
1686	call iso_verif_traceur(xt(1,i,k),'recolorise 637')
1687	enddo !do i=1,klon
1688	enddo !do k=1,klev
1689	#endif
1690
1691	return
1692	end subroutine isotrac_recolorise_saturation
1693
1694	subroutine isotrac_recolorise_boite(xt,boite_map)
1695	USE dimphy, only: klon,klev
1696	#ifdef ISOVERIF
1697	USE isotopes_verif_mod
1698	#endif
1699	implicit none
1700
1701	! subroutine écrite à la base pour tagguer 3 boites AMMA.
1702	! Mais ça peut être générique, selon comment est initialisée boite_map
1703
1704	! inout
1705	real xt(ntraciso,klon,klev)
1706	! input
1707	integer boite_map(klon,klev)
1708	! locals
1709	integer i,k
1710	integer izone_recoit,izone,iiso
1711	integer ixt,ixt_recoit
1712
1713	do k=1,klev
1714	do i=1,klon
1715	izone_recoit=boite_map(i,k)
1716	if (izone_recoit.gt.0) then
1717	! on est dans une boite connue
1718	! toutes les molécules sont converties en cete couleur
1719	do izone=1,ntraceurs_zone
1720	if (izone.ne.izone_recoit) then
1721	! on met les traceurs izone dans izone_recoit
1722	do iiso=1,niso
1723	ixt=index_trac(izone,iiso)
1724	ixt_recoit=index_trac(izone_recoit,iiso)
1725	xt(ixt_recoit,i,k)=xt(ixt_recoit,i,k) &
1726	& +xt(ixt,i,k)
1727	xt(ixt,i,k)=0.0
1728	enddo
1729	endif !if (izone.ne.izone_recoit) then
1730	enddo !do izone=2,ntraceurs_zone
1731	endif !if (izone_recoit.gt.0) then
1732	enddo !do i=1,klon
1733	enddo !do k=1,klev
1734
1735	#ifdef ISOVERIF
1736	do k=1,klev
1737	do i=1,klon
1738	call iso_verif_traceur(xt(1,i,k),'recolorise 514')
1739	enddo !do i=1,klon
1740	enddo !do k=1,klev
1741	#endif
1742
1743	return
1744	end subroutine isotrac_recolorise_boite
1745
1746	subroutine isotrac_recolorise_extra(xt,rlat)
1747	USE dimphy, only: klon,klev
1748	usE isotrac_mod, only: lim_tag20,izone_trop,izone_extra
1749	#ifdef ISOVERIF
1750	USE isotopes_verif_mod
1751	#endif
1752	implicit none
1753
1754	! subroutine écrite pour l'option de taggage 20
1755	! permet de retagguer la vapeur tropicale en vapeur
1756	! extratropicale dès qu'elle atteint 35° de latitude
1757
1758	! inout
1759	real xt(ntraciso,klon,klev)
1760	! input
1761	real rlat(klon)
1762	! locals
1763	integer i,k
1764	integer iiso,ixt,ixt_recoit
1765
1766	! write(,) 'iso_traceurs_routines 723: lim_tag20=',lim_tag20
1767	do k=1,klev
1768	do i=1,klon
1769	if (abs(rlat(i)).gt.lim_tag20) then
1770	! on met les traceurs izone_trop dans izone_extra
1771	do iiso=1,niso
1772	ixt=index_trac(izone_trop,iiso)
1773	ixt_recoit=index_trac(izone_extra,iiso)
1774	xt(ixt_recoit,i,k)=xt(ixt_recoit,i,k) &
1775	& +xt(ixt,i,k)
1776	xt(ixt,i,k)=0.0
1777	enddo
1778	endif ! if (abs(rlat(i)).lt.35.0) then
1779	enddo !do i=1,klon
1780	enddo !do k=1,klev
1781
1782	#ifdef ISOVERIF
1783	do k=1,klev
1784	do i=1,klon
1785	call iso_verif_traceur(xt(1,i,k),'recolorise 741')
1786	enddo !do i=1,klon
1787	enddo !do k=1,klev
1788	#endif
1789
1790	return
1791	end subroutine isotrac_recolorise_extra
1792
1793	subroutine isotrac_recolorise_conv(xt,rlat,presnivs,rain_con)
1794	USE dimphy, only: klon,klev
1795	use isotrac_mod, only: lim_precip_tag22, &
1796	& izone_conv_BT,izone_conv_UT
1797	#ifdef ISOVERIF
1798	USE isotopes_verif_mod
1799	#endif
1800	implicit none
1801
1802	! subroutine écrite pour l'option de taggage 20
1803	! permet de retagguer la vapeur tropicale en vapeur
1804	! extratropicale dès qu'elle atteint 35° de latitude
1805
1806	! inout
1807	real xt(ntraciso,klon,klev)
1808	! input
1809	real rlat(klon)
1810	real presnivs(klev)
1811	real rain_con(klon)
1812	! locals
1813	integer i,k
1814	integer iiso,ixt,ixt_recoit,izone
1815
1816	! write(,) 'iso_traceurs_routines 723: lim_tag20=',lim_tag20
1817	! write(,) 'presnivs=',presnivs
1818	! stop
1819	do k=1,klev
1820	do i=1,klon
1821	#ifdef ISOVERIF
1822	if ((abs(rlat(i)).lt.30.0).and.(k.eq.1)) then
1823	endif
1824	#endif
1825	if ((abs(rlat(i)).lt.30.0).and. &
1826	& (rain_con(i)*86400.gt.lim_precip_tag22)) then
1827
1828	! on met les traceurs izone_trop dans izone_conv fn z
1829	do iiso=1,niso
1830	if (presnivs(k).gt.650.0*100.0) then
1831	ixt_recoit=index_trac(izone_conv_BT,iiso)
1832	else
1833	ixt_recoit=index_trac(izone_conv_UT,iiso)
1834	endif
1835	do izone=1,ntraceurs_zone
1836	ixt=index_trac(izone,iiso)
1837	if (ixt.ne.ixt_recoit) then
1838	xt(ixt_recoit,i,k)=xt(ixt_recoit,i,k) &
1839	& +xt(ixt,i,k)
1840	xt(ixt,i,k)=0.0
1841	endif !if (ixt.ne.ixt_recoit) then
1842	enddo !do izone=1,ntraceurs_zone
1843	enddo !do iiso=1,niso
1844	! write(,) 'k,presnivs,ixt,ixt_recoit=',k,presnivs(k),
1845	! : ixt,ixt_recoit
1846	! write(,) 'xt(:,i,k)=',xt(:,i,k)
1847	endif ! if (abs(rlat(i)).lt.35.0) then
1848	enddo !do i=1,klon
1849	enddo !do k=1,klev
1850
1851	#ifdef ISOVERIF
1852	do k=1,klev
1853	do i=1,klon
1854	call iso_verif_traceur(xt(1,i,k),'recolorise 741')
1855	enddo !do i=1,klon
1856	enddo !do k=1,klev
1857	#endif
1858
1859	return
1860	end subroutine isotrac_recolorise_conv
1861
1862
1863	subroutine boite_AMMA_init(lat,lon,presnivs,boite_map)
1864	USE dimphy, only: klon,klev
1865	#ifdef ISOVERIF
1866	USE isotopes_verif_mod
1867	#endif
1868	USE isotrac_mod, only: izone_aej,izone_mousson,izone_harmattan
1869	implicit none
1870
1871	real lat(klon),lon(klon)
1872	real presnivs(klev)
1873	! output
1874	integer boite_map(klon,klev)
1875	! locals
1876	integer i,k
1877
1878	! write(,) 'izone_aej,izone_mousson,izone_harmattan=',
1879	! : izone_aej,izone_mousson,izone_harmattan
1880	do k=1,klev
1881	do i=1,klon
1882	boite_map(i,k)=0.0
1883	! write(,) 'i,k,lat,lon,pres=',
1884	! : i,k,lat(i),lon(i),presnivs(k)
1885	if ((presnivs(k).le.700.0*100.0).and. &
1886	& (presnivs(k).gt.400.0*100.0).and. &
1887	& (lat(i).gt.8.0).and. &
1888	& (lat(i).lt.20.0).and. &
1889	& (lon(i).gt.10.0).and. &
1890	& (lon(i).lt.30.0)) then
1891	boite_map(i,k)=izone_aej
1892	! write(,) ' -> zone AEJ'
1893	else if ((presnivs(k).ge.850.0*100.0).and. &
1894	& (lat(i).gt.-5.0).and. &
1895	& (lat(i).le.8.0).and. &
1896	& (lon(i).gt.-40.0).and. &
1897	& (lon(i).lt.15.0)) then
1898	boite_map(i,k)=izone_mousson
1899	! write(,) ' -> zone flux de mousson'
1900	else if ((presnivs(k).gt.700.0*100.0).and. &
1901	& (lat(i).ge.20.0).and. &
1902	& (lat(i).lt.30.0).and. &
1903	& (lon(i).gt.-10.0).and. &
1904	& (lon(i).lt.40.0)) then
1905	boite_map(i,k)=izone_harmattan
1906	! write(,) ' -> zone Harmattan'
1907	endif
1908	! write(,) ' ** boite_map=',boite_map(i,k)
1909	enddo
1910	enddo
1911
1912	return
1913	end subroutine boite_AMMA_init
1914
1915
1916	subroutine boite_UT_extra_init(lat,lon,presnivs,boite_map)
1917	USE dimphy, only: klon,klev
1918	use isotrac_mod, only: izone_extra,izone_trop
1919	#ifdef ISOVERIF
1920	USE isotopes_verif_mod
1921	#endif
1922	implicit none
1923
1924	real lat(klon),lon(klon)
1925	real presnivs(klev)
1926	! output
1927	integer boite_map(klon,klev)
1928	! locals
1929	integer i,k
1930
1931	! write(,) 'izone_trop,izone_extra=',
1932	! : izone_trop,izone_extra
1933	do k=1,klev
1934	do i=1,klon
1935	boite_map(i,k)=0.0
1936	! write(,) 'i,k,lat,lon,pres=',
1937	! : i,k,lat(i),lon(i),presnivs(k)
1938	if ((presnivs(k).le.500.0*100.0) &
1939	& .and.(abs(lat(i)).lt.15.0)) then
1940	boite_map(i,k)=izone_trop
1941	! write(,) ' -> zone trop'
1942	else if (abs(lat(i)).gt.35.0) then
1943	boite_map(i,k)=izone_extra
1944	! write(,) ' -> zone extratropiques'
1945	endif
1946	! write(,) ' ** boite_map=',boite_map(i,k)
1947	enddo
1948	enddo
1949
1950	return
1951	end subroutine boite_UT_extra_init
1952
1953
1954	function index_zone_lat(lat)
1955	use isotrac_mod, only: lattag_min,dlattag,nzone_lat
1956	implicit none
1957
1958	! inputs
1959	real lat,pres
1960	! output
1961	integer index_zone_lat
1962
1963	if (lat.lt.lattag_min) then
1964	index_zone_lat=1
1965	else
1966	index_zone_lat=int((lat-lattag_min)/dlattag)+2
1967	index_zone_lat=min(index_zone_lat,nzone_lat)
1968	endif
1969
1970	return
1971	end function index_zone_lat
1972
1973
1974	function index_zone_pres(pres)
1975	use isotrac_mod, only: nzone_pres,zone_pres
1976	implicit none
1977
1978	! inputs
1979	real lat,pres
1980	! output
1981	integer index_zone_pres
1982	!integer find_index
1983
1984
1985	index_zone_pres=find_index(pres,nzone_pres,zone_pres)
1986	write(,) 'iso_traceurs_routines 802: pres,index_zone_pres=', &
1987	& pres,index_zone_pres
1988	write(,) 'zone_pres=',zone_pres(1:nzone_pres-1)
1989
1990	return
1991	end function index_zone_pres
1992
1993	function find_index(pres,nzone_pres,zone_pres)
1994	implicit none
1995
1996	! inputs
1997	real pres
1998	integer nzone_pres
1999	real zone_pres(nzone_pres)
2000	! output
2001	integer find_index
2002	logical continu
2003
2004	if (nzone_pres.gt.1) then
2005	if (pres.ge.zone_pres(1)) then
2006	find_index=1
2007	else if (pres.lt.zone_pres(nzone_pres-1)) then
2008	find_index=nzone_pres
2009	else !if (t(i,k).ge.zone_temp1) then
2010	continu=.true.
2011	find_index=2
2012	do while (continu)
2013	if (pres.ge.zone_pres(find_index)) then
2014	continu=.false.
2015	! c'est izone_temp, zone trouvée
2016	else
2017	find_index=find_index+1
2018	endif
2019	enddo !do while (continu)
2020	endif !if (t(i,k).ge.zone_temp1) then
2021
2022	else !if (nzone_pres.gt.1) then
2023	find_index=1
2024	endif !if (nzone_pres.gt.1) then
2025
2026	return
2027	end function find_index
2028
2029	function index_zone_latpres(lat,pres)
2030	use isotrac_mod, only: nzone_lat
2031	implicit none
2032
2033	! inputs
2034	real lat,pres
2035	! output
2036	integer index_zone_latpres
2037	! locals
2038	integer index_lat
2039	integer index_pres
2040	!integer index_zone_lat
2041	!integer index_zone_pres
2042
2043	index_lat=index_zone_lat(lat)
2044	index_pres=index_zone_pres(pres)
2045	index_zone_latpres=index_lat+(index_pres-1)*nzone_lat
2046
2047	return
2048	end function index_zone_latpres
2049
2050	subroutine iso_recolorise_condensation(qt,cond, &
2051	& xt,zxtcond,tcond,ep,xtres, &
2052	& seuil_in)
2053	USE dimphy, only: klon,klev
2054	USE isotopes_mod, only: bidouille_anti_divergence,iso_eau
2055	use isotrac_mod, only: option_seuil_tag_tmin,izone_cond, &
2056	& nzone_temp,zone_temp
2057	#ifdef ISOVERIF
2058	USE isotopes_verif_mod
2059	#endif
2060	implicit none
2061
2062	! on recolorise la vapeur résiduelle selon la température de condensation
2063	! on supose qu'une vapeur xt,q condense en cond,zxtcond, à une
2064	! température tcond. A ce stade, la vapeur initiale n'est pas
2065	! retranchée de son condensat. On calcule les tags dans la vepru
2066	! résiduelle xtres qu'on aurait si on retranchait un fraction ep du
2067	! condensat
2068
2069	! inputs
2070	real qt
2071	real cond
2072	real tcond
2073	real ep
2074	real xt(ntraciso)
2075	real zxtcond(ntraciso)
2076	real seuil_in
2077	! outputs
2078	real xtres(ntraciso)
2079	! locals
2080	integer izone_temp,izone
2081	integer ixt,ixt_recoit
2082	integer iiso
2083	!integer find_index
2084	real fcond, qmicro
2085	! real f
2086
2087	if ((cond.gt.0.0).and.(qt.gt.0.0)) then
2088
2089	izone_temp=find_index(tcond,nzone_temp,zone_temp)
2090	! WRITE(,) 'pgm 901 tmp: izone_temp=',izone_temp
2091	#ifdef ISOVERIF
2092	do ixt=1,ntraciso
2093	call iso_verif_positif(xt(ixt)-zxtcond(ixt), &
2094	& 'iso_trac 898')
2095	enddo !do ixt=1,ntraciso
2096	call iso_verif_traceur_justmass(xt, &
2097	& 'iso_trac_routines 906')
2098	call iso_verif_traceur_justmass(zxtcond, &
2099	& 'iso_trac_routines 908')
2100	#endif
2101	! bidouille
2102	if (bidouille_anti_divergence) then
2103	call iso_verif_traceur_jbidouille(xt)
2104	call iso_verif_traceur_jbidouille(zxtcond)
2105	endif
2106
2107	do ixt=1,niso
2108	xtres(ixt)=xt(ixt)-ep*zxtcond(ixt)
2109	enddo
2110	do ixt=1+niso,ntraciso
2111	xtres(ixt)=0.0
2112	enddo
2113	! write(,) 'iso_trac_routines tmp 916: xtres=',xtres
2114
2115	#ifdef ISOVERIF
2116	do ixt=1,ntraciso
2117	call iso_verif_positif(xtres(ixt), &
2118	& 'iso_trac_routines 921')
2119	enddo
2120	#endif
2121
2122	! cas de izone sfc et izone precip et izone cond et izone< izone_temp
2123
2124	! write(,) 'iso_trac 940: cond/qt,seuil_in,izone_temp=',
2125	! : cond/qt,seuil_in,izone_temp
2126
2127	if (option_seuil_tag_tmin.eq.2) then
2128	qmicro=0.0
2129	do izone=nzone_temp+1,ntraceurs_zone
2130	ixt= index_trac(izone,iso_eau)
2131	qmicro=qmicro+xt(ixt)
2132	enddo !do izone=nzone_temp+1,ntraceurs_zone
2133	if (qt-qmicro.gt.0.0) then
2134	fcond=(cond-qmicro)/(qt-qmicro)
2135	else
2136	fcond=0.0
2137	endif
2138	else
2139	fcond=cond/qt
2140	endif
2141
2142	if (fcond.gt.seuil_in) then
2143
2144	! on les transfert à izone_temp
2145	do izone=1,ntraceurs_zone
2146	if ((izone.gt.nzone_temp).or.(izone.lt.izone_temp)) then
2147	! ieau=index_trac(izone,iso_eau)
2148	do iiso=1,niso
2149	ixt= index_trac(izone,iiso)
2150	ixt_recoit=index_trac(izone_temp,iiso) ! recepteur
2151	xtres(ixt_recoit)=xtres(ixt_recoit) &
2152	& +(xt(ixt)-zxtcond(ixt))
2153	xtres(ixt)=0.0
2154	! write(,) 'iso_trac 920: izone,ixt,',
2155	! : 'ixt_recoit=',
2156	! : izone,ixt,ixt_recoit
2157	! write(,) 'isotrac 924: xt=',xt
2158	! write(,) 'isotrac 925: zxtcond=',zxtcond
2159	enddo !do iiso=1,niso
2160	! write(,) 'iso_trac tmp 944: izone,xtres=',izone,xtres
2161	endif !if (izone.ne.izone_cond) then
2162	enddo !do izone=nzones_temp+1,ntraceurs_zone
2163
2164	else !if (cond/qt.gt.seuil_in) then
2165
2166	! on les laisse sur place
2167	do izone=1,ntraceurs_zone
2168	if ((izone.gt.nzone_temp).or.(izone.lt.izone_temp)) then
2169	do iiso=1,niso
2170	ixt= index_trac(izone,iiso)
2171	xtres(ixt)=(xt(ixt)-zxtcond(ixt))
2172	enddo !do iiso=1,niso
2173	endif !if (izone.ne.izone_cond) then
2174	enddo !do izone=nzones_temp+1,ntraceurs_zone
2175
2176	endif !if (cond/qt.gt.seuil_in) then
2177
2178	! izone_temp est conservé, on lui enlève juste son
2179	! condesat
2180	do iiso=1,niso
2181	ixt_recoit=index_trac(izone_temp,iiso) ! recepteur
2182	xtres(ixt_recoit)=xtres(ixt_recoit) &
2183	& +(xt(ixt_recoit)-zxtcond(ixt_recoit))
2184	enddo !do iiso=1,niso
2185
2186	#ifdef ISOVERIF
2187	do ixt=1,ntraciso
2188	call iso_verif_positif(xtres(ixt), &
2189	& 'iso_trac_routines 940')
2190	enddo
2191	#endif
2192
2193	! cas des zones > izone temp
2194	! on conserve le condensat résiduel
2195	do izone=izone_temp+1,nzone_temp
2196	do iiso=1,niso
2197	ixt= index_trac(izone,iiso)
2198	xtres(ixt)=xt(ixt)-zxtcond(ixt)
2199	! write(,) 'iso_trac 931: izone,ixt,ixt_recoit=',
2200	! : izone,ixt,ixt_recoit
2201	! write(,) 'isotrac 934: xt=',xt
2202	! write(,) 'isotrac 935: zxtice=',zxtice
2203	enddo !do iiso=1,niso
2204	! write(,) 'iso_trac tmp 965: izone,xtres=',izone,xtres
2205	enddo !do izone=izone_temp+1,nzones_temp
2206
2207	! on rajoute le condensat qui ne precipite pas
2208	if (ep.lt.1.0) then
2209	do iiso=1,niso
2210	ixt= index_trac(izone_cond,iiso)
2211	xtres(ixt)=xtres(ixt)+(1.0-ep)*zxtcond(iiso)
2212	! write(,) 'iso_trac 940: izone,ixt,ixt_recoit=',
2213	! : izone,ixt,ixt_recoit
2214	! write(,) 'isotrac 1014: xt=',xt
2215	! write(,) 'isotrac 945: zxtice=',zxtice
2216	enddo !do iiso=1,niso
2217	endif !if (ep.lt.0.0) then
2218
2219	else
2220	! si cond=0 ou qt=0, tot reste pareil
2221	do ixt=1,ntraciso
2222	xtres(ixt)=xt(ixt)
2223	enddo !do ixt=1,ntraciso
2224
2225	endif ! if (qt.gt.0.0) then
2226
2227	#ifdef ISOVERIF
2228	if (iso_verif_traceur_jm_nostop(xtres, &
2229	& 'iso_trac_routines 166').eq.1) then
2230	write(,) 'isotrac 1024: xt=',xt
2231	write(,) 'zxtcond=',zxtcond
2232	write(,) 'xtres=',xtres
2233	write(,) 'ep=',ep
2234	stop
2235	endif
2236	#endif
2237	#ifdef ISOVERIF
2238	do ixt=1,ntraciso
2239	call iso_verif_positif(xtres(ixt),'iso_trac_routines 953')
2240	enddo
2241	if (nzone_temp.ge.5) then
2242	if (iso_verif_tag17_q_deltaD_chns(xtres, &
2243	& 'iso_trac_routines 1025').eq.1) then
2244	write(,) 'xt=',xt
2245	write(,) 'zxtcond=',zxtcond
2246	write(,) 'xtres=',xtres
2247	write(,) 'ep=',ep
2248	write(,) 'tcond=',tcond
2249	write(,) 'izone_temp=',izone_temp
2250	stop
2251	endif
2252	endif
2253	! write(,) 'isotrac 1048: sortie de iso_recolorise_condensation'
2254	#endif
2255
2256	return
2257	end subroutine iso_recolorise_condensation
2258
2259	subroutine bassin_map_init_opt20(lat,bassin_map)
2260	USE dimphy, only: klon
2261	use isotrac_mod, only: izone_cont,izone_trop,lim_tag20
2262	#ifdef ISOVERIF
2263	USE isotopes_verif_mod
2264	#endif
2265	implicit none
2266
2267	! inputs
2268	real lat(klon)
2269	! output
2270	integer bassin_map(klon)
2271	! locals
2272	integer i
2273
2274	write(,) 'iso_traceurs_routines 1142: lim_tag20=',lim_tag20
2275	do i=1,klon
2276	if (abs(lat(i)).gt.lim_tag20) then
2277	bassin_map(i)=izone_cont
2278	else
2279	bassin_map(i)=izone_trop
2280	endif
2281	enddo !do i=1,klon
2282
2283	return
2284	end subroutine bassin_map_init_opt20
2285
2286	subroutine isotrac_recolorise_general(xt_seri,t_seri,zx_rh,presnivs)
2287	USE geometry_mod, ONLY : latitude_deg
2288	USE dimphy, only: klon,klev
2289	use isotrac_mod, only: option_traceurs,boite_map
2290	implicit none
2291
2292	! inputs
2293	real, dimension(ntraciso,klon,klev), intent(in) :: xt_seri
2294	real, dimension(klon,klev), intent(in) :: t_seri
2295	real, dimension(klon,klev), intent(in) :: zx_rh
2296	real, dimension(klev), intent(in) :: presnivs
2297
2298	if (option_traceurs.eq.4) then
2299	call isotrac_recolorise_tmin(xt_seri,t_seri)
2300	elseif ((option_traceurs.eq.5).or. &
2301	& (option_traceurs.eq.21)) then
2302	call isotrac_recolorise_boite(xt_seri,boite_map)
2303	elseif (option_traceurs.eq.13) then
2304	call isotrac_recolorise_tmin_sfrev(xt_seri,t_seri)
2305	elseif (option_traceurs.eq.14) then
2306	call isotrac_recolorise_saturation(xt_seri,zx_rh,latitude_deg,presnivs)
2307	elseif (option_traceurs.eq.20) then
2308	call isotrac_recolorise_extra(xt_seri,latitude_deg)
2309	endif !if (option_traceurs.eq.4) then
2310
2311	return
2312	end subroutine isotrac_recolorise_general
2313
2314
2315
2316
2317	subroutine iso_verif_traceur_jbid_vect(x,n,m)
2318	USE isotopes_mod, ONLY: bidouille_anti_divergence,iso_eau,ridicule
2319	!use isotrac_mod, only: ntraceurs_zone=>nzone
2320	USE infotrac_phy, ONLY: ntraceurs_zone=>nzone
2321	implicit none
2322
2323	! version vectrisée de iso_verif_traceur_jbidouille
2324	! inputs
2325	integer n,m
2326	real x(ntraciso,n,m)
2327
2328	! locals
2329	integer iiso,izone,ixt,i,j
2330	real xtractot(n,m)
2331
2332	if (bidouille_anti_divergence) then
2333	do iiso=1,niso
2334
2335	do j=1,m
2336	do i=1,n
2337	xtractot(i,j)=0.0
2338	enddo !do j=1,m
2339	enddo !do j=1,m
2340
2341	do izone=1,ntraceurs_zone
2342	ixt=index_trac(izone,iiso)
2343	do j=1,m
2344	do i=1,n
2345	xtractot(i,j)=xtractot(i,j)+x(ixt,i,j)
2346	enddo !do j=1,m
2347	enddo !do j=1,m
2348	enddo !do izone=1,ntraceurs_zone
2349
2350	! on réajuste pour que les traceurs fasses bien la somme
2351	! des traceurs
2352	do izone=1,ntraceurs_zone
2353	ixt=index_trac(izone,iiso)
2354	do j=1,m
2355	do i=1,n
2356	! if (abs(xtractot(i,j)).gt.ridicule*10) then
2357	if (abs(xtractot(i,j)).gt.ridicule) then
2358	! modif le 19 fev 2011
2359	x(ixt,i,j)=x(ixt,i,j)/xtractot(i,j)*x(iiso,i,j)
2360	endif !if (abs(xtractot(i,j)).gt.ridicule*10) then
2361	enddo !do i=1,n
2362	enddo !do j=1,m
2363	enddo !do izone=1,ntraceurs_zone
2364
2365	! ! ajout le 19 fev 2011
2366	! ! on rend plutot les vérifs plus strictes
2367	! ixt=index_trac(izone_poubelle,iiso)
2368	! do j=1,m
2369	! do i=1,n
2370	! if ((abs(xtractot(i,j)).lt.1e-18).and.
2371	! : (x(iiso,i,j).gt.ridicule)) then
2372	! x(ixt,i,j)=x(iiso,i,j)
2373	! endif !if (abs(xtractot(i,j)).gt.ridicule*10) then
2374	! enddo ! do i=1,n
2375	! enddo !do j=1,m
2376
2377	enddo !do iiso=1,ntraceurs_iso
2378	endif !if (bidouille_anti_divergence) then
2379
2380	end subroutine iso_verif_traceur_jbid_vect
2381
2382	subroutine iso_verif_traceur_jbidouille(x)
2383	USE isotopes_mod, ONLY: bidouille_anti_divergence,iso_eau,ridicule
2384	implicit none
2385
2386	! on réajuste aussi les valeurs des traceurs pour la
2387	! conservation de la masse, dans le cas bidouille
2388
2389	! inputs
2390	real x(ntraciso)
2391
2392	! locals
2393	integer iiso,izone,ixt
2394	real xtractot
2395
2396	if (bidouille_anti_divergence) then
2397
2398	do iiso=1,niso
2399
2400	xtractot=0.0
2401	do izone=1,ntraceurs_zone
2402	ixt=index_trac(izone,iiso)
2403	xtractot=xtractot+x(ixt)
2404	enddo !do izone=1,ntraceurs_zone
2405
2406	! on réajuste pour que les traceurs fasses bien la somme
2407	! des traceurs
2408	if (abs(xtractot).gt.ridicule*10) then
2409	do izone=1,ntraceurs_zone
2410	ixt=index_trac(izone,iiso)
2411	x(ixt)=x(ixt)/xtractot*x(iiso)
2412	enddo !do izone=1,ntraceurs_zone
2413	endif
2414
2415	enddo !do iiso=1,ntraceurs_iso
2416	endif !if (bidouille_anti_divergence) then
2417
2418	end subroutine iso_verif_traceur_jbidouille
2419
2420
2421	subroutine iso_verif_traceur_jbid_pos(x)
2422	USE isotopes_mod, ONLY: bidouille_anti_divergence,iso_eau,ridicule
2423	!#ifdef ISOVERIF
2424	! use isotopes_verif_mod, only: iso_verif_traceur_pbidouille
2425	!#endif
2426	implicit none
2427
2428	! on réajuste les valeurs des traceurs pour qu'il n'y ai pas de
2429	! valeurs négatives. Si valeurs négatives -> on pompe les autres
2430	! traceurs
2431	! attention: fait la même chose pour tous les isos -> peut
2432	! induire des fractionnements.
2433	! Pour ne pas induire des fractionnements, prendre
2434	! iso_verif_traceur_jbid_pos2
2435	! avantage de cette subroutine: conserve la masse en isotopes
2436	! légers, ce qui nest pas le cas de pos2
2437
2438	! inputs
2439	real x(ntraciso)
2440
2441	! locals
2442	integer iiso,izone,ixt
2443	real xtractot,xtractotprec
2444
2445	if (bidouille_anti_divergence) then
2446
2447	! write(,) 'pgm 532 tmp: x=',x
2448	do iiso=1,niso
2449
2450	xtractot=0.0
2451	xtractotprec=0.0
2452	do izone=1,ntraceurs_zone
2453	ixt=index_trac(izone,iiso)
2454	xtractotprec=xtractotprec+x(ixt)
2455	x(ixt)=max(x(ixt),0.0)
2456	xtractot=xtractot+x(ixt)
2457	enddo !do izone=1,ntraceurs_zone
2458	! write(,) 'iiso,xtractotprec,xtractot=',
2459	! : iiso,xtractotprec,xtractot
2460
2461	if (xtractot.gt.xtractotprec) then
2462	! on réajuste pour que les traceurs fasses bien la somme
2463	! des traceurs
2464	if (abs(xtractot).gt.ridicule) then
2465	do izone=1,ntraceurs_zone
2466	ixt=index_trac(izone,iiso)
2467	x(ixt)=x(ixt)*xtractotprec/xtractot
2468	enddo !do izone=1,ntraceurs_zone
2469	! on modifie aussi l'isotope de base si lui aussi était
2470	! négatif
2471	! x(iiso)=xtractot
2472	else !if (abs(xtractot).gt.ridicule) then
2473	! normallement, valeurs restantes très faibles
2474	! on ne fait rien.
2475	! on met juste un max
2476	x(iiso)=max(x(iiso),0.0)
2477	do izone=1,ntraceurs_zone
2478	ixt=index_trac(izone,iiso)
2479	x(ixt)=max(x(ixt),0.0)
2480	enddo !do izone=1,ntraceurs_zone
2481	endif !if (abs(xtractot).gt.ridicule) then
2482	endif !if (xtractot.gt.xtractotprec) then
2483
2484	enddo !do iiso=1,ntraceurs_iso
2485	#ifdef ISOVERIF
2486	call iso_verif_traceur_pbidouille(x,'iso_verif_trac 558')
2487	#else
2488	call iso_verif_traceur_jbidouille(x)
2489	#endif
2490	endif !if (bidouille_anti_divergence) then
2491
2492	end subroutine iso_verif_traceur_jbid_pos
2493
2494	subroutine iso_verif_traceur_jbid_pos_vect(n,m,x)
2495	USE isotopes_mod, ONLY: bidouille_anti_divergence,iso_eau,ridicule
2496	#ifdef ISOVERIF
2497	USE isotopes_verif_mod
2498	#endif
2499	implicit none
2500
2501	! inputs
2502	integer n,m
2503	real x(ntraciso,n,m)
2504
2505	! locals
2506	integer iiso,izone,ixt
2507	real xtractot(n,m),xtractotprec(n,m)
2508	integer i,j
2509
2510	if (bidouille_anti_divergence) then
2511
2512	! write(,) 'pgm 532 tmp: x=',x
2513
2514	do iiso=1,niso
2515	do j=1,m
2516	do i=1,n
2517	xtractot(i,j)=0.0
2518	xtractotprec(i,j)=0.0
2519	enddo !do j=1,m
2520	enddo !do i=1,n
2521	do izone=1,ntraceurs_zone
2522	ixt=index_trac(izone,iiso)
2523
2524	do j=1,m
2525	do i=1,n
2526	xtractotprec(i,j)=xtractotprec(i,j)+x(ixt,i,j)
2527	x(ixt,i,j)=max(x(ixt,i,j),0.0)
2528	xtractot(i,j)=xtractot(i,j)+x(ixt,i,j)
2529	enddo !do i=1,n
2530	enddo !do j=1,m
2531
2532	enddo !do izone=1,ntraceurs_zone
2533	! write(,) 'iiso,xtractotprec,xtractot=',
2534	! : iiso,xtractotprec,xtractot
2535
2536	do j=1,m
2537	do i=1,n
2538	if (xtractot(i,j).gt.xtractotprec(i,j)) then
2539	! on réajuste pour que les traceurs fasses bien la somme
2540	! des traceurs
2541	if (abs(xtractot(i,j)).gt.ridicule) then
2542	do izone=1,ntraceurs_zone
2543	ixt=index_trac(izone,iiso)
2544	x(ixt,i,j)=x(ixt,i,j)*xtractotprec(i,j)/xtractot(i,j)
2545	enddo !do izone=1,ntraceurs_zone
2546	! on modifie aussi l'isotope de base si lui aussi était
2547	! négatif
2548	! x(iiso)=xtractot
2549	else !if (abs(xtractot).gt.ridicule) then
2550	! normallement, valeurs restantes très faibles
2551	! on ne fait rien.
2552	! on met juste un max
2553	x(iiso,i,j)=max(x(iiso,i,j),0.0)
2554	do izone=1,ntraceurs_zone
2555	ixt=index_trac(izone,iiso)
2556	x(ixt,i,j)=max(x(ixt,i,j),0.0)
2557	enddo !do izone=1,ntraceurs_zone
2558	endif !if (abs(xtractot).gt.ridicule) then
2559	endif !if (xtractot.gt.xtractotprec) then
2560	enddo !do i=1,n
2561	enddo !do j=1,m
2562
2563
2564	enddo !do iiso=1,ntraceurs_iso
2565	#ifdef ISOVERIF
2566	call iso_verif_traceur_pbid_vect(x,n,m,'iso_verif_trac 558')
2567	#else
2568	call iso_verif_traceur_jbid_vect(x,n,m)
2569	#endif
2570	endif !if (bidouille_anti_divergence) then
2571
2572	end subroutine iso_verif_traceur_jbid_pos_vect
2573
2574	subroutine iso_verif_traceur_jbid_pos2(x,q)
2575	USE isotopes_mod, ONLY: bidouille_anti_divergence,iso_eau,ridicule
2576	#ifdef ISOVERIF
2577	use isotopes_verif_mod
2578	#endif
2579	implicit none
2580
2581	! même but que iso_verif_traceur_jbid_pos, mais n'induit
2582	! pas de fractionnement.
2583	! on regarde si xteau est positif. S'il ne l'est pas, on pompe
2584	! dans les autres tags pour le mettre à 0. On conserve la compo
2585	! iso.
2586	! Pb: ne conserve pas la masse d'isotopes légers.
2587
2588	! inputs
2589	real x(ntraciso),q
2590
2591	! locals
2592	integer iiso,izone,ixt,ieau
2593	real dqtmp,factmp
2594
2595	if (bidouille_anti_divergence) then
2596	! write(,) 'iso_verif_trac 578 tmp: q,xt=',
2597	! : q,x(1:ntraciso)
2598	if (q.gt.0.0) then
2599	dqtmp=0.0
2600	do izone=1,ntraceurs_zone
2601	ieau=index_trac(izone,iso_eau)
2602	if (x(ieau).lt.0.0) then
2603	! write(,) 'local_x<0 pour izone=',izone
2604	dqtmp=dqtmp-x(ieau)
2605	do iiso=1,niso
2606	ixt=index_trac(izone,iiso)
2607	x(ixt) =0.0
2608	enddo !do iiso=1,niso
2609	endif !if (local_xt(ieau,i,k).lt.0.0) then
2610	enddo !do izone=1,ntraceurs_zone
2611	! write(,) 'dqtmp=',dqtmp
2612	if (dqtmp.gt.0.0) then
2613	! write(,) 'iso_verif_trac 593 warning: q,dqtmp,xt=',
2614	! : q,dqtmp,x(1:ntraciso)
2615	! on redistribue la négativité des traceurs dans les
2616	! traceurs positifs
2617	! factmp=(1.0-dqtmp/(local_q(i,k)+dqtmp))
2618	! correction janv 2010
2619	factmp=(q/(q+dqtmp))
2620	! write(,) 'factmp=',factmp
2621	do izone=1,ntraceurs_zone
2622	ieau=index_trac(izone,iso_eau)
2623	if (x(ieau).gt.0.0) then
2624	do iiso=1,niso
2625	ixt=index_trac(izone,iiso)
2626	x(ixt)=x(ixt)*factmp
2627	enddo !do iiso=1,niso
2628	endif !if (local_xt(ieau,i,k).gt.0.0) then
2629	enddo ! do izone=1,ntraceurs_zone
2630	! write(,) 'apres bidouille: xt=',x(1:ntraciso)
2631	endif !if (dqtmp.gt.0.0) then
2632	#ifdef ISOVERIF
2633	call iso_verif_traceur(x,'iso_verif_traceurs 612')
2634	#endif
2635	endif !if (local_q(i,k).lt.0.0) then
2636
2637	#ifdef ISOVERIF
2638	call iso_verif_traceur_pbidouille(x,'iso_verif_trac 625')
2639	#endif
2640	endif ! if (bidouille_anti_divergence) then
2641
2642	end subroutine iso_verif_traceur_jbid_pos2
2643
2644	subroutine iso_verif_traceur_jbid_vect1D(x,n)
2645	USE isotopes_mod, ONLY: bidouille_anti_divergence,iso_eau,ridicule
2646	implicit none
2647
2648	! version vectrisée de iso_verif_traceur_jbidouille
2649
2650	! inputs
2651	integer n
2652	real x(ntraciso,n)
2653
2654	! locals
2655	integer iiso,izone,ixt,i
2656	real xtractot
2657
2658	if (bidouille_anti_divergence) then
2659
2660	do i=1,n
2661	do iiso=1,niso
2662
2663	xtractot=0.0
2664	do izone=1,ntraceurs_zone
2665	ixt=index_trac(izone,iiso)
2666	xtractot=xtractot+x(ixt,i)
2667	enddo !do izone=1,ntraceurs_zone
2668
2669	! on réajuste pour que les traceurs fasses bien la somme
2670	! des traceurs
2671	if (abs(xtractot).gt.ridicule*10) then
2672	do izone=1,ntraceurs_zone
2673	ixt=index_trac(izone,iiso)
2674	x(ixt,i)=x(ixt,i)/xtractot*x(iiso,i)
2675	enddo !do izone=1,ntraceurs_zone
2676	endif
2677
2678	enddo !do iiso=1,ntraceurs_iso
2679	enddo !do i=1,n
2680	endif !if (bidouille_anti_divergence) then
2681
2682	end subroutine iso_verif_traceur_jbid_vect1D
2683
2684	! on met ces routines ici pour éviter dépendances circulaires
2685	#ifdef ISOVERIF
2686
2687	subroutine iso_verif_traceur_pbidouille(x,err_msg)
2688	use isotopes_verif_mod
2689	implicit none
2690	! vérifier des choses sur les traceurs
2691	! * toutes les zones donne t l'istope total
2692	! * pas de deltaD aberrant
2693	! on réajuste aussi les valeurs des traceurs pour la
2694	! conservation de la masse, dans le cas bidouille
2695
2696	! on prend les valeurs pas défaut pour
2697	! errmax,errmaxrel,ridicule_trac,deltalimtrac
2698
2699	! inputs
2700	real x(ntraciso)
2701	character() err_msg ! message d''erreur à afficher
2702
2703	! local
2704	!integer iso_verif_traceur_pbid_ns
2705
2706	if (iso_verif_traceur_pbid_ns(x,err_msg).eq.1) then
2707	stop
2708	endif
2709
2710	end subroutine iso_verif_traceur_pbidouille
2711
2712	function iso_verif_traceur_pbid_ns(x,err_msg)
2713	use isotopes_mod, ONLY: iso_HDO,bidouille_anti_divergence
2714	use isotrac_mod, only: ridicule_trac
2715	use isotopes_verif_mod
2716	implicit none
2717	! vérifier des choses sur les traceurs
2718	! * toutes les zones donne t l'istope total
2719	! * pas de deltaD aberrant
2720	! on réajuste aussi les valeurs des traceurs pour la
2721	! conservation de la masse, dans le cas bidouille
2722
2723	! on prend les valeurs pas défaut pour
2724	! errmax,errmaxrel,ridicule_trac,deltalimtrac
2725
2726	! inputs
2727	real x(ntraciso)
2728	character() err_msg ! message d''erreur à afficher
2729
2730	! output
2731	integer iso_verif_traceur_pbid_ns
2732
2733	! locals
2734	!integer iso_verif_traceur_noNaN_nostop
2735	!integer iso_verif_tracm_choix_nostop
2736	!integer iso_verif_tracdD_choix_nostop
2737	integer iiso,izone,ixt
2738	real xtractot
2739
2740	! verif noNaN
2741	iso_verif_traceur_pbid_ns=0
2742
2743	if (iso_verif_traceur_noNaN_nostop(x,err_msg).eq.1) then
2744	! stop
2745	iso_verif_traceur_pbid_ns=1
2746	endif
2747
2748	! verif masse
2749	if (iso_verif_tracm_choix_nostop(x,err_msg, &
2750	& errmax10,errmaxrel50).eq.1) then
2751	! on est plus laxiste car ça vient en général après une
2752	! bidouille pour iso_eau normal
2753	! stop
2754	iso_verif_traceur_pbid_ns=1
2755	endif
2756
2757	if (bidouille_anti_divergence) then
2758	! on réajuste pour que les traceurs fasses bien la somme
2759	! des traceurs
2760	call iso_verif_traceur_jbidouille(x)
2761	endif !if (bidouille_anti_divergence) then
2762
2763	! verif deltaD
2764	if (iso_HDO.gt.0) then
2765	if (iso_verif_tracdD_choix_nostop(x,err_msg, &
2766	& ridicule_trac,deltalimtrac).eq.1) then
2767	! stop
2768	iso_verif_traceur_pbid_ns=1
2769	endif
2770	endif !if (iso_HDO.gt.0) then
2771
2772	end function iso_verif_traceur_pbid_ns
2773
2774	subroutine iso_verif_traceur_pbid_vect(x,n,m,err_msg)
2775	use isotopes_mod, ONLY: iso_HDO,bidouille_anti_divergence
2776	use isotopes_verif_mod
2777	implicit none
2778
2779	! inputs
2780	integer n,m
2781	real x(ntraciso,n,m)
2782	character() err_msg ! message d''erreur à afficher
2783
2784	! locals
2785	integer iiso,izone,ixt
2786	real xtractot
2787
2788	! verif noNaN
2789	call iso_verif_traceur_noNaN_vect(x,n,m,err_msg)
2790
2791	! verif masse
2792	call iso_verif_trac_masse_vect(x,n,m,err_msg,errmax*10, &
2793	& errmaxrel*50)
2794
2795	if (bidouille_anti_divergence) then
2796	! on réajuste pour que les traceurs fasses bien la somme
2797	! des traceurs
2798	call iso_verif_traceur_jbid_vect(x,n,m)
2799	endif !if (bidouille_anti_divergence) then
2800
2801	! verif deltaD
2802	if (iso_HDO.gt.0) then
2803	call iso_verif_tracdd_vect(x,n,m,err_msg)
2804	endif
2805
2806	end subroutine iso_verif_traceur_pbid_vect
2807	#endif
2808
2809	END MODULE isotrac_routines_mod
2810	#endif
2811	#endif

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

Download in other formats:

Original Format