source: LMDZ6/trunk/libf/phylmdiso/isotrac_routines_mod.F90 @ 4657

Last change on this file since 4657 was 4491, checked in by crisi, 21 months ago

Bug corrections in LMDZiso, especially for water tagging

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