Index: /LMDZ6/trunk/libf/phylmd/calltherm_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/calltherm_mod.F90	(revision 5513)
+++ /LMDZ6/trunk/libf/phylmd/calltherm_mod.F90	(revision 5513)
@@ -0,0 +1,544 @@
+!
+! $Id: calltherm.F90 5390 2024-12-05 16:09:25Z ymeurdesoif $
+!
+MODULE calltherm_mod
+
+CONTAINS
+
+      subroutine calltherm(itap, dtime  &
+     &      ,pplay,paprs,pphi,weak_inversion  &
+     &      ,u_seri_,v_seri_,t_seri_,q_seri_,t_env,q_env,zqsat,debut  &
+     &      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs  &
+     &    ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,lmax,ratqscth,&
+     &       ratqsdiff,zqsatth,ale_bl,alp_bl,lalim_conv,wght_th, &
+     &       zmax0,f0,zw2,fraca,ztv,zpspsk,ztla,zthl &
+!!! nrlmd le 10/04/2012
+     &      ,pbl_tke,pctsrf,omega,airephy &
+     &      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
+     &      ,n2,s2,strig,zcong,ale_bl_stat &
+     &      ,therm_tke_max,env_tke_max &
+     &      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
+     &      ,alp_bl_conv,alp_bl_stat &
+!!! fin nrlmd le 10/04/2012
+     &      ,zqla,ztva &
+#ifdef ISO         
+     &      ,xt_seri,d_xt_ajs &
+#ifdef DIAGISO         
+     &      ,q_the,xt_the &
+#endif
+#endif         
+     &   )
+
+      USE dimphy
+      USE indice_sol_mod
+      USE print_control_mod, ONLY: prt_level,lunout
+      USE lmdz_thermcell_alp, ONLY: thermcell_alp
+      USE lmdz_thermcell_main, ONLY: thermcell_main
+      USE lmdz_thermcell_old, ONLY: thermcell_old=>thermcell, thermcell_2002, thermcell_eau, calcul_sec, thermcell_sec
+#ifdef ISO
+      use infotrac_phy, ONLY: ntiso
+#ifdef ISOVERIF
+      USE isotopes_mod, ONLY: iso_eau,iso_HDO
+      USE isotopes_verif_mod, ONLY: iso_verif_aberrant_enc_vect2D, &
+        iso_verif_egalite_vect2D
+#endif   
+#endif
+      USE clesphys_mod_h
+      USE thermcell_old_mod_h, ONLY: r_aspect_thermals, l_mix_thermals, w2di_thermals
+      implicit none
+
+      INTEGER, INTENT(IN) :: itap
+!IM 140508
+      REAL dtime
+      LOGICAL debut
+      LOGICAL logexpr0, logexpr2(klon,klev), logexpr1(klon)
+      REAL fact(klon)
+      INTEGER nbptspb
+
+      REAL, DIMENSION(klon,klev), INTENT(IN) :: u_seri_,v_seri_
+      REAL, DIMENSION(klon,klev), INTENT(IN) :: t_seri_,q_seri_
+      REAL, DIMENSION(klon,klev), INTENT(IN) :: t_env,q_env
+      REAL, DIMENSION(klon,klev) :: u_seri,v_seri
+      REAL, DIMENSION(klon,klev) :: t_seri,q_seri
+      REAL, DIMENSION(klon,klev) :: qmemoire
+      REAL weak_inversion(klon)
+      REAL paprs(klon,klev+1)
+      REAL pplay(klon,klev)
+      REAL pphi(klon,klev)
+      real zlev(klon,klev+1) 
+!test: on sort lentr et a* pour alimenter KE
+      REAL wght_th(klon,klev)
+      INTEGER lalim_conv(klon)
+      REAL zw2(klon,klev+1),fraca(klon,klev+1)
+
+!FH Update Thermiques
+      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
+      REAL d_u_ajs(klon,klev),d_v_ajs(klon,klev)
+      real fm_therm(klon,klev+1)
+      real entr_therm(klon,klev),detr_therm(klon,klev)
+
+!********************************************************
+!     declarations
+      LOGICAL flag_bidouille_stratocu
+      real fmc_therm(klon,klev+1),zqasc(klon,klev)
+      real zqla(klon,klev)
+      real ztv(klon,klev),ztva(klon,klev)
+      real zpspsk(klon,klev)
+      real ztla(klon,klev)
+      real zthl(klon,klev)
+      real wmax_sec(klon)
+      real zcong(klon)
+      real zmax_sec(klon)
+      real f_sec(klon)
+      real detrc_therm(klon,klev)
+! FH WARNING : il semble que ces save ne servent a rien
+!     save fmc_therm, detrc_therm
+      real clwcon0(klon,klev)
+      real zqsat(klon,klev)
+      real zw_sec(klon,klev+1)
+      integer lmix_sec(klon)
+      integer lmax(klon)
+      real ratqscth(klon,klev)
+      real ratqsdiff(klon,klev)
+      real zqsatth(klon,klev)  
+!nouvelles variables pour la convection
+      real ale_bl(klon)
+      real alp_bl(klon)
+      real ale(klon)
+      real alp(klon)
+!RC
+      !on garde le zmax du pas de temps precedent
+      real zmax0(klon), f0(klon)
+
+!!! nrlmd le 10/04/2012
+      real pbl_tke(klon,klev+1,nbsrf)
+      real pctsrf(klon,nbsrf)
+      real omega(klon,klev)
+      real airephy(klon)
+      real zlcl_th(klon),fraca0(klon),w0(klon),w_conv(klon)
+      real therm_tke_max0(klon),env_tke_max0(klon)
+      real n2(klon),s2(klon),strig(klon)
+      real ale_bl_stat(klon)
+      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
+      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
+!!! fin nrlmd le 10/04/2012
+
+!********************************************************
+
+      real, dimension(klon) :: pcon
+      real, dimension(klon,klev) :: rhobarz,wth3
+      integer,dimension(klon) :: lalim
+      real, dimension(klon,klev+1) :: fm
+      real, dimension(klon,klev) :: alim_star
+      real, dimension(klon) :: zmax
+
+
+
+
+! variables locales
+      REAL d_t_the(klon,klev), d_q_the(klon,klev)
+      REAL d_u_the(klon,klev),d_v_the(klon,klev)
+!
+      real zfm_therm(klon,klev+1),zdt
+      real zentr_therm(klon,klev),zdetr_therm(klon,klev)
+! FH A VERIFIER : SAVE INUTILES
+!      save zentr_therm,zfm_therm
+
+      character (len=20),PARAMETER :: modname='calltherm'
+      character (len=80) :: abort_message
+
+      integer i,k,isplit
+      logical :: new_thermcell
+
+#ifdef ISO
+      REAL xt_seri(ntiso,klon,klev),xtmemoire(ntiso,klon,klev)
+      REAL d_xt_ajs(ntiso,klon,klev)
+      real d_xt_the(ntiso,klon,klev)
+#ifdef DIAGISO
+      real q_the(klon,klev)
+      real xt_the(ntiso,klon,klev)
+#endif
+      real qprec(klon,klev)
+      integer ixt
+#endif
+
+
+      u_seri(:,:)=u_seri_(:,:)
+      v_seri(:,:)=v_seri_(:,:)
+      t_seri(:,:)=t_seri_(:,:)
+      q_seri(:,:)=q_seri_(:,:)
+
+
+!  Modele du thermique
+!  ===================
+!         print*,'thermiques: WARNING on passe t au lieu de t_seri'
+
+
+! On prend comme valeur initiale des thermiques la valeur du pas
+! de temps precedent
+         zfm_therm(:,:)=fm_therm(:,:)
+         zdetr_therm(:,:)=detr_therm(:,:)
+         zentr_therm(:,:)=entr_therm(:,:)
+
+! On reinitialise les flux de masse a zero pour le cumul en
+! cas de splitting
+         fm_therm(:,:)=0.
+         entr_therm(:,:)=0.
+         detr_therm(:,:)=0.
+
+         ale_bl(:)=0.
+         alp_bl(:)=0.
+         if (prt_level.ge.10) then
+          print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
+         endif
+
+!   tests sur les valeurs negatives de l'eau
+         logexpr0=prt_level.ge.10
+         nbptspb=0
+         do k=1,klev
+            do i=1,klon
+! Attention teste abderr 19-03-09
+!               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
+                logexpr2(i,k)=.not.q_seri(i,k).ge.1.e-15
+               if (logexpr2(i,k)) then
+#ifdef ISO
+                qprec(i,k)=q_seri(i,k)
+#endif
+                q_seri(i,k)=1.e-15
+                nbptspb=nbptspb+1
+#ifdef ISO
+                do ixt=1,ntiso
+                  xt_seri(ixt,i,k)=1.e-15*(xt_seri(ixt,i,k)/qprec(i,k))
+                  ! xt_seri(ixt,i,k)=1.e-15*(Rdefault(index_iso(ixt)))
+                enddo
+#endif
+               endif
+!               if (logexpr0) &
+!    &             print*,'WARN eau<0 avant therm i=',i,'  k=',k  &
+!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k)
+            enddo
+         enddo
+         if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
+
+
+         new_thermcell=iflag_thermals>=15.and.iflag_thermals<=18
+#ifdef ISO
+      if (.not.new_thermcell) then
+           CALL abort_physic('calltherm 234','isos pas prevus ici',1)
+      endif
+#ifdef ISOVERIF
+      if (iso_eau.gt.0) then
+       call iso_verif_egalite_vect2D( &
+     &           xt_seri,q_seri, &
+     &           'calltherm 174',ntiso,klon,klev) 
+      endif !if (iso_eau.gt.0) then
+#endif    
+#endif
+         zdt=dtime/REAL(nsplit_thermals)
+
+
+         do isplit=1,nsplit_thermals
+
+          if (iflag_thermals>=1000) then
+            CALL thermcell_2002(klon,klev,zdt,iflag_thermals   &
+     &      ,pplay,paprs,pphi  &
+     &      ,u_seri,v_seri,t_seri,q_seri  &
+     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
+     &      ,zfm_therm,zentr_therm,fraca,zw2  &
+     &      ,r_aspect_thermals,30.,w2di_thermals  &
+     &      ,tau_thermals)
+          else if (iflag_thermals.eq.2) then
+            CALL thermcell_sec(klon,klev,zdt  &
+     &      ,pplay,paprs,pphi,zlev  &
+     &      ,u_seri,v_seri,t_seri,q_seri  &
+     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
+     &      ,zfm_therm,zentr_therm  &
+     &      ,r_aspect_thermals,30.,w2di_thermals  &
+     &      ,tau_thermals)
+          else if (iflag_thermals.eq.3) then
+            CALL thermcell_old(klon,klev,zdt  &
+     &      ,pplay,paprs,pphi  &
+     &      ,u_seri,v_seri,t_seri,q_seri  &
+     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
+     &      ,zfm_therm,zentr_therm  &
+     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
+     &      ,tau_thermals)
+          else if (iflag_thermals.eq.10) then
+            CALL thermcell_eau(klon,klev,zdt  &
+     &      ,pplay,paprs,pphi  &
+     &      ,u_seri,v_seri,t_seri,q_seri  &
+     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
+     &      ,zfm_therm,zentr_therm  &
+     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
+     &      ,tau_thermals)
+          else if (iflag_thermals.eq.11) then
+              abort_message = 'cas non prevu dans calltherm'
+              CALL abort_physic (modname,abort_message,1)
+          else if (iflag_thermals.eq.12) then
+            CALL calcul_sec(klon,klev,zdt  &
+     &      ,pplay,paprs,pphi,zlev  &
+     &      ,u_seri,v_seri,t_seri,q_seri  &
+     &      ,zmax_sec,wmax_sec,zw_sec,lmix_sec  &
+     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
+     &      ,tau_thermals)
+          else if (iflag_thermals==13.or.iflag_thermals==14) then
+              abort_message = 'thermcellV0_main enleve svn>2084'
+              CALL abort_physic (modname,abort_message,1)
+          else if (new_thermcell) then
+            CALL thermcell_main(itap,klon,klev,zdt  &
+     &      ,pplay,paprs,pphi,debut  &
+     &      ,u_seri,v_seri,t_seri,q_seri,t_env,q_env  &
+     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
+     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
+     &      ,ratqscth,ratqsdiff,zqsatth  &
+     &      ,zmax0,f0,zw2,fraca,ztv,zpspsk &
+     &      ,ztla,zthl,ztva &
+     &      ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,zcong &
+#ifdef ISO         
+     &      ,xt_seri,d_xt_the &
+#endif         
+     &   ) 
+
+            CALL thermcell_alp(klon,klev,zdt  &                      ! in
+     &        ,pplay,paprs  &                                        ! in
+     &        ,zfm_therm,zentr_therm,lmax  &                         ! in
+     &        ,pbl_tke,pctsrf,omega,airephy &                        ! in
+     &        ,zw2,fraca &                                           ! in
+     &        ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &  ! in
+     &        ,zcong,ale,alp,lalim_conv,wght_th &                          ! out
+     &        ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &! out
+     &        ,n2,s2,strig,ale_bl_stat &                                   ! out
+     &        ,therm_tke_max,env_tke_max &                           ! out
+     &        ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &          ! out
+     &        ,alp_bl_conv,alp_bl_stat &                             ! out
+     &        )
+
+           if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
+         else
+           abort_message = 'Cas des thermiques non prevu'
+           CALL abort_physic (modname,abort_message,1)
+         endif
+
+! Attention : les noms sont contre intuitif.
+! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille.
+! Il aurait mieux valu avoir un nobidouille_stratocu
+! Et pour simplifier :
+! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15)
+! Ce serait bien de changer, mai en prenant le temps de vérifier que ca
+! fait bien ce qu'on croit.
+
+       flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18
+
+! Calcul a posteriori du niveau max des thermiques pour les schémas qui
+! ne la sortent pas.
+      if (iflag_thermals<=12.or.iflag_thermals>=1000) then
+         lmax(:)=1
+         do k=1,klev-1
+            zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
+         enddo
+         do k=1,klev-1
+            do i=1,klon
+               if (zfm_therm(i,k+1)>0.) lmax(i)=k
+            enddo
+         enddo
+      endif
+
+      fact(:)=0.
+      DO i=1,klon
+       logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
+       IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
+      ENDDO
+
+     DO k=1,klev
+!  transformation de la derivee en tendance
+            d_t_the(:,k)=d_t_the(:,k)*dtime*fact(:)
+            d_u_the(:,k)=d_u_the(:,k)*dtime*fact(:)
+            d_v_the(:,k)=d_v_the(:,k)*dtime*fact(:)
+            d_q_the(:,k)=d_q_the(:,k)*dtime*fact(:)
+            fm_therm(:,k)=fm_therm(:,k)  &
+     &      +zfm_therm(:,k)*fact(:)
+            entr_therm(:,k)=entr_therm(:,k)  &
+     &       +zentr_therm(:,k)*fact(:)
+            detr_therm(:,k)=detr_therm(:,k)  &
+     &       +zdetr_therm(:,k)*fact(:)
+#ifdef ISO
+            do ixt=1,ntiso
+              d_xt_the(ixt,:,k)=d_xt_the(ixt,:,k)*dtime*fact(:)
+            enddo
+#endif
+      ENDDO
+       fm_therm(:,klev+1)=0.
+
+
+
+!  accumulation de la tendance
+            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
+            d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
+            d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
+            d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
+#ifdef ISO
+            d_xt_ajs(:,:,:)=d_xt_ajs(:,:,:)+d_xt_the(:,:,:)
+#endif
+
+!  incrementation des variables meteo
+            t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
+            u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
+            v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
+            qmemoire(:,:)=q_seri(:,:)
+            q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
+#ifdef ISO
+            xtmemoire(:,:,:)=xt_seri(:,:,:)
+            xt_seri(:,:,:) = xt_seri(:,:,:) + d_xt_the(:,:,:)
+#ifdef ISOVERIF
+!      write(*,*) 'calltherm 350 tmp: ajout d_xt_the'
+      if (iso_HDO.gt.0) then
+!      i=479
+!      k=4
+!      write(*,*) 'xt_seri(iso_hdo,i,k),q_seri(i,k)=', &
+!     &   xt_seri(iso_hdo,i,k),q_seri(i,k)
+!      write(*,*) 'd_xt_the(iso_hdo,i,k),d_q_the(i,k)=', &
+!     &   d_xt_the(iso_hdo,i,k),d_q_the(i,k)
+      call iso_verif_aberrant_enc_vect2D( &
+     &        xt_seri,q_seri, &
+     &        'calltherm 353, apres ajout d_xt_the',ntiso,klon,klev)
+      endif      
+#endif
+#endif 
+           if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
+
+       DO i=1,klon
+            fm_therm(i,klev+1)=0.
+            ale_bl(i)=ale_bl(i)+ale(i)/REAL(nsplit_thermals)
+!            write(22,*)'ALE CALLTHERM',ale_bl(i),ale(i)
+            alp_bl(i)=alp_bl(i)+alp(i)/REAL(nsplit_thermals)
+!            write(23,*)'ALP CALLTHERM',alp_bl(i),alp(i)
+        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)
+       ENDDO
+
+!IM 060508 marche pas comme cela !!!        enddo ! isplit
+
+!   tests sur les valeurs negatives de l'eau
+         nbptspb=0
+            DO k = 1, klev
+            DO i = 1, klon
+               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
+               if (logexpr2(i,k)) then
+                q_seri(i,k)=1.e-15
+                nbptspb=nbptspb+1
+#ifdef ISO
+                do ixt=1,ntiso
+                  xt_seri(ixt,i,k)=1.e-15*(xtmemoire(ixt,i,k)/qmemoire(i,k))
+                enddo
+#endif
+!                if (prt_level.ge.10) then
+!                  print*,'WARN eau<0 apres therm i=',i,'  k=',k  &
+!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k),  &
+!    &         'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
+                 endif
+            ENDDO
+            ENDDO
+#ifdef ISO
+#ifdef ISOVERIF
+      if (iso_HDO.gt.0) then
+      call iso_verif_aberrant_enc_vect2D( &
+     &        xt_seri,q_seri, &
+     &        'calltherm 393, apres bidouille q<0',ntiso,klon,klev)
+      endif      
+#endif
+#endif
+
+        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
+! tests sur les valeurs de la temperature
+        nbptspb=0
+            DO k = 1, klev
+            DO i = 1, klon
+               logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
+               if (logexpr2(i,k)) nbptspb=nbptspb+1
+!              if ((t_seri(i,k).lt.50.) .or.  &
+!    &              (t_seri(i,k).gt.370.)) then
+!                 print*,'WARN temp apres therm i=',i,'  k=',k  &
+!    &         ,' t_seri',t_seri(i,k)
+!              CALL abort
+!              endif
+            ENDDO
+            ENDDO
+        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
+         enddo ! isplit
+
+!
+!***************************************************************
+!     calcul du flux ascencant conservatif
+!            print*,'<<<<calcul flux ascendant conservatif'
+
+      fmc_therm=0.
+               do k=1,klev
+            do i=1,klon
+                  if (entr_therm(i,k).gt.0.) then
+                     fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
+                  else
+                     fmc_therm(i,k+1)=fmc_therm(i,k)
+                  endif
+                  detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1))  &
+     &                 -(fmc_therm(i,k)-fm_therm(i,k))
+               enddo
+            enddo
+      
+     
+!****************************************************************
+!     calcul de l'humidite dans l'ascendance
+!      print*,'<<<<calcul de lhumidite dans thermique'
+!CR:on ne le calcule que pour le cas sec
+      if (iflag_thermals.le.11) then      
+      do i=1,klon
+         zqasc(i,1)=q_seri(i,1)
+         do k=2,klev
+            if (fmc_therm(i,k+1).gt.1.e-6) then
+               zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1)  &
+     &              +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
+!CR:test on asseche le thermique
+!               zqasc(i,k)=zqasc(i,k)/2.
+!            else
+!               zqasc(i,k)=q_seri(i,k)
+            endif
+         enddo
+       enddo
+      
+
+!     calcul de l'eau condensee dans l'ascendance
+!             print*,'<<<<calcul de leau condensee dans thermique'
+             do i=1,klon
+                do k=1,klev
+                   clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
+                   if (clwcon0(i,k).lt.0. .or.   &
+     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
+                      clwcon0(i,k)=0.
+                   endif
+                enddo
+             enddo
+       else
+              do i=1,klon
+                do k=1,klev
+                   clwcon0(i,k)=zqla(i,k)  
+                   if (clwcon0(i,k).lt.0. .or.   &
+     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
+                   clwcon0(i,k)=0. 
+                   endif
+                enddo
+             enddo
+       endif
+!*******************************************************************    
+
+
+!jyg  Protection contre les temperatures nulles
+          do i=1,klon
+             do k=1,klev
+                if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
+             enddo
+          enddo
+
+
+      return
+
+      end subroutine calltherm
+
+END MODULE calltherm_mod
