source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/calltherm.F90 @ 3663

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

Add modification for isotopes

  • Property svn:executable set to *
File size: 18.6 KB
Line 
1!
2! $Id: calltherm.F90 2346 2015-08-21 15:13:46Z emillour $
3!
4      subroutine calltherm(dtime  &
5     &      ,pplay,paprs,pphi,weak_inversion  &
6     &      ,u_seri,v_seri,t_seri,q_seri,zqsat,debut  &
7     &      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs  &
8     &      ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,lmax,ratqscth,  &
9     &       ratqsdiff,zqsatth,Ale_bl,Alp_bl,lalim_conv,wght_th, &
10     &       zmax0,f0,zw2,fraca,ztv,zpspsk,ztla,zthl &
11!!! nrlmd le 10/04/2012
12     &      ,pbl_tke,pctsrf,omega,airephy &
13     &      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
14     &      ,n2,s2,ale_bl_stat &
15     &      ,therm_tke_max,env_tke_max &
16     &      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
17     &      ,alp_bl_conv,alp_bl_stat &
18!!! fin nrlmd le 10/04/2012
19     &      ,zqla,ztva &
20#ifdef ISO         
21     &      ,xt_seri,d_xt_ajs &
22#ifdef DIAGISO         
23     &      ,q_the,xt_the &
24#endif
25#endif         
26     &   )
27
28      USE dimphy
29      USE indice_sol_mod
30      USE print_control_mod, ONLY: prt_level,lunout
31#ifdef ISO
32      use infotrac_phy, ONLY: ntraciso
33#ifdef ISOVERIF
34      USE isotopes_mod, ONLY: iso_eau,iso_HDO
35      USE isotopes_verif_mod, ONLY: iso_verif_aberrant_enc_vect2D, &
36        iso_verif_egalite_vect2D
37#endif   
38#endif
39
40      implicit none
41      include "thermcell.h"
42
43
44!IM 140508
45      INTEGER, SAVE ::  itap
46!$OMP THREADPRIVATE(itap)
47      REAL dtime
48      LOGICAL debut
49      LOGICAL logexpr0, logexpr2(klon,klev), logexpr1(klon)
50      REAL fact(klon)
51      INTEGER nbptspb
52
53      REAL u_seri(klon,klev),v_seri(klon,klev)
54      REAL t_seri(klon,klev),q_seri(klon,klev),qmemoire(klon,klev)
55      REAL weak_inversion(klon)
56      REAL paprs(klon,klev+1)
57      REAL pplay(klon,klev)
58      REAL pphi(klon,klev)
59      real zlev(klon,klev+1)
60!test: on sort lentr et a* pour alimenter KE
61      REAL wght_th(klon,klev)
62      INTEGER lalim_conv(klon)
63      REAL zw2(klon,klev+1),fraca(klon,klev+1)
64
65!FH Update Thermiques
66      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
67      REAL d_u_ajs(klon,klev),d_v_ajs(klon,klev)
68      real fm_therm(klon,klev+1)
69      real entr_therm(klon,klev),detr_therm(klon,klev)
70
71!********************************************************
72!     declarations
73      LOGICAL flag_bidouille_stratocu
74      real fmc_therm(klon,klev+1),zqasc(klon,klev)
75      real zqla(klon,klev)
76      real zqta(klon,klev)
77      real ztv(klon,klev),ztva(klon,klev)
78      real zpspsk(klon,klev)
79      real ztla(klon,klev)
80      real zthl(klon,klev)
81      real wmax_sec(klon)
82      real zmax_sec(klon)
83      real f_sec(klon)
84      real detrc_therm(klon,klev)
85! FH WARNING : il semble que ces save ne servent a rien
86!     save fmc_therm, detrc_therm
87      real clwcon0(klon,klev)
88      real zqsat(klon,klev)
89      real zw_sec(klon,klev+1)
90      integer lmix_sec(klon)
91      integer lmax(klon)
92      real ratqscth(klon,klev)
93      real ratqsdiff(klon,klev)
94      real zqsatth(klon,klev) 
95!nouvelles variables pour la convection
96      real Ale_bl(klon)
97      real Alp_bl(klon)
98      real Ale(klon)
99      real Alp(klon)
100!RC
101      !on garde le zmax du pas de temps precedent
102      real zmax0(klon), f0(klon)
103
104!!! nrlmd le 10/04/2012
105      real pbl_tke(klon,klev+1,nbsrf)
106      real pctsrf(klon,nbsrf)
107      real omega(klon,klev)
108      real airephy(klon)
109      real zlcl_th(klon),fraca0(klon),w0(klon),w_conv(klon)
110      real therm_tke_max0(klon),env_tke_max0(klon)
111      real n2(klon),s2(klon)
112      real ale_bl_stat(klon)
113      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
114      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
115!!! fin nrlmd le 10/04/2012
116
117!********************************************************
118
119
120! variables locales
121      REAL d_t_the(klon,klev), d_q_the(klon,klev)
122      REAL d_u_the(klon,klev),d_v_the(klon,klev)
123!
124      real zfm_therm(klon,klev+1),zdt
125      real zentr_therm(klon,klev),zdetr_therm(klon,klev)
126! FH A VERIFIER : SAVE INUTILES
127!      save zentr_therm,zfm_therm
128
129      character (len=20) :: modname='calltherm'
130      character (len=80) :: abort_message
131
132      integer i,k
133      logical, save :: first=.true.
134
135#ifdef ISO
136      REAL xt_seri(ntraciso,klon,klev),xtmemoire(ntraciso,klon,klev)
137      REAL d_xt_ajs(ntraciso,klon,klev)
138      real d_xt_the(ntraciso,klon,klev)
139#ifdef DIAGISO
140      real q_the(klon,klev)
141      real xt_the(ntraciso,klon,klev)
142#endif
143      real qprec(klon,klev)
144      integer ixt
145#endif
146
147
148!$OMP THREADPRIVATE(first)
149!********************************************************
150      if (first) then
151        itap=0
152        first=.false.
153      endif
154
155! Incrementer le compteur de la physique
156     itap   = itap + 1
157
158!  Modele du thermique
159!  ===================
160!         print*,'thermiques: WARNING on passe t au lieu de t_seri'
161
162
163! On prend comme valeur initiale des thermiques la valeur du pas
164! de temps precedent
165         zfm_therm(:,:)=fm_therm(:,:)
166         zdetr_therm(:,:)=detr_therm(:,:)
167         zentr_therm(:,:)=entr_therm(:,:)
168
169! On reinitialise les flux de masse a zero pour le cumul en
170! cas de splitting
171         fm_therm(:,:)=0.
172         entr_therm(:,:)=0.
173         detr_therm(:,:)=0.
174
175         Ale_bl(:)=0.
176         Alp_bl(:)=0.
177         if (prt_level.ge.10) then
178          print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
179         endif
180
181!   tests sur les valeurs negatives de l'eau
182         logexpr0=prt_level.ge.10
183         nbptspb=0
184         do k=1,klev
185            do i=1,klon
186! Attention teste abderr 19-03-09
187!               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
188                logexpr2(i,k)=.not.q_seri(i,k).ge.1.e-15
189               if (logexpr2(i,k)) then
190#ifdef ISO
191                qprec(i,k)=q_seri(i,k)
192#endif               
193                q_seri(i,k)=1.e-15
194                nbptspb=nbptspb+1
195#ifdef ISO
196                do ixt=1,ntraciso
197                  xt_seri(ixt,i,k)=1.e-15*(xt_seri(ixt,i,k)/qprec(i,k))
198                  ! xt_seri(ixt,i,k)=1.e-15*(Rdefault(index_iso(ixt)))
199                enddo
200#endif             
201               endif
202!               if (logexpr0) &
203!    &             print*,'WARN eau<0 avant therm i=',i,'  k=',k  &
204!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k)
205            enddo
206         enddo
207         if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
208
209#ifdef ISO
210#ifdef ISOVERIF
211      if (iso_eau.gt.0) then
212       call iso_verif_egalite_vect2D( &
213     &           xt_seri,q_seri, &
214     &           'calltherm 174',ntraciso,klon,klev)
215      endif !if (iso_eau.gt.0) then
216#endif   
217#endif
218         zdt=dtime/REAL(nsplit_thermals)
219         do isplit=1,nsplit_thermals
220
221          if (iflag_thermals>=1000) then
222#ifdef ISO
223              write(*,*) 'calltherm 173: isos pas prévus ici'
224              stop
225#endif       
226            CALL thermcell_2002(klon,klev,zdt,iflag_thermals   &
227     &      ,pplay,paprs,pphi  &
228     &      ,u_seri,v_seri,t_seri,q_seri  &
229     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
230     &      ,zfm_therm,zentr_therm,fraca,zw2  &
231     &      ,r_aspect_thermals,30.,w2di_thermals  &
232     &      ,tau_thermals)
233          else if (iflag_thermals.eq.2) then
234#ifdef ISO
235              write(*,*) 'calltherm 186: isos pas prévus ici'
236              stop
237#endif           
238            CALL thermcell_sec(klon,klev,zdt  &
239     &      ,pplay,paprs,pphi,zlev  &
240     &      ,u_seri,v_seri,t_seri,q_seri  &
241     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
242     &      ,zfm_therm,zentr_therm  &
243     &      ,r_aspect_thermals,30.,w2di_thermals  &
244     &      ,tau_thermals)
245          else if (iflag_thermals.eq.3) then
246#ifdef ISO
247              write(*,*) 'calltherm 199: isos pas prévus ici'
248              stop
249#endif       
250            CALL thermcell(klon,klev,zdt  &
251     &      ,pplay,paprs,pphi  &
252     &      ,u_seri,v_seri,t_seri,q_seri  &
253     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
254     &      ,zfm_therm,zentr_therm  &
255     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
256     &      ,tau_thermals)
257          else if (iflag_thermals.eq.10) then
258#ifdef ISO
259              write(*,*) 'calltherm 212: isos pas prévus ici'
260              stop
261#endif         
262            CALL thermcell_eau(klon,klev,zdt  &
263     &      ,pplay,paprs,pphi  &
264     &      ,u_seri,v_seri,t_seri,q_seri  &
265     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
266     &      ,zfm_therm,zentr_therm  &
267     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
268     &      ,tau_thermals)
269#ifdef ISO
270            write(*,*) 'calltherm 211: pas prévu pour isotopes'
271            stop
272#endif
273
274          else if (iflag_thermals.eq.11) then
275              abort_message = 'cas non prevu dans calltherm'
276              CALL abort_physic (modname,abort_message,1)
277
278!           CALL thermcell_pluie(klon,klev,zdt  &
279!   &      ,pplay,paprs,pphi,zlev  &
280!    &      ,u_seri,v_seri,t_seri,q_seri  &
281!    &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
282!    &      ,zfm_therm,zentr_therm,zqla  &
283!    &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
284!    &      ,tau_thermals,3)
285          else if (iflag_thermals.eq.12) then
286#ifdef ISO
287              write(*,*) 'calltherm 240: isos pas prévus ici'
288              stop
289#endif             
290            CALL calcul_sec(klon,klev,zdt  &
291     &      ,pplay,paprs,pphi,zlev  &
292     &      ,u_seri,v_seri,t_seri,q_seri  &
293     &      ,zmax_sec,wmax_sec,zw_sec,lmix_sec  &
294     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
295     &      ,tau_thermals)
296          else if (iflag_thermals==13.or.iflag_thermals==14) then
297#ifdef ISO
298              write(*,*) 'calltherm 252: isos pas prévus ici'
299              stop
300#endif             
301            CALL thermcellV0_main(itap,klon,klev,zdt  &
302     &      ,pplay,paprs,pphi,debut  &
303     &      ,u_seri,v_seri,t_seri,q_seri  &
304     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
305     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
306     &      ,ratqscth,ratqsdiff,zqsatth  &
307     &      ,r_aspect_thermals,l_mix_thermals  &
308     &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
309     &      ,zmax0,f0,zw2,fraca)
310          else if (iflag_thermals>=15.and.iflag_thermals<=18) then
311
312!            print*,'THERM iflag_thermas_ed=',iflag_thermals_ed
313            CALL thermcell_main(itap,klon,klev,zdt  &
314     &      ,pplay,paprs,pphi,debut  &
315     &      ,u_seri,v_seri,t_seri,q_seri  &
316     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
317     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
318     &      ,ratqscth,ratqsdiff,zqsatth  &
319!    &      ,r_aspect_thermals,l_mix_thermals &
320!    &      ,tau_thermals,iflag_thermals_ed,iflag_coupl &
321     &      ,Ale,Alp,lalim_conv,wght_th &
322     &      ,zmax0,f0,zw2,fraca,ztv,zpspsk &
323     &      ,ztla,zthl &
324!!! nrlmd le 10/04/2012
325     &      ,pbl_tke,pctsrf,omega,airephy &
326     &      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
327     &      ,n2,s2,ale_bl_stat &
328     &      ,therm_tke_max,env_tke_max &
329     &      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
330     &      ,alp_bl_conv,alp_bl_stat &
331!!! fin nrlmd le 10/04/2012
332     &      ,ztva &
333#ifdef ISO         
334     &      ,xt_seri,d_xt_the &
335#ifdef DIAGISO         
336     &      ,q_the,xt_the &
337#endif         
338#endif         
339     &   )
340           if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
341         else
342           abort_message = 'Cas des thermiques non prevu'
343           CALL abort_physic (modname,abort_message,1)
344         endif
345
346! Attention : les noms sont contre intuitif.
347! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille.
348! Il aurait mieux valu avoir un nobidouille_stratocu
349! Et pour simplifier :
350! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15)
351! Ce serait bien de changer, mai en prenant le temps de vérifier que ca
352! fait bien ce qu'on croit.
353
354       flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18
355
356! Calcul a posteriori du niveau max des thermiques pour les schémas qui
357! ne la sortent pas.
358      if (iflag_thermals<=12.or.iflag_thermals>=1000) then
359         lmax(:)=1
360         do k=1,klev-1
361            zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
362         enddo
363         do k=1,klev-1
364            do i=1,klon
365               if (zfm_therm(i,k+1)>0.) lmax(i)=k
366            enddo
367         enddo
368      endif
369
370      fact(:)=0.
371      DO i=1,klon
372       logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
373       IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
374      ENDDO
375
376     DO k=1,klev
377!  transformation de la derivee en tendance
378            d_t_the(:,k)=d_t_the(:,k)*dtime*fact(:)
379            d_u_the(:,k)=d_u_the(:,k)*dtime*fact(:)
380            d_v_the(:,k)=d_v_the(:,k)*dtime*fact(:)
381            d_q_the(:,k)=d_q_the(:,k)*dtime*fact(:)
382            fm_therm(:,k)=fm_therm(:,k)  &
383     &      +zfm_therm(:,k)*fact(:)
384            entr_therm(:,k)=entr_therm(:,k)  &
385     &       +zentr_therm(:,k)*fact(:)
386            detr_therm(:,k)=detr_therm(:,k)  &
387     &       +zdetr_therm(:,k)*fact(:)
388#ifdef ISO
389            do ixt=1,ntraciso
390              d_xt_the(ixt,:,k)=d_xt_the(ixt,:,k)*dtime*fact(:)
391            enddo
392#endif     
393      ENDDO
394       fm_therm(:,klev+1)=0.
395
396
397
398!  accumulation de la tendance
399            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
400            d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
401            d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
402            d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
403#ifdef ISO
404            d_xt_ajs(:,:,:)=d_xt_ajs(:,:,:)+d_xt_the(:,:,:)
405#endif       
406
407!  incrementation des variables meteo
408            t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
409            u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
410            v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
411            qmemoire(:,:)=q_seri(:,:)
412            q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
413#ifdef ISO
414            xtmemoire(:,:,:)=xt_seri(:,:,:)
415            xt_seri(:,:,:) = xt_seri(:,:,:) + d_xt_the(:,:,:)
416#ifdef ISOVERIF
417      write(*,*) 'calltherm 350 tmp: ajout d_xt_the'
418      if (iso_HDO.gt.0) then
419!      i=479
420!      k=4
421!      write(*,*) 'xt_seri(iso_hdo,i,k),q_seri(i,k)=', &
422!     &   xt_seri(iso_hdo,i,k),q_seri(i,k)
423!      write(*,*) 'd_xt_the(iso_hdo,i,k),d_q_the(i,k)=', &
424!     &   d_xt_the(iso_hdo,i,k),d_q_the(i,k)
425      call iso_verif_aberrant_enc_vect2D( &
426     &        xt_seri,q_seri, &
427     &        'calltherm 353, apres ajout d_xt_the',ntraciso,klon,klev)
428      endif     
429#endif
430#endif         
431
432           if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
433
434       DO i=1,klon
435            fm_therm(i,klev+1)=0.
436            Ale_bl(i)=Ale_bl(i)+Ale(i)/REAL(nsplit_thermals)
437!            write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
438            Alp_bl(i)=Alp_bl(i)+Alp(i)/REAL(nsplit_thermals)
439!            write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
440        if(prt_level.GE.10) print*,'calltherm i Alp_bl Alp Ale_bl Ale',i,Alp_bl(i),Alp(i),Ale_bl(i),Ale(i)
441       ENDDO
442
443!IM 060508 marche pas comme cela !!!        enddo ! isplit
444
445!   tests sur les valeurs negatives de l'eau
446         nbptspb=0
447            DO k = 1, klev
448            DO i = 1, klon
449               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
450               if (logexpr2(i,k)) then
451                q_seri(i,k)=1.e-15
452                nbptspb=nbptspb+1
453#ifdef ISO
454                do ixt=1,ntraciso
455                  xt_seri(ixt,i,k)=1.e-15*(xtmemoire(ixt,i,k)/qmemoire(i,k))
456                enddo
457#endif         
458!                if (prt_level.ge.10) then
459!                  print*,'WARN eau<0 apres therm i=',i,'  k=',k  &
460!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k),  &
461!    &         'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
462                 endif
463            ENDDO
464            ENDDO
465#ifdef ISO
466#ifdef ISOVERIF
467      if (iso_HDO.gt.0) then
468      call iso_verif_aberrant_enc_vect2D( &
469     &        xt_seri,q_seri, &
470     &        'calltherm 393, apres bidouille q<0',ntraciso,klon,klev)
471      endif     
472#endif
473#endif
474
475        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
476! tests sur les valeurs de la temperature
477        nbptspb=0
478            DO k = 1, klev
479            DO i = 1, klon
480               logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
481               if (logexpr2(i,k)) nbptspb=nbptspb+1
482!              if ((t_seri(i,k).lt.50.) .or.  &
483!    &              (t_seri(i,k).gt.370.)) then
484!                 print*,'WARN temp apres therm i=',i,'  k=',k  &
485!    &         ,' t_seri',t_seri(i,k)
486!              CALL abort
487!              endif
488            ENDDO
489            ENDDO
490        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
491         enddo ! isplit
492
493!
494!***************************************************************
495!     calcul du flux ascencant conservatif
496!            print*,'<<<<calcul flux ascendant conservatif'
497
498      fmc_therm=0.
499               do k=1,klev
500            do i=1,klon
501                  if (entr_therm(i,k).gt.0.) then
502                     fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
503                  else
504                     fmc_therm(i,k+1)=fmc_therm(i,k)
505                  endif
506                  detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1))  &
507     &                 -(fmc_therm(i,k)-fm_therm(i,k))
508               enddo
509            enddo
510     
511     
512!****************************************************************
513!     calcul de l'humidite dans l'ascendance
514!      print*,'<<<<calcul de lhumidite dans thermique'
515!CR:on ne le calcule que pour le cas sec
516      if (iflag_thermals.le.11) then     
517      do i=1,klon
518         zqasc(i,1)=q_seri(i,1)
519         do k=2,klev
520            if (fmc_therm(i,k+1).gt.1.e-6) then
521               zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1)  &
522     &              +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
523!CR:test on asseche le thermique
524!               zqasc(i,k)=zqasc(i,k)/2.
525!            else
526!               zqasc(i,k)=q_seri(i,k)
527            endif
528         enddo
529       enddo
530     
531
532!     calcul de l'eau condensee dans l'ascendance
533!             print*,'<<<<calcul de leau condensee dans thermique'
534             do i=1,klon
535                do k=1,klev
536                   clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
537                   if (clwcon0(i,k).lt.0. .or.   &
538     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
539                      clwcon0(i,k)=0.
540                   endif
541                enddo
542             enddo
543       else
544              do i=1,klon
545                do k=1,klev
546                   clwcon0(i,k)=zqla(i,k) 
547                   if (clwcon0(i,k).lt.0. .or.   &
548     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
549                   clwcon0(i,k)=0.
550                   endif
551                enddo
552             enddo
553       endif
554!*******************************************************************   
555
556
557!jyg  Protection contre les temperatures nulles
558          do i=1,klon
559             do k=1,klev
560                if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
561             enddo
562          enddo
563
564
565      return
566
567      end
Note: See TracBrowser for help on using the repository browser.