source: LMDZ6/trunk/libf/phylmd/calltherm_mod.F90 @ 6099

Last change on this file since 6099 was 6099, checked in by fhourdin, 6 weeks ago

phylmdiso compil. + sep init/compute in physiq_mod

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