source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/isotrac_routines_mod.F90 @ 3642

Last change on this file since 3642 was 3331, checked in by acozic, 6 years ago

Add modification for isotopes

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