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

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