source: LMDZ6/branches/Ocean_skin/libf/phylmdiso/isotrac_routines_mod.F90 @ 5454

Last change on this file since 5454 was 4368, checked in by lguez, 2 years ago

Sync latest trunk changes to Ocean_skin

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