source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/calltherm.F90 @ 5097

Last change on this file since 5097 was 4727, checked in by idelkadi, 13 months ago

Merged trunk changes -r4488:4726 LMDZ_ECRad branch

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