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

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

Initial import of the physics wih isotopes from Camille Risi
CR

File size: 18.4 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              CALL abort_gcm('calltherm 173','isos pas prevus ici',1)
224#endif
225            CALL thermcell_2002(klon,klev,zdt,iflag_thermals   &
226     &      ,pplay,paprs,pphi  &
227     &      ,u_seri,v_seri,t_seri,q_seri  &
228     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
229     &      ,zfm_therm,zentr_therm,fraca,zw2  &
230     &      ,r_aspect_thermals,30.,w2di_thermals  &
231     &      ,tau_thermals)
232          else if (iflag_thermals.eq.2) then
233#ifdef ISO
234              CALL abort_gcm('calltherm 186','isos pas prevus ici',1) 
235#endif
236            CALL thermcell_sec(klon,klev,zdt  &
237     &      ,pplay,paprs,pphi,zlev  &
238     &      ,u_seri,v_seri,t_seri,q_seri  &
239     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
240     &      ,zfm_therm,zentr_therm  &
241     &      ,r_aspect_thermals,30.,w2di_thermals  &
242     &      ,tau_thermals)
243          else if (iflag_thermals.eq.3) then
244#ifdef ISO
245              write(*,*) 'calltherm 199: isos pas prévus ici'
246              stop
247#endif
248            CALL thermcell(klon,klev,zdt  &
249     &      ,pplay,paprs,pphi  &
250     &      ,u_seri,v_seri,t_seri,q_seri  &
251     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
252     &      ,zfm_therm,zentr_therm  &
253     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
254     &      ,tau_thermals)
255          else if (iflag_thermals.eq.10) then
256#ifdef ISO
257              CALL abort_gcm('calltherm 212','isos pas prevus ici',1) 
258#endif
259            CALL thermcell_eau(klon,klev,zdt  &
260     &      ,pplay,paprs,pphi  &
261     &      ,u_seri,v_seri,t_seri,q_seri  &
262     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
263     &      ,zfm_therm,zentr_therm  &
264     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
265     &      ,tau_thermals)
266#ifdef ISO
267              CALL abort_gcm('calltherm 267','isos pas prevus ici',1) 
268#endif
269          else if (iflag_thermals.eq.11) then
270              abort_message = 'cas non prevu dans calltherm'
271              CALL abort_physic (modname,abort_message,1)
272
273!           CALL thermcell_pluie(klon,klev,zdt  &
274!   &      ,pplay,paprs,pphi,zlev  &
275!    &      ,u_seri,v_seri,t_seri,q_seri  &
276!    &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
277!    &      ,zfm_therm,zentr_therm,zqla  &
278!    &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
279!    &      ,tau_thermals,3)
280          else if (iflag_thermals.eq.12) then
281#ifdef ISO
282              CALL abort_gcm('calltherm 282','isos pas prevus ici',1) 
283#endif
284            CALL calcul_sec(klon,klev,zdt  &
285     &      ,pplay,paprs,pphi,zlev  &
286     &      ,u_seri,v_seri,t_seri,q_seri  &
287     &      ,zmax_sec,wmax_sec,zw_sec,lmix_sec  &
288     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
289     &      ,tau_thermals)
290          else if (iflag_thermals==13.or.iflag_thermals==14) then
291#ifdef ISO
292              CALL abort_gcm('calltherm 292','isos pas prevus ici',1) 
293#endif
294            CALL thermcellV0_main(itap,klon,klev,zdt  &
295     &      ,pplay,paprs,pphi,debut  &
296     &      ,u_seri,v_seri,t_seri,q_seri  &
297     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
298     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
299     &      ,ratqscth,ratqsdiff,zqsatth  &
300     &      ,r_aspect_thermals,l_mix_thermals  &
301     &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
302     &      ,zmax0,f0,zw2,fraca)
303          else if (iflag_thermals>=15.and.iflag_thermals<=18) then
304
305!            print*,'THERM iflag_thermas_ed=',iflag_thermals_ed
306            CALL thermcell_main(itap,klon,klev,zdt  &
307     &      ,pplay,paprs,pphi,debut  &
308     &      ,u_seri,v_seri,t_seri,q_seri  &
309     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
310     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
311     &      ,ratqscth,ratqsdiff,zqsatth  &
312!    &      ,r_aspect_thermals,l_mix_thermals &
313!    &      ,tau_thermals,iflag_thermals_ed,iflag_coupl &
314     &      ,Ale,Alp,lalim_conv,wght_th &
315     &      ,zmax0,f0,zw2,fraca,ztv,zpspsk &
316     &      ,ztla,zthl &
317!!! nrlmd le 10/04/2012
318     &      ,pbl_tke,pctsrf,omega,airephy &
319     &      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
320     &      ,n2,s2,ale_bl_stat &
321     &      ,therm_tke_max,env_tke_max &
322     &      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
323     &      ,alp_bl_conv,alp_bl_stat &
324!!! fin nrlmd le 10/04/2012
325     &      ,ztva &
326#ifdef ISO         
327     &      ,xt_seri,d_xt_the &
328#endif         
329     &   )
330           if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
331         else
332           abort_message = 'Cas des thermiques non prevu'
333           CALL abort_physic (modname,abort_message,1)
334         endif
335
336! Attention : les noms sont contre intuitif.
337! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille.
338! Il aurait mieux valu avoir un nobidouille_stratocu
339! Et pour simplifier :
340! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15)
341! Ce serait bien de changer, mai en prenant le temps de vérifier que ca
342! fait bien ce qu'on croit.
343
344       flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18
345
346! Calcul a posteriori du niveau max des thermiques pour les schémas qui
347! ne la sortent pas.
348      if (iflag_thermals<=12.or.iflag_thermals>=1000) then
349         lmax(:)=1
350         do k=1,klev-1
351            zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
352         enddo
353         do k=1,klev-1
354            do i=1,klon
355               if (zfm_therm(i,k+1)>0.) lmax(i)=k
356            enddo
357         enddo
358      endif
359
360      fact(:)=0.
361      DO i=1,klon
362       logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
363       IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
364      ENDDO
365
366     DO k=1,klev
367!  transformation de la derivee en tendance
368            d_t_the(:,k)=d_t_the(:,k)*dtime*fact(:)
369            d_u_the(:,k)=d_u_the(:,k)*dtime*fact(:)
370            d_v_the(:,k)=d_v_the(:,k)*dtime*fact(:)
371            d_q_the(:,k)=d_q_the(:,k)*dtime*fact(:)
372            fm_therm(:,k)=fm_therm(:,k)  &
373     &      +zfm_therm(:,k)*fact(:)
374            entr_therm(:,k)=entr_therm(:,k)  &
375     &       +zentr_therm(:,k)*fact(:)
376            detr_therm(:,k)=detr_therm(:,k)  &
377     &       +zdetr_therm(:,k)*fact(:)
378#ifdef ISO
379            do ixt=1,ntraciso
380              d_xt_the(ixt,:,k)=d_xt_the(ixt,:,k)*dtime*fact(:)
381            enddo
382#endif
383      ENDDO
384       fm_therm(:,klev+1)=0.
385
386
387
388!  accumulation de la tendance
389            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
390            d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
391            d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
392            d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
393#ifdef ISO
394            d_xt_ajs(:,:,:)=d_xt_ajs(:,:,:)+d_xt_the(:,:,:)
395#endif
396
397!  incrementation des variables meteo
398            t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
399            u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
400            v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
401            qmemoire(:,:)=q_seri(:,:)
402            q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
403#ifdef ISO
404            xtmemoire(:,:,:)=xt_seri(:,:,:)
405            xt_seri(:,:,:) = xt_seri(:,:,:) + d_xt_the(:,:,:)
406#ifdef ISOVERIF
407      write(*,*) 'calltherm 350 tmp: ajout d_xt_the'
408      if (iso_HDO.gt.0) then
409!      i=479
410!      k=4
411!      write(*,*) 'xt_seri(iso_hdo,i,k),q_seri(i,k)=', &
412!     &   xt_seri(iso_hdo,i,k),q_seri(i,k)
413!      write(*,*) 'd_xt_the(iso_hdo,i,k),d_q_the(i,k)=', &
414!     &   d_xt_the(iso_hdo,i,k),d_q_the(i,k)
415      call iso_verif_aberrant_enc_vect2D( &
416     &        xt_seri,q_seri, &
417     &        'calltherm 353, apres ajout d_xt_the',ntraciso,klon,klev)
418      endif     
419#endif
420#endif
421           if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
422
423       DO i=1,klon
424            fm_therm(i,klev+1)=0.
425            Ale_bl(i)=Ale_bl(i)+Ale(i)/REAL(nsplit_thermals)
426!            write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
427            Alp_bl(i)=Alp_bl(i)+Alp(i)/REAL(nsplit_thermals)
428!            write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
429        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)
430       ENDDO
431
432!IM 060508 marche pas comme cela !!!        enddo ! isplit
433
434!   tests sur les valeurs negatives de l'eau
435         nbptspb=0
436            DO k = 1, klev
437            DO i = 1, klon
438               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
439               if (logexpr2(i,k)) then
440                q_seri(i,k)=1.e-15
441                nbptspb=nbptspb+1
442#ifdef ISO
443                do ixt=1,ntraciso
444                  xt_seri(ixt,i,k)=1.e-15*(xtmemoire(ixt,i,k)/qmemoire(i,k))
445                enddo
446#endif
447!                if (prt_level.ge.10) then
448!                  print*,'WARN eau<0 apres therm i=',i,'  k=',k  &
449!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k),  &
450!    &         'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
451                 endif
452            ENDDO
453            ENDDO
454#ifdef ISO
455#ifdef ISOVERIF
456      if (iso_HDO.gt.0) then
457      call iso_verif_aberrant_enc_vect2D( &
458     &        xt_seri,q_seri, &
459     &        'calltherm 393, apres bidouille q<0',ntraciso,klon,klev)
460      endif     
461#endif
462#endif
463
464        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
465! tests sur les valeurs de la temperature
466        nbptspb=0
467            DO k = 1, klev
468            DO i = 1, klon
469               logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
470               if (logexpr2(i,k)) nbptspb=nbptspb+1
471!              if ((t_seri(i,k).lt.50.) .or.  &
472!    &              (t_seri(i,k).gt.370.)) then
473!                 print*,'WARN temp apres therm i=',i,'  k=',k  &
474!    &         ,' t_seri',t_seri(i,k)
475!              CALL abort
476!              endif
477            ENDDO
478            ENDDO
479        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
480         enddo ! isplit
481
482!
483!***************************************************************
484!     calcul du flux ascencant conservatif
485!            print*,'<<<<calcul flux ascendant conservatif'
486
487      fmc_therm=0.
488               do k=1,klev
489            do i=1,klon
490                  if (entr_therm(i,k).gt.0.) then
491                     fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
492                  else
493                     fmc_therm(i,k+1)=fmc_therm(i,k)
494                  endif
495                  detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1))  &
496     &                 -(fmc_therm(i,k)-fm_therm(i,k))
497               enddo
498            enddo
499     
500     
501!****************************************************************
502!     calcul de l'humidite dans l'ascendance
503!      print*,'<<<<calcul de lhumidite dans thermique'
504!CR:on ne le calcule que pour le cas sec
505      if (iflag_thermals.le.11) then     
506      do i=1,klon
507         zqasc(i,1)=q_seri(i,1)
508         do k=2,klev
509            if (fmc_therm(i,k+1).gt.1.e-6) then
510               zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1)  &
511     &              +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
512!CR:test on asseche le thermique
513!               zqasc(i,k)=zqasc(i,k)/2.
514!            else
515!               zqasc(i,k)=q_seri(i,k)
516            endif
517         enddo
518       enddo
519     
520
521!     calcul de l'eau condensee dans l'ascendance
522!             print*,'<<<<calcul de leau condensee dans thermique'
523             do i=1,klon
524                do k=1,klev
525                   clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
526                   if (clwcon0(i,k).lt.0. .or.   &
527     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
528                      clwcon0(i,k)=0.
529                   endif
530                enddo
531             enddo
532       else
533              do i=1,klon
534                do k=1,klev
535                   clwcon0(i,k)=zqla(i,k) 
536                   if (clwcon0(i,k).lt.0. .or.   &
537     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
538                   clwcon0(i,k)=0.
539                   endif
540                enddo
541             enddo
542       endif
543!*******************************************************************   
544
545
546!jyg  Protection contre les temperatures nulles
547          do i=1,klon
548             do k=1,klev
549                if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
550             enddo
551          enddo
552
553
554      return
555
556      end
Note: See TracBrowser for help on using the repository browser.