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

Last change on this file since 3927 was 3927, checked in by Laurent Fairhead, 3 years ago

Initial import of the physics wih isotopes from Camille Risi
CR

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