Index: DZ6/trunk/libf/phylmdiso/calltherm.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/calltherm.F90	(revision 4090)
+++ 	(revision )
@@ -1,539 +1,0 @@
-!
-! $Id$
-!
-      subroutine calltherm(dtime  &
-     &      ,pplay,paprs,pphi,weak_inversion  &
-     &      ,u_seri,v_seri,t_seri,q_seri,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,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
-#ifdef ISO
-      use infotrac_phy, ONLY: ntraciso
-#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
-
-      implicit none
-      include "clesphys.h"
-      include "thermcell_old.h"
-
-
-!IM 140508
-      INTEGER, SAVE ::  itap
-!$OMP THREADPRIVATE(itap)
-      REAL dtime
-      LOGICAL debut
-      LOGICAL logexpr0, logexpr2(klon,klev), logexpr1(klon)
-      REAL fact(klon)
-      INTEGER nbptspb
-
-      REAL u_seri(klon,klev),v_seri(klon,klev)
-      REAL t_seri(klon,klev),q_seri(klon,klev),qmemoire(klon,klev)
-      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 zqta(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 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)
-      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) :: modname='calltherm'
-      character (len=80) :: abort_message
-
-      integer i,k,isplit
-      logical, save :: first=.true.
-      logical :: new_thermcell
-
-#ifdef ISO
-      REAL xt_seri(ntraciso,klon,klev),xtmemoire(ntraciso,klon,klev)
-      REAL d_xt_ajs(ntraciso,klon,klev)
-      real d_xt_the(ntraciso,klon,klev)
-#ifdef DIAGISO
-      real q_the(klon,klev)
-      real xt_the(ntraciso,klon,klev)
-#endif
-      real qprec(klon,klev)
-      integer ixt
-#endif
-
-
-!$OMP THREADPRIVATE(first)
-!********************************************************
-      if (first) then
-        itap=0
-        first=.false.
-      endif
-
-! Incrementer le compteur de la physique
-     itap   = itap + 1
-
-!  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,ntraciso
-                  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_gcm('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',ntraciso,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(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  &
-     &      ,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 &
-#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
-     &        ,ale,alp,lalim_conv,wght_th &                          ! out
-     &        ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &! out
-     &        ,n2,s2,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,ntraciso
-              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',ntraciso,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,ntraciso
-                  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',ntraciso,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
Index: DZ6/trunk/libf/phylmdiso/phys_state_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phys_state_var_mod.F90	(revision 4090)
+++ 	(revision )
@@ -1,878 +1,0 @@
-!
-! $Id: phys_state_var_mod.F90 4088 2022-03-10 07:03:20Z fhourdin $
-!
-      MODULE phys_state_var_mod
-! Variables sauvegardees pour le startphy.nc
-!======================================================================
-!
-!
-!======================================================================
-! Declaration des variables
-      USE dimphy
-      USE netcdf, only: nf90_fill_real
-      INTEGER, PARAMETER :: nlevSTD=17
-      INTEGER, PARAMETER :: nlevSTD8=8
-      INTEGER, PARAMETER :: nlevSTD3=3
-      INTEGER, PARAMETER :: nout=3
-      INTEGER, PARAMETER :: napisccp=1
-      INTEGER, SAVE :: radpas  ! radiation is called every "radpas" step
-      INTEGER, SAVE :: cvpas   ! convection is called every "cvpas" step
-      INTEGER, SAVE :: cvpas_0 = 1 ! reference value for cvpas
-      INTEGER, SAVE :: wkpas   ! wake scheme is called every "wkpas" step
-      REAL, PARAMETER :: missing_val_nf90=nf90_fill_real
-!$OMP THREADPRIVATE(radpas)
-!$OMP THREADPRIVATE(cvpas)
-!$OMP THREADPRIVATE(cvpas_0)
-!$OMP THREADPRIVATE(wkpas)
-      REAL, SAVE :: phys_tstep=0, solaire_etat0
-!$OMP THREADPRIVATE(phys_tstep, solaire_etat0)
-
-      REAL, ALLOCATABLE, SAVE :: pctsrf(:,:)
-!$OMP THREADPRIVATE(pctsrf)
-      REAL, ALLOCATABLE, SAVE :: ftsol(:,:)
-!$OMP THREADPRIVATE(ftsol)
-      REAL, ALLOCATABLE, SAVE :: beta_aridity(:,:)
-!$OMP THREADPRIVATE(beta_aridity)
-      REAL,ALLOCATABLE,SAVE :: qsol(:),fevap(:,:),z0m(:,:),z0h(:,:),agesno(:,:)
-!$OMP THREADPRIVATE(qsol,fevap,z0m,z0h,agesno)
-!FC drag des arbres
-      REAL, ALLOCATABLE, SAVE :: treedrg(:,:,:)
-!$OMP THREADPRIVATE(treedrg)
-
-!      character(len=6), SAVE :: ocean
-!!!!!!$OMP THREADPRIVATE(ocean)
-!      logical, SAVE :: ok_veget 
-!!!!!!$OMP THREADPRIVATE(ok_veget)
-      REAL, ALLOCATABLE, SAVE :: falb1(:,:), falb2(:,:)
-!$OMP THREADPRIVATE(falb1, falb2)
-
-!albedo SB >>>
-      REAL, ALLOCATABLE, SAVE :: falb_dif(:,:,:), falb_dir(:,:,:)
-      REAL, ALLOCATABLE, SAVE :: chl_con(:)
-!$OMP THREADPRIVATE(falb_dir,falb_dif,chl_con)
-!albedo SB <<<
-
-
-      REAL, ALLOCATABLE, SAVE :: rain_fall(:), snow_fall(:)
-!$OMP THREADPRIVATE( rain_fall, snow_fall)
-      REAL, ALLOCATABLE, SAVE :: solsw(:), solswfdiff(:), sollw(:)
-!$OMP THREADPRIVATE(solsw, solswfdiff, sollw)
-      REAL, ALLOCATABLE, SAVE :: radsol(:)
-!$OMP THREADPRIVATE(radsol)
-      REAL, ALLOCATABLE, SAVE :: swradcorr(:)
-!$OMP THREADPRIVATE(swradcorr)
-#ifdef ISO
-      REAL,ALLOCATABLE,SAVE :: xtsol(:,:),fxtevap(:,:,:)
-!$OMP THREADPRIVATE(xtsol,fxtevap)
-      REAL, ALLOCATABLE, SAVE :: xtrain_fall(:,:), xtsnow_fall(:,:)
-!$OMP THREADPRIVATE(xtrain_fall,xtsnow_fall)
-#endif
-
-!clesphy0 param physiq
-!
-! Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
-!
-      REAL, ALLOCATABLE, SAVE :: zmea(:), zstd(:), zsig(:), zgam(:)
-!$OMP THREADPRIVATE(zmea, zstd, zsig, zgam)
-      REAL, ALLOCATABLE, SAVE :: zthe(:), zpic(:), zval(:)
-!$OMP THREADPRIVATE(zthe, zpic, zval)
-!     REAL tabcntr0(100)
-      REAL, ALLOCATABLE, SAVE :: rugoro(:)
-!$OMP THREADPRIVATE(rugoro)
-      REAL, ALLOCATABLE, SAVE :: t_ancien(:,:), q_ancien(:,:)
-!$OMP THREADPRIVATE(t_ancien, q_ancien)
-      REAL, ALLOCATABLE, SAVE :: ql_ancien(:,:), qs_ancien(:,:)
-!$OMP THREADPRIVATE(ql_ancien, qs_ancien)
-      REAL, ALLOCATABLE, SAVE :: prw_ancien(:), prlw_ancien(:), prsw_ancien(:)
-!$OMP THREADPRIVATE(prw_ancien, prlw_ancien, prsw_ancien)
-#ifdef ISO
-      REAL, ALLOCATABLE, SAVE :: xt_ancien(:,:,:),xtl_ancien(:,:,:),xts_ancien(:,:,:)
-!$OMP THREADPRIVATE(xt_ancien,xtl_ancien,xts_ancien)
-#endif
-      REAL, ALLOCATABLE, SAVE :: u_ancien(:,:), v_ancien(:,:)
-!$OMP THREADPRIVATE(u_ancien, v_ancien)
-!!! RomP >>>
-      REAL, ALLOCATABLE, SAVE :: tr_ancien(:,:,:)
-!$OMP THREADPRIVATE(tr_ancien)
-!!! RomP <<<
-      LOGICAL, SAVE :: ancien_ok
-!$OMP THREADPRIVATE(ancien_ok)
-      REAL, ALLOCATABLE, SAVE :: clwcon(:,:),rnebcon(:,:)
-!$OMP THREADPRIVATE(clwcon,rnebcon)
-      REAL, ALLOCATABLE, SAVE :: rneb_ancien(:,:)
-!$OMP THREADPRIVATE(rneb_ancien)
-      REAL, ALLOCATABLE, SAVE :: qtc_cv(:,:),sigt_cv(:,:)
-!$OMP THREADPRIVATE(qtc_cv,sigt_cv)
-      REAL, ALLOCATABLE, SAVE :: ratqs(:,:)
-!$OMP THREADPRIVATE(ratqs)
-      REAL, ALLOCATABLE, SAVE :: pbl_tke(:,:,:) ! turb kinetic energy
-      REAL, ALLOCATABLE, SAVE :: coefh(:,:,:) ! Kz enthalpie
-      REAL, ALLOCATABLE, SAVE :: coefm(:,:,:) ! Kz momentum
-!$OMP THREADPRIVATE(pbl_tke, coefh,coefm)
-      REAL, ALLOCATABLE, SAVE :: zmax0(:), f0(:) ! 
-!$OMP THREADPRIVATE(zmax0,f0)
-      REAL, ALLOCATABLE, SAVE :: sig1(:,:), w01(:,:)
-!$OMP THREADPRIVATE(sig1,w01)
-      REAL, ALLOCATABLE, SAVE :: entr_therm(:,:), fm_therm(:,:)
-!$OMP THREADPRIVATE(entr_therm,fm_therm)
-      REAL, ALLOCATABLE, SAVE :: detr_therm(:,:)
-!$OMP THREADPRIVATE(detr_therm)
-!IM 150408
-!     pour phsystoke avec thermiques
-      REAL,ALLOCATABLE,SAVE :: clwcon0th(:,:),rnebcon0th(:,:)
-!$OMP THREADPRIVATE(clwcon0th,rnebcon0th)
-! radiation outputs
-      REAL,ALLOCATABLE,SAVE :: swdnc0(:,:), swdn0(:,:), swdn(:,:)
-!$OMP THREADPRIVATE(swdnc0,swdn0,swdn)
-      REAL,ALLOCATABLE,SAVE :: swupc0(:,:), swup0(:,:), swup(:,:)
-!$OMP THREADPRIVATE(swupc0, swup0,swup)
-      REAL,ALLOCATABLE,SAVE :: SWdn200clr(:), SWdn200(:)
-!$OMP THREADPRIVATE(SWdn200clr,SWdn200)
-      REAL,ALLOCATABLE,SAVE :: SWup200clr(:), SWup200(:)
-!$OMP THREADPRIVATE(SWup200clr,SWup200)
-      REAL,ALLOCATABLE,SAVE :: lwdnc0(:,:), lwdn0(:,:), lwdn(:,:)
-!$OMP THREADPRIVATE(lwdnc0,lwdn0,lwdn)
-      REAL,ALLOCATABLE,SAVE :: lwupc0(:,:), lwup0(:,:), lwup(:,:)
-!$OMP THREADPRIVATE(lwupc0,lwup0,lwup)
-      REAL,ALLOCATABLE,SAVE :: LWdn200clr(:), LWdn200(:)
-!$OMP THREADPRIVATE(LWdn200clr,LWdn200)
-      REAL,ALLOCATABLE,SAVE :: LWup200clr(:), LWup200(:)
-!$OMP THREADPRIVATE(LWup200clr,LWup200)
-      REAL,ALLOCATABLE,SAVE :: LWdnTOA(:), LWdnTOAclr(:)
-!$OMP THREADPRIVATE(LWdnTOA,LWdnTOAclr)
-! pressure level
-      REAL,ALLOCATABLE,SAVE :: tsumSTD(:,:,:)
-!$OMP THREADPRIVATE(tsumSTD)
-      REAL,ALLOCATABLE,SAVE :: usumSTD(:,:,:), vsumSTD(:,:,:)
-!$OMP THREADPRIVATE(usumSTD,vsumSTD)
-      REAL,ALLOCATABLE,SAVE :: wsumSTD(:,:,:), phisumSTD(:,:,:)
-!$OMP THREADPRIVATE(wsumSTD,phisumSTD)
-      REAL,ALLOCATABLE,SAVE :: qsumSTD(:,:,:), rhsumSTD(:,:,:)
-!$OMP THREADPRIVATE(qsumSTD,rhsumSTD)
-      REAL,ALLOCATABLE,SAVE :: tnondef(:,:,:) 
-!$OMP THREADPRIVATE(tnondef)
-      REAL,ALLOCATABLE,SAVE :: uvsumSTD(:,:,:)
-!$OMP THREADPRIVATE(uvsumSTD)
-      REAL,ALLOCATABLE,SAVE :: vqsumSTD(:,:,:)
-!$OMP THREADPRIVATE(vqsumSTD)
-      REAL,ALLOCATABLE,SAVE :: vTsumSTD(:,:,:)
-!$OMP THREADPRIVATE(vTsumSTD)
-      REAL,ALLOCATABLE,SAVE :: wqsumSTD(:,:,:)
-!$OMP THREADPRIVATE(wqsumSTD)
-      REAL,ALLOCATABLE,SAVE :: vphisumSTD(:,:,:)
-!$OMP THREADPRIVATE(vphisumSTD)
-      REAL,ALLOCATABLE,SAVE :: wTsumSTD(:,:,:)
-!$OMP THREADPRIVATE(wTsumSTD)
-      REAL,ALLOCATABLE,SAVE :: u2sumSTD(:,:,:)
-!$OMP THREADPRIVATE(u2sumSTD)
-      REAL,ALLOCATABLE,SAVE :: v2sumSTD(:,:,:)
-!$OMP THREADPRIVATE(v2sumSTD)
-      REAL,ALLOCATABLE,SAVE :: T2sumSTD(:,:,:)
-!$OMP THREADPRIVATE(T2sumSTD)
-      REAL,ALLOCATABLE,SAVE :: O3sumSTD(:,:,:), O3daysumSTD(:,:,:)
-!$OMP THREADPRIVATE(O3sumSTD,O3daysumSTD) 
-!IM begin
-      REAL,ALLOCATABLE,SAVE :: wlevSTD(:,:), ulevSTD(:,:), vlevSTD(:,:)
-!$OMP THREADPRIVATE(wlevSTD,ulevSTD,vlevSTD)
-      REAL,ALLOCATABLE,SAVE :: tlevSTD(:,:), qlevSTD(:,:), rhlevSTD(:,:)
-!$OMP THREADPRIVATE(tlevSTD,qlevSTD,rhlevSTD)
-      REAL,ALLOCATABLE,SAVE :: philevSTD(:,:)
-!$OMP THREADPRIVATE(philevSTD)
-      REAL,ALLOCATABLE,SAVE :: uvSTD(:,:)
-!$OMP THREADPRIVATE(uvSTD)
-      REAL,ALLOCATABLE,SAVE :: vqSTD(:,:)
-!$OMP THREADPRIVATE(vqSTD)
-      REAL,ALLOCATABLE,SAVE :: vTSTD(:,:)
-!$OMP THREADPRIVATE(vTSTD)
-      REAL,ALLOCATABLE,SAVE :: wqSTD(:,:)
-!$OMP THREADPRIVATE(wqSTD)
-      REAL,ALLOCATABLE,SAVE :: vphiSTD(:,:)
-!$OMP THREADPRIVATE(vphiSTD)
-      REAL,ALLOCATABLE,SAVE :: wTSTD(:,:)
-!$OMP THREADPRIVATE(wTSTD)
-      REAL,ALLOCATABLE,SAVE :: u2STD(:,:)
-!$OMP THREADPRIVATE(u2STD)
-      REAL,ALLOCATABLE,SAVE :: v2STD(:,:) 
-!$OMP THREADPRIVATE(v2STD)
-      REAL,ALLOCATABLE,SAVE :: T2STD(:,:)
-!$OMP THREADPRIVATE(T2STD)
-      REAL,ALLOCATABLE,SAVE :: O3STD(:,:), O3daySTD(:,:)
-!$OMP THREADPRIVATE(O3STD,O3daySTD)
-!IM end
-      INTEGER,ALLOCATABLE,SAVE :: seed_old(:,:)
-!$OMP THREADPRIVATE(seed_old)
-      REAL,ALLOCATABLE,SAVE :: zuthe(:),zvthe(:)
-!$OMP THREADPRIVATE(zuthe,zvthe)
-      REAL,ALLOCATABLE,SAVE :: alb_neig(:)
-!$OMP THREADPRIVATE(alb_neig)
-!cloud base mass flux
-      REAL,ALLOCATABLE,SAVE :: ema_cbmf(:)
-!$OMP THREADPRIVATE(ema_cbmf)
-!cloud base pressure & cloud top pressure
-      REAL,ALLOCATABLE,SAVE :: ema_pcb(:), ema_pct(:)
-!$OMP THREADPRIVATE(ema_pcb,ema_pct)
-      REAL,ALLOCATABLE,SAVE :: Mipsh(:,:)     ! mass flux shed from  adiab. ascents
-!$OMP THREADPRIVATE(Mipsh)
-      REAL,ALLOCATABLE,SAVE :: Ma(:,:)       ! undilute upward mass flux
-!$OMP THREADPRIVATE(Ma)
-      REAL,ALLOCATABLE,SAVE :: qcondc(:,:)    ! in-cld water content from convect
-!$OMP THREADPRIVATE(qcondc)
-      REAL,ALLOCATABLE,SAVE :: wd(:) ! sb
-!$OMP THREADPRIVATE(wd)
-      REAL,ALLOCATABLE,SAVE :: sigd(:)
-!$OMP THREADPRIVATE(sigd)
-!
-      REAL,ALLOCATABLE,SAVE :: cin(:)
-!$OMP THREADPRIVATE(cin)
-! ftd : convective heating due to unsaturated downdraughts
-      REAL,ALLOCATABLE,SAVE :: ftd(:,:)
-!$OMP THREADPRIVATE(ftd)
-! fqd : convective moistening due to unsaturated downdraughts
-      REAL,ALLOCATABLE,SAVE :: fqd(:,:)     
-!$OMP THREADPRIVATE(fqd)
-#ifdef ISO
-      REAL, ALLOCATABLE, SAVE :: fxtd(:,:,:)
-!$OMP THREADPRIVATE(fxtd)
-#endif
-!34EK
-! -- Variables de controle de ALE et ALP
-!ALE : Energie disponible pour soulevement : utilisee par la 
-!      convection d'Emanuel pour le declenchement et la regulation
-      REAL,ALLOCATABLE,SAVE :: ALE(:)
-!$OMP THREADPRIVATE(ALE)
-!ALP : Puissance  disponible pour soulevement
-      REAL,ALLOCATABLE,SAVE :: ALP(:)
-!$OMP THREADPRIVATE(ALP)
-!
-! nouvelles variables pour le couplage convection-couche limite
-      REAL,ALLOCATABLE,SAVE :: Ale_bl(:)
-!$OMP THREADPRIVATE(Ale_bl)
-      REAL,ALLOCATABLE,SAVE :: Alp_bl(:)
-!$OMP THREADPRIVATE(Alp_bl)
-      INTEGER,ALLOCATABLE,SAVE :: lalim_conv(:)
-!$OMP THREADPRIVATE(lalim_conv)
-      REAL,ALLOCATABLE,SAVE :: wght_th(:,:)
-!$OMP THREADPRIVATE(wght_th)
-      REAL,ALLOCATABLE,SAVE    :: ale_wake(:)
-!$OMP THREADPRIVATE(ale_wake)
-      REAL,ALLOCATABLE,SAVE    :: ale_bl_stat(:)
-!$OMP THREADPRIVATE(ale_bl_stat)
-!
-! variables de la wake
-! wake_deltat : ecart de temperature avec la zone non perturbee
-! wake_deltaq : ecart d'humidite avec la zone non perturbee
-! wake_s      : fraction surfacique occupee par la poche froide
-! awake_dens  : number of active wakes per unit area
-! wake_dens   : number of wakes per unit area
-! cv_gen      : birth rate of cumulonimbus per unit area.
-! wake_occ    : occurence of wakes (= 1 if wakes occur, =0 otherwise)
-! wake_Cstar  : vitesse d'etalement de la poche
-! wake_pe     : wake potential energy - WAPE
-! wake_fip    : Gust Front Impinging power - ALP
-      REAL,ALLOCATABLE,SAVE :: wake_deltat(:,:)
-!$OMP THREADPRIVATE(wake_deltat)
-      REAL,ALLOCATABLE,SAVE :: wake_deltaq(:,:)
-!$OMP THREADPRIVATE(wake_deltaq)
-#ifdef ISO
-      REAL, ALLOCATABLE, SAVE :: wake_deltaxt(:,:,:)
-!$OMP THREADPRIVATE(wake_deltaxt)
-#endif
-      REAL,ALLOCATABLE,SAVE :: wake_s(:)
-!$OMP THREADPRIVATE(wake_s)
-      REAL,ALLOCATABLE,SAVE :: awake_dens(:), wake_dens(:)
-!$OMP THREADPRIVATE(awake_dens, wake_dens)
-      REAL,ALLOCATABLE,SAVE :: cv_gen(:)
-!$OMP THREADPRIVATE(cv_gen)
-      REAL,ALLOCATABLE,SAVE :: wake_Cstar(:)
-!$OMP THREADPRIVATE(wake_Cstar)
-      REAL,ALLOCATABLE,SAVE :: wake_pe(:)
-!$OMP THREADPRIVATE(wake_pe)
-      REAL,ALLOCATABLE,SAVE :: wake_fip(:)
-!$OMP THREADPRIVATE(wake_fip)
-!
-!jyg<
-! variables related to the spitting of the PBL between wake and 
-! off-wake regions.
-! wake_delta_pbl_TKE : difference TKE_w - TKE_x
-      REAL,ALLOCATABLE,SAVE :: wake_delta_pbl_TKE(:,:,:)
-!$OMP THREADPRIVATE(wake_delta_pbl_TKE)
-!nrlmd<
-      REAL, ALLOCATABLE, SAVE :: delta_tsurf(:,:) ! Surface temperature difference inside-outside cold pool
-!$OMP THREADPRIVATE(delta_tsurf)
-!>nrlmd
-!>jyg
-!
-! pfrac_impa : Produits des coefs lessivage impaction
-! pfrac_nucl : Produits des coefs lessivage nucleation
-! pfrac_1nucl: Produits des coefs lessi nucl (alpha = 1) 
-      REAL,ALLOCATABLE,SAVE :: pfrac_impa(:,:), pfrac_nucl(:,:)
-!$OMP THREADPRIVATE(pfrac_impa,pfrac_nucl)
-      REAL,ALLOCATABLE,SAVE :: pfrac_1nucl(:,:)
-!$OMP THREADPRIVATE(pfrac_1nucl)
-!
-      REAL,ALLOCATABLE,SAVE :: total_rain(:), nday_rain(:)  
-!$OMP THREADPRIVATE(total_rain,nday_rain)
-      REAL,ALLOCATABLE,SAVE :: paire_ter(:)
-!$OMP THREADPRIVATE(paire_ter)
-! albsol1: albedo du sol total pour SW visible
-! albsol2: albedo du sol total pour SW proche IR
-      REAL,ALLOCATABLE,SAVE :: albsol1(:), albsol2(:)
-!$OMP THREADPRIVATE(albsol1,albsol2)
-
-!albedo SB >>>
-      REAL,ALLOCATABLE,SAVE :: albsol_dif(:,:),albsol_dir(:,:)
-!$OMP THREADPRIVATE(albsol_dif,albsol_dir)
-!albedo SB <<<
-
-
-      REAL, ALLOCATABLE, SAVE:: wo(:, :, :)
-      ! column-density of ozone in a layer, in kilo-Dobsons
-      ! Third dimension has size 1 or 2.
-      ! "wo(:, :, 1)" is for the average day-night field, 
-      ! "wo(:, :, 2)" is for daylight time.
-      !$OMP THREADPRIVATE(wo)
-
-! heat : chauffage solaire
-! heat0: chauffage solaire ciel clair
-! cool : refroidissement infrarouge
-! cool0 : refroidissement infrarouge ciel clair
-! sollwdown : downward LW flux at surface
-! sollwdownclr : downward CS LW flux at surface
-! toplwdown : downward CS LW flux at TOA
-! toplwdownclr : downward CS LW flux at TOA
-! heat_volc : chauffage solaire du au volcanisme
-! cool_volc : refroidissement infrarouge du au volcanisme
-      REAL,ALLOCATABLE,SAVE :: clwcon0(:,:),rnebcon0(:,:)
-!$OMP THREADPRIVATE(clwcon0,rnebcon0)
-      REAL,ALLOCATABLE,SAVE :: heat(:,:)   
-!$OMP THREADPRIVATE(heat)
-      REAL,ALLOCATABLE,SAVE :: heat0(:,:)
-!$OMP THREADPRIVATE(heat0)
-      REAL,ALLOCATABLE,SAVE :: cool(:,:)
-!$OMP THREADPRIVATE(cool)
-      REAL,ALLOCATABLE,SAVE :: cool0(:,:)
-!$OMP THREADPRIVATE(cool0)
-      REAL,ALLOCATABLE,SAVE :: heat_volc(:,:)   
-!$OMP THREADPRIVATE(heat_volc)
-      REAL,ALLOCATABLE,SAVE :: cool_volc(:,:)
-!$OMP THREADPRIVATE(cool_volc)
-      REAL,ALLOCATABLE,SAVE :: topsw(:), toplw(:)
-!$OMP THREADPRIVATE(topsw,toplw)
-      REAL,ALLOCATABLE,SAVE :: sollwdown(:)
-!$OMP THREADPRIVATE(sollwdown)
-      REAL,ALLOCATABLE,SAVE :: gustiness(:)
-!$OMP THREADPRIVATE(gustiness)
-      REAL,ALLOCATABLE,SAVE :: sollwdownclr(:)
-!$OMP THREADPRIVATE(sollwdownclr)
-      REAL,ALLOCATABLE,SAVE :: toplwdown(:)
-!$OMP THREADPRIVATE(toplwdown)
-      REAL,ALLOCATABLE,SAVE :: toplwdownclr(:)
-!$OMP THREADPRIVATE(toplwdownclr)
-      REAL,ALLOCATABLE,SAVE :: topsw0(:),toplw0(:),solsw0(:),sollw0(:)
-!$OMP THREADPRIVATE(topsw0,toplw0,solsw0,sollw0)
-      REAL,ALLOCATABLE,SAVE :: albpla(:)
-!$OMP THREADPRIVATE(albpla)
-
-!IM ajout variables CFMIP2/CMIP5
-      REAL,ALLOCATABLE,SAVE :: heatp(:,:), coolp(:,:)
-!$OMP THREADPRIVATE(heatp, coolp)
-      REAL,ALLOCATABLE,SAVE :: heat0p(:,:), cool0p(:,:)
-!$OMP THREADPRIVATE(heat0p, cool0p)
-      REAL,ALLOCATABLE,SAVE :: radsolp(:), topswp(:), toplwp(:)
-!$OMP THREADPRIVATE(radsolp, topswp, toplwp)
-      REAL,ALLOCATABLE,SAVE :: albplap(:)
-!$OMP THREADPRIVATE(albplap)
-      REAL,ALLOCATABLE,SAVE :: solswp(:), solswfdiffp(:), sollwp(:)
-!$OMP THREADPRIVATE(solswp, solswfdiffp, sollwp)
-      REAL,ALLOCATABLE,SAVE :: sollwdownp(:)
-!$OMP THREADPRIVATE(sollwdownp)
-      REAL,ALLOCATABLE,SAVE :: topsw0p(:),toplw0p(:)
-      REAL,ALLOCATABLE,SAVE :: solsw0p(:),sollw0p(:)
-!$OMP THREADPRIVATE(topsw0p,toplw0p,solsw0p,sollw0p)
-      REAL,ALLOCATABLE,SAVE :: lwdnc0p(:,:), lwdn0p(:,:), lwdnp(:,:)
-      REAL,ALLOCATABLE,SAVE :: lwupc0p(:,:), lwup0p(:,:), lwupp(:,:)
-!$OMP THREADPRIVATE(lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp)
-      REAL,ALLOCATABLE,SAVE :: swdnc0p(:,:), swdn0p(:,:), swdnp(:,:)
-      REAL,ALLOCATABLE,SAVE :: swupc0p(:,:), swup0p(:,:), swupp(:,:)
-!$OMP THREADPRIVATE(swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp)
-
-! pbase : cloud base pressure
-! bbase : cloud base buoyancy
-      REAL,ALLOCATABLE,SAVE :: cape(:)
-!$OMP THREADPRIVATE(cape)
-      REAL,ALLOCATABLE,SAVE :: pbase(:)
-!$OMP THREADPRIVATE(pbase)
-      REAL,ALLOCATABLE,SAVE :: bbase(:)
-!$OMP THREADPRIVATE(bbase)
-!
-      REAL,SAVE,ALLOCATABLE :: zqasc(:,:)
-!$OMP THREADPRIVATE( zqasc)
-      INTEGER,ALLOCATABLE,SAVE :: ibas_con(:), itop_con(:)
-!$OMP THREADPRIVATE(ibas_con,itop_con)
-      REAL,SAVE,ALLOCATABLE :: rain_con(:)
-!$OMP THREADPRIVATE(rain_con)
-      REAL,SAVE,ALLOCATABLE :: snow_con(:)
-!$OMP THREADPRIVATE(snow_con)
-!
-#ifdef ISO
-      REAL,SAVE,ALLOCATABLE :: xtrain_con(:,:)
-!$OMP THREADPRIVATE(xtrain_con)
-      REAL,SAVE,ALLOCATABLE :: xtsnow_con(:,:)
-!$OMP THREADPRIVATE(xtsnow_con)
-#endif
-      REAL,SAVE,ALLOCATABLE :: rlonPOS(:)
-!$OMP THREADPRIVATE(rlonPOS)
-      REAL,SAVE,ALLOCATABLE :: newsst(:)
-!$OMP THREADPRIVATE(newsst)
-      REAL,SAVE,ALLOCATABLE :: ustar(:,:),u10m(:,:), v10m(:,:),wstar(:,:)
-!$OMP THREADPRIVATE(ustar,u10m,v10m,wstar)
-!
-! ok_ade=T -ADE=topswad-topsw
-! ok_aie=T ->
-!       ok_ade=T -AIE=topswai-topswad
-!       ok_ade=F -AIE=topswai-topsw
-!
-!topswad, solswad : Aerosol direct effect
-      REAL,SAVE,ALLOCATABLE :: topswad(:), solswad(:)
-!$OMP THREADPRIVATE(topswad,solswad)
-!topswai, solswai : Aerosol indirect effect
-      REAL,SAVE,ALLOCATABLE :: topswai(:), solswai(:)
-!$OMP THREADPRIVATE(topswai,solswai)
-
-      REAL,SAVE,ALLOCATABLE :: tau_aero(:,:,:,:), piz_aero(:,:,:,:), cg_aero(:,:,:,:)
-!$OMP THREADPRIVATE(tau_aero, piz_aero, cg_aero)
-      REAL,SAVE,ALLOCATABLE :: tau_aero_sw_rrtm(:,:,:,:), piz_aero_sw_rrtm(:,:,:,:), cg_aero_sw_rrtm(:,:,:,:)
-!$OMP THREADPRIVATE(tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm)
-      REAL,SAVE,ALLOCATABLE :: tau_aero_lw_rrtm(:,:,:,:), piz_aero_lw_rrtm(:,:,:,:), cg_aero_lw_rrtm(:,:,:,:)
-!$OMP THREADPRIVATE(tau_aero_lw_rrtm, piz_aero_lw_rrtm, cg_aero_lw_rrtm)
-      REAL,SAVE,ALLOCATABLE :: ccm(:,:,:)
-!$OMP THREADPRIVATE(ccm)
-
-      REAL,SAVE,ALLOCATABLE :: ale_bl_trig(:)
-!$OMP THREADPRIVATE(ale_bl_trig)
-
-      REAL,SAVE,ALLOCATABLE :: ratqs_inter(:,:)
-!$OMP THREADPRIVATE(ratqs_inter)
-
-#ifdef ISO
-#ifdef ISOTRAC
-      INTEGER,SAVE,ALLOCATABLE :: bassin_map(:)
-!$OMP THREADPRIVATE(bassin_map)
-      INTEGER,SAVE,ALLOCATABLE :: boite_map(:,:)
-!$OMP THREADPRIVATE(boite_map)
-#endif   
-#endif
-      REAL, ALLOCATABLE, SAVE:: du_gwd_rando(:, :), du_gwd_front(:, :)
-      !$OMP THREADPRIVATE(du_gwd_rando, du_gwd_front)
-      ! tendencies on wind due to gravity waves
-
-      LOGICAL,SAVE :: is_initialized=.FALSE.
-!$OMP THREADPRIVATE(is_initialized)    
-
-      ! Ocean-atmosphere interface:
-
-      REAL, ALLOCATABLE, SAVE:: ds_ns(:) ! (klon)
-      ! "delta salinity near surface". Salinity variation in the
-      ! near-surface turbulent layer. That is subskin salinity minus
-      ! foundation salinity. In ppt.
-
-      REAL, ALLOCATABLE, SAVE:: dt_ns(:) ! (klon)
-      ! "delta temperature near surface". Temperature variation in the
-      ! near-surface turbulent layer. That is subskin temperature
-      ! minus foundation temperature. (Can be negative.) In K.
-      
-      REAL, ALLOCATABLE, SAVE:: delta_sst(:) ! (klon)
-      ! Ocean-air interface temperature minus bulk SST, in
-      ! K. Allocated and defined only if activate_ocean_skin >= 1.
-
-      REAL, ALLOCATABLE, SAVE:: delta_sal(:) ! (klon)
-      ! Ocean-air interface salinity minus bulk salinity, in ppt
-      
-      !$OMP THREADPRIVATE(delta_sal, ds_ns, dt_ns, delta_sst)
-
-    CONTAINS
-
-!======================================================================
-SUBROUTINE phys_state_var_init(read_climoz)
-USE dimphy
-USE aero_mod
-USE infotrac_phy, ONLY : nbtr
-#ifdef ISO
-USE infotrac_phy, ONLY : ntraciso,niso
-#endif
-USE indice_sol_mod
-use config_ocean_skin_m, only: activate_ocean_skin
-IMPLICIT NONE
-
-integer, intent(in)::  read_climoz
-! read ozone climatology
-! Allowed values are 0, 1 and 2
-! 0: do not read an ozone climatology
-! 1: read a single ozone climatology that will be used day and night
-! 2: read two ozone climatologies, the average day and night
-! climatology and the daylight climatology
-
-include "clesphys.h"
-
-      print*, 'is_initialized', is_initialized
-      IF (is_initialized) RETURN
-      is_initialized=.TRUE.
-      ALLOCATE(pctsrf(klon,nbsrf))
-      ALLOCATE(ftsol(klon,nbsrf))
-      ALLOCATE(beta_aridity(klon,nbsrf))
-      ALLOCATE(qsol(klon),fevap(klon,nbsrf))
-      ALLOCATE(z0m(klon,nbsrf+1),z0h(klon,nbsrf+1),agesno(klon,nbsrf))
-!FC
-      ALLOCATE(treedrg(klon,klev,nbsrf))
-      ALLOCATE(falb1(klon,nbsrf))
-      ALLOCATE(falb2(klon,nbsrf))
-!albedo SB >>> 
-      print*, 'allocate falb'
-      ALLOCATE(falb_dir(klon,nsw,nbsrf),falb_dif(klon,nsw,nbsrf))
-!!      print*, 'allocate falb good', falb_dir(1,1,1)
-      ALLOCATE(chl_con(klon))
-!albedo SB <<<
-      ALLOCATE(rain_fall(klon))
-      ALLOCATE(snow_fall(klon))
-      ALLOCATE(solsw(klon), solswfdiff(klon), sollw(klon))
-      sollw=0.0
-      ALLOCATE(radsol(klon))
-      ALLOCATE(swradcorr(klon))
-      ALLOCATE(zmea(klon), zstd(klon), zsig(klon), zgam(klon))
-      ALLOCATE(zthe(klon), zpic(klon), zval(klon))
-
-      ALLOCATE(rugoro(klon))
-      ALLOCATE(t_ancien(klon,klev), q_ancien(klon,klev))
-      ALLOCATE(ql_ancien(klon,klev), qs_ancien(klon,klev))
-      ALLOCATE(prw_ancien(klon), prlw_ancien(klon), prsw_ancien(klon))
-      ALLOCATE(u_ancien(klon,klev), v_ancien(klon,klev))
-!!! Rom P >>>
-      ALLOCATE(tr_ancien(klon,klev,nbtr))
-!!! Rom P <<<
-      ALLOCATE(clwcon(klon,klev),rnebcon(klon,klev))
-      ALLOCATE(rneb_ancien(klon,klev))
-      ALLOCATE(qtc_cv(klon,klev),sigt_cv(klon,klev))
-      ALLOCATE(ratqs(klon,klev))
-      ALLOCATE(pbl_tke(klon,klev+1,nbsrf+1))
-!nrlmd<
-      ALLOCATE(delta_tsurf(klon,nbsrf))
-!>nrlmd
-      ALLOCATE(coefh(klon,klev+1,nbsrf+1))
-      ALLOCATE(coefm(klon,klev+1,nbsrf+1))
-      ALLOCATE(zmax0(klon), f0(klon))
-      ALLOCATE(sig1(klon,klev), w01(klon,klev))
-      ALLOCATE(entr_therm(klon,klev), fm_therm(klon,klev+1))
-      ALLOCATE(detr_therm(klon,klev))
-!     pour phsystoke avec thermiques
-      ALLOCATE(clwcon0th(klon,klev),rnebcon0th(klon,klev))
-! radiation outputs
-      ALLOCATE(swdnc0(klon,klevp1), swdn0(klon,klevp1), swdn(klon,klevp1))
-      ALLOCATE(swupc0(klon,klevp1), swup0(klon,klevp1), swup(klon,klevp1))
-      ALLOCATE(lwdnc0(klon,klevp1), lwdn0(klon,klevp1), lwdn(klon,klevp1))
-      ALLOCATE(lwupc0(klon,klevp1), lwup0(klon,klevp1), lwup(klon,klevp1))
-      ALLOCATE(SWdn200clr(klon), SWdn200(klon))
-      ALLOCATE(SWup200clr(klon), SWup200(klon))
-      ALLOCATE(LWdn200clr(klon), LWdn200(klon))
-      ALLOCATE(LWup200clr(klon), LWup200(klon))
-      ALLOCATE(LWdnTOA(klon), LWdnTOAclr(klon))
-! pressure level
-      ALLOCATE(tsumSTD(klon,nlevSTD,nout))
-      ALLOCATE(usumSTD(klon,nlevSTD,nout), vsumSTD(klon,nlevSTD,nout))
-      ALLOCATE(wsumSTD(klon,nlevSTD,nout), phisumSTD(klon,nlevSTD,nout))
-      ALLOCATE(qsumSTD(klon,nlevSTD,nout), rhsumSTD(klon,nlevSTD,nout))
-      ALLOCATE(tnondef(klon,nlevSTD,nout))
-      ALLOCATE(uvsumSTD(klon,nlevSTD,nout))
-      ALLOCATE(vqsumSTD(klon,nlevSTD,nout))
-      ALLOCATE(vTsumSTD(klon,nlevSTD,nout))
-      ALLOCATE(wqsumSTD(klon,nlevSTD,nout))
-      ALLOCATE(vphisumSTD(klon,nlevSTD,nout))
-      ALLOCATE(wTsumSTD(klon,nlevSTD,nout))
-      ALLOCATE(u2sumSTD(klon,nlevSTD,nout))
-      ALLOCATE(v2sumSTD(klon,nlevSTD,nout))
-      ALLOCATE(T2sumSTD(klon,nlevSTD,nout))
-      ALLOCATE(O3sumSTD(klon,nlevSTD,nout))
-      ALLOCATE(O3daysumSTD(klon,nlevSTD,nout))
-!IM beg
-      ALLOCATE(wlevSTD(klon,nlevSTD), ulevSTD(klon,nlevSTD), vlevSTD(klon,nlevSTD))
-      ALLOCATE(tlevSTD(klon,nlevSTD), qlevSTD(klon,nlevSTD), rhlevSTD(klon,nlevSTD))
-      ALLOCATE(philevSTD(klon,nlevSTD))
-      ALLOCATE(uvSTD(klon,nlevSTD),vqSTD(klon,nlevSTD))
-      ALLOCATE(vTSTD(klon,nlevSTD),wqSTD(klon,nlevSTD))
-      ALLOCATE(vphiSTD(klon,nlevSTD),wTSTD(klon,nlevSTD))
-      ALLOCATE(u2STD(klon,nlevSTD),v2STD(klon,nlevSTD))
-      ALLOCATE(T2STD(klon,nlevSTD))
-      ALLOCATE(O3STD(klon,nlevSTD))
-      ALLOCATE(O3daySTD(klon,nlevSTD))
-!IM end
-      ALLOCATE(seed_old(klon,napisccp))
-      ALLOCATE(zuthe(klon),zvthe(klon))
-      ALLOCATE(alb_neig(klon))
-!cloud base mass flux
-      ALLOCATE(ema_cbmf(klon))
-!cloud base pressure & cloud top pressure
-      ALLOCATE(ema_pcb(klon), ema_pct(klon))
-!
-      ALLOCATE(Mipsh(klon,klev))
-      ALLOCATE(Ma(klon,klev))
-      ALLOCATE(qcondc(klon,klev))
-      ALLOCATE(wd(klon))
-      ALLOCATE(sigd(klon))
-      ALLOCATE(cin(klon), ALE(klon), ALP(klon))
-      ALLOCATE(ftd(klon,klev), fqd(klon,klev))
-      ALLOCATE(Ale_bl(klon))
-      ALLOCATE(ale_wake(klon))
-      ALLOCATE(ale_bl_stat(klon))
-      ale_bl_stat(:)=0
-      ALLOCATE(Alp_bl(klon))
-      ALLOCATE(lalim_conv(klon))
-      ALLOCATE(wght_th(klon,klev))
-      ALLOCATE(wake_deltat(klon,klev), wake_deltaq(klon,klev))
-      ALLOCATE(wake_s(klon), awake_dens(klon), wake_dens(klon))
-!!      awake_dens = 0.  ! initialized in phyetat0
-      ALLOCATE(cv_gen(klon))
-      ALLOCATE(wake_Cstar(klon))
-      ALLOCATE(wake_pe(klon), wake_fip(klon))
-!jyg<
-      ALLOCATE(wake_delta_pbl_TKE(klon,klev+1,nbsrf+1))
-!>jyg
-      ALLOCATE(pfrac_impa(klon,klev), pfrac_nucl(klon,klev))
-      ALLOCATE(pfrac_1nucl(klon,klev))
-      ALLOCATE(total_rain(klon), nday_rain(klon))
-      ALLOCATE(paire_ter(klon))
-      ALLOCATE(albsol1(klon), albsol2(klon))
-!albedo SB >>>
-      ALLOCATE(albsol_dir(klon,nsw),albsol_dif(klon,nsw))
-!albedo SB <<<
-
-      if (read_climoz <= 1) then
-         ALLOCATE(wo(klon,klev, 1))
-      else
-         ! read_climoz == 2
-         ALLOCATE(wo(klon,klev, 2))
-      end if
-      
-      ALLOCATE(clwcon0(klon,klev),rnebcon0(klon,klev))
-      ALLOCATE(heat(klon,klev), heat0(klon,klev)) 
-      ALLOCATE(cool(klon,klev), cool0(klon,klev))
-      ALLOCATE(heat_volc(klon,klev), cool_volc(klon,klev)) 
-      ALLOCATE(topsw(klon), toplw(klon))
-      ALLOCATE(sollwdown(klon), sollwdownclr(klon))
-      sollwdown = 0.
-      ALLOCATE(toplwdown(klon), toplwdownclr(klon))
-      ALLOCATE(topsw0(klon),toplw0(klon),solsw0(klon),sollw0(klon))
-      sollw0 = 0.
-      ALLOCATE(albpla(klon))
-!IM ajout variables CFMIP2/CMIP5
-      ALLOCATE(heatp(klon,klev), coolp(klon,klev))
-      ALLOCATE(heat0p(klon,klev), cool0p(klon,klev))
-      ALLOCATE(radsolp(klon), topswp(klon), toplwp(klon))
-      ALLOCATE(albplap(klon))
-      ALLOCATE(solswp(klon), solswfdiffp(klon), sollwp(klon))
-      ALLOCATE(gustiness(klon))
-      ALLOCATE(sollwdownp(klon))
-      ALLOCATE(topsw0p(klon),toplw0p(klon))
-      ALLOCATE(solsw0p(klon),sollw0p(klon))
-      ALLOCATE(lwdnc0p(klon,klevp1), lwdn0p(klon,klevp1), lwdnp(klon,klevp1))
-      ALLOCATE(lwupc0p(klon,klevp1), lwup0p(klon,klevp1), lwupp(klon,klevp1))
-      ALLOCATE(swdnc0p(klon,klevp1), swdn0p(klon,klevp1), swdnp(klon,klevp1))
-      ALLOCATE(swupc0p(klon,klevp1), swup0p(klon,klevp1), swupp(klon,klevp1))
-
-      ALLOCATE(cape(klon))
-      ALLOCATE(pbase(klon),bbase(klon))
-      ALLOCATE(zqasc(klon,klev))
-      ALLOCATE(ibas_con(klon), itop_con(klon))
-      ALLOCATE(rain_con(klon), snow_con(klon))
-      ALLOCATE(rlonPOS(klon))
-      ALLOCATE(newsst(klon))
-      ALLOCATE(ustar(klon,nbsrf),u10m(klon,nbsrf), v10m(klon,nbsrf),wstar(klon,nbsrf+1))
-      ALLOCATE(topswad(klon), solswad(klon))
-      ALLOCATE(topswai(klon), solswai(klon))
-      ALLOCATE(tau_aero(klon,klev,naero_grp,nbands),piz_aero(klon,klev,naero_grp,nbands),cg_aero(klon,klev,naero_grp,nbands))
-      ALLOCATE(tau_aero_sw_rrtm(klon,klev,2,nbands_sw_rrtm),piz_aero_sw_rrtm(klon,klev,2,nbands_sw_rrtm))
-      ALLOCATE(cg_aero_sw_rrtm(klon,klev,2,nbands_sw_rrtm))
-      ALLOCATE(tau_aero_lw_rrtm(klon,klev,2,nbands_lw_rrtm),piz_aero_lw_rrtm(klon,klev,2,nbands_lw_rrtm))
-      ALLOCATE(cg_aero_lw_rrtm(klon,klev,2,nbands_lw_rrtm))
-      ALLOCATE(ccm(klon,klev,nbands))
-
-#ifdef ISO
-      ALLOCATE(xtsol(niso,klon),fxtevap(ntraciso,klon,nbsrf))
-      ALLOCATE(fxtd(ntraciso,klon,klev))
-      ALLOCATE(wake_deltaxt(ntraciso,klon,klev))
-      ALLOCATE(xt_ancien(ntraciso,klon,klev))
-      ALLOCATE(xtl_ancien(ntraciso,klon,klev))
-      ALLOCATE(xts_ancien(ntraciso,klon,klev))
-      ALLOCATE(xtrain_fall(ntraciso,klon))
-      ALLOCATE(xtsnow_fall(ntraciso,klon))
-      ALLOCATE(xtrain_con(ntraciso,klon))
-      ALLOCATE(xtsnow_con(ntraciso,klon))
-#ifdef ISOTRAC
-      ALLOCATE(bassin_map(klon))
-      ALLOCATE(boite_map(klon,klev))  
-#endif      
-#endif
-
-      ALLOCATE(ale_bl_trig(klon))
-      ALLOCATE(ratqs_inter(klon,klev))
-      IF (ok_gwd_rando) THEN
-        ALLOCATE(du_gwd_rando(klon, klev))
-        du_gwd_rando(:,:)=0.
-      ENDIF
-      IF (.not. ok_hines .and. ok_gwd_rando) THEN
-        ALLOCATE(du_gwd_front(klon, klev))
-        du_gwd_front(:,:) = 0 !ym missing init   
-      ENDIF
-      if (activate_ocean_skin >= 1) ALLOCATE(delta_sal(klon), ds_ns(klon), &
-           dt_ns(klon), delta_sst(klon))
-
-    END SUBROUTINE phys_state_var_init
-
-!======================================================================
-    SUBROUTINE phys_state_var_end
-      ! Useful only for lmdz1d.
-!USE dimphy
-USE indice_sol_mod
-use config_ocean_skin_m, only: activate_ocean_skin
-IMPLICIT NONE
-include "clesphys.h"
-
-      DEALLOCATE(pctsrf, ftsol, falb1, falb2)
-      DEALLOCATE(beta_aridity)
-      DEALLOCATE(qsol,fevap,z0m,z0h,agesno)
-!FC
-      DEALLOCATE(treedrg)
-      DEALLOCATE(rain_fall, snow_fall, solsw, solswfdiff, sollw, radsol, swradcorr)
-      DEALLOCATE(zmea, zstd, zsig, zgam)
-      DEALLOCATE(zthe, zpic, zval)
-      DEALLOCATE(rugoro, t_ancien, q_ancien, clwcon, rnebcon)
-      DEALLOCATE(qs_ancien, ql_ancien, rneb_ancien)
-      DEALLOCATE(prw_ancien, prlw_ancien, prsw_ancien)
-      DEALLOCATE(qtc_cv,sigt_cv)
-      DEALLOCATE(u_ancien, v_ancien)
-      DEALLOCATE(tr_ancien)                           !RomP
-      DEALLOCATE(ratqs, pbl_tke,coefh,coefm)
-      DEALLOCATE(zmax0, f0)
-      DEALLOCATE(sig1, w01)
-      DEALLOCATE(entr_therm, fm_therm)
-      DEALLOCATE(detr_therm)
-      DEALLOCATE(clwcon0th, rnebcon0th)
-! radiation outputs
-      DEALLOCATE(swdnc0, swdn0, swdn)
-      DEALLOCATE(swupc0, swup0, swup)
-      DEALLOCATE(lwdnc0, lwdn0, lwdn)
-      DEALLOCATE(lwupc0, lwup0, lwup)
-      DEALLOCATE(SWdn200clr, SWdn200)
-      DEALLOCATE(SWup200clr, SWup200)
-      DEALLOCATE(LWdn200clr, LWdn200)
-      DEALLOCATE(LWup200clr, LWup200)
-      DEALLOCATE(LWdnTOA, LWdnTOAclr)
-! pressure level
-      DEALLOCATE(tsumSTD)
-      DEALLOCATE(usumSTD, vsumSTD)
-      DEALLOCATE(wsumSTD, phisumSTD)
-      DEALLOCATE(tnondef)
-      DEALLOCATE(qsumSTD, rhsumSTD)
-      DEALLOCATE(uvsumSTD)
-      DEALLOCATE(vqsumSTD)
-      DEALLOCATE(vTsumSTD)
-      DEALLOCATE(wqsumSTD)
-      DEALLOCATE(vphisumSTD)
-      DEALLOCATE(wTsumSTD)
-      DEALLOCATE(u2sumSTD)
-      DEALLOCATE(v2sumSTD)
-      DEALLOCATE(T2sumSTD)
-      DEALLOCATE(O3sumSTD)
-      DEALLOCATE(O3daysumSTD)
-!IM beg
-      DEALLOCATE(wlevSTD,ulevSTD,vlevSTD,tlevSTD,qlevSTD,rhlevSTD,philevSTD)
-      DEALLOCATE(uvSTD,vqSTD,vTSTD,wqSTD,vphiSTD,wTSTD,u2STD,v2STD,T2STD,O3STD,O3daySTD)
-!IM end
-      DEALLOCATE(seed_old)
-      DEALLOCATE(zuthe, zvthe)
-      DEALLOCATE(alb_neig)
-      DEALLOCATE(ema_cbmf)
-      DEALLOCATE(ema_pcb, ema_pct)
-      DEALLOCATE(Mipsh, Ma, qcondc)
-      DEALLOCATE(wd, sigd)
-      DEALLOCATE(cin, ALE, ALP)
-      DEALLOCATE(ftd, fqd)
-      DEALLOCATE(Ale_bl, Alp_bl)
-      DEALLOCATE(ale_wake)
-      DEALLOCATE(ale_bl_stat)
-      DEALLOCATE(lalim_conv, wght_th)
-      DEALLOCATE(wake_deltat, wake_deltaq)
-      DEALLOCATE(wake_s, awake_dens, wake_dens)
-      DEALLOCATE(cv_gen)
-      DEALLOCATE(wake_Cstar, wake_pe, wake_fip)
-!jyg<
-      DEALLOCATE(wake_delta_pbl_TKE)
-!nrlmd<
-      DEALLOCATE(delta_tsurf)
-!>nrlmd
-!>jyg
-      DEALLOCATE(pfrac_impa, pfrac_nucl)
-      DEALLOCATE(pfrac_1nucl)
-      DEALLOCATE(total_rain, nday_rain)
-      DEALLOCATE(paire_ter)
-      DEALLOCATE(albsol1, albsol2)
-!albedo SB >>>
-      DEALLOCATE(albsol_dir,albsol_dif,falb_dir,falb_dif,chl_con)
-!albedo SB <<<
-      DEALLOCATE(wo)
-      DEALLOCATE(clwcon0,rnebcon0)
-      DEALLOCATE(heat, heat0) 
-      DEALLOCATE(cool, cool0)
-      DEALLOCATE(heat_volc, cool_volc) 
-      DEALLOCATE(topsw, toplw)
-      DEALLOCATE(sollwdown, sollwdownclr)
-      DEALLOCATE(gustiness)
-      DEALLOCATE(toplwdown, toplwdownclr)
-      DEALLOCATE(topsw0,toplw0,solsw0,sollw0)
-      DEALLOCATE(albpla)
-!IM ajout variables CFMIP2/CMIP5
-      DEALLOCATE(heatp, coolp)
-      DEALLOCATE(heat0p, cool0p)
-      DEALLOCATE(radsolp, topswp, toplwp)
-      DEALLOCATE(albplap)
-      DEALLOCATE(solswp, solswfdiffp, sollwp)
-      DEALLOCATE(sollwdownp)
-      DEALLOCATE(topsw0p,toplw0p)
-      DEALLOCATE(solsw0p,sollw0p)
-      DEALLOCATE(lwdnc0p, lwdn0p, lwdnp)
-      DEALLOCATE(lwupc0p, lwup0p, lwupp)
-      DEALLOCATE(swdnc0p, swdn0p, swdnp)
-      DEALLOCATE(swupc0p, swup0p, swupp)
-      DEALLOCATE(cape)
-      DEALLOCATE(pbase,bbase)
-      DEALLOCATE(zqasc)
-      DEALLOCATE(ibas_con, itop_con)
-      DEALLOCATE(rain_con, snow_con)
-      DEALLOCATE(rlonPOS)
-      DEALLOCATE(newsst)
-      DEALLOCATE(ustar,u10m, v10m,wstar)
-      DEALLOCATE(topswad, solswad)
-      DEALLOCATE(topswai, solswai)
-      DEALLOCATE(tau_aero,piz_aero,cg_aero)
-      DEALLOCATE(tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
-      DEALLOCATE(tau_aero_lw_rrtm,piz_aero_lw_rrtm,cg_aero_lw_rrtm)
-      DEALLOCATE(ccm)
-      if (ok_gwd_rando) DEALLOCATE(du_gwd_rando)
-      if (.not. ok_hines .and. ok_gwd_rando) DEALLOCATE(du_gwd_front)
-      DEALLOCATE(ale_bl_trig)
-      DEALLOCATE(ratqs_inter)
-
-      if (activate_ocean_skin >= 1) deALLOCATE(delta_sal, ds_ns, dt_ns, &
-           delta_sst)
-
-#ifdef ISO    
-      DEALLOCATE(xtsol,fxtevap)  
-      DEALLOCATE(xt_ancien,xtl_ancien,xts_ancien, fxtd, wake_deltaxt)
-      DEALLOCATE(xtrain_fall, xtsnow_fall, xtrain_con, xtsnow_con)
-#ifdef ISOTRAC 
-      DEALLOCATE(bassin_map,boite_map) 
-#endif        
-#endif
-      is_initialized=.FALSE.
-      
-END SUBROUTINE phys_state_var_end
-
-      END MODULE phys_state_var_mod
Index: DZ6/trunk/libf/phylmdiso/thermcell_main.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/thermcell_main.F90	(revision 4090)
+++ 	(revision )
@@ -1,857 +1,0 @@
-
-! $Id: thermcell_main.F90 3451 2019-01-27 11:07:30Z fhourdin $
-!
-      subroutine thermcell_main(itap,ngrid,nlay,ptimestep  &
-     &                  ,pplay,pplev,pphi,debut  &
-     &                  ,pu,pv,pt,po  &
-     &                  ,pduadj,pdvadj,pdtadj,pdoadj  &
-     &                  ,fm0,entr0,detr0,zqta,zqla,lmax  &
-     &                  ,ratqscth,ratqsdiff,zqsatth  &
-     &                  ,zmax0, f0,zw2,fraca,ztv &
-     &                  ,zpspsk,ztla,zthl,ztva &
-     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
-#ifdef ISO         
-     &      ,xtpo,xtpdoadj &
-#endif         
-     &   )
-
-
-      USE thermcell_ini_mod, ONLY: thermcell_ini,dqimpl,dvdq,prt_level,lunout,prt_level
-      USE thermcell_ini_mod, ONLY: iflag_thermals_closure,iflag_thermals_ed,tau_thermals,r_aspect_thermals
-      USE thermcell_ini_mod, ONLY: RD,RG
-
-#ifdef ISO
-  USE infotrac_phy, ONLY : ntraciso
-#ifdef ISOVERIF
-  USE isotopes_mod, ONLY : iso_eau,iso_HDO
-  USE isotopes_verif_mod, ONLY: iso_verif_egalite, &
-        iso_verif_aberrant_encadre
-#endif
-#endif
-
-
-      IMPLICIT NONE
-
-!=======================================================================
-!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
-!   Version du 09.02.07
-!   Calcul du transport vertical dans la couche limite en presence
-!   de "thermiques" explicitement representes avec processus nuageux
-!
-!   Reecriture a partir d'un listing papier a Habas, le 14/02/00
-!
-!   le thermique est suppose homogene et dissipe par melange avec
-!   son environnement. la longueur l_mix controle l'efficacite du
-!   melange
-!
-!   Le calcul du transport des differentes especes se fait en prenant
-!   en compte:
-!     1. un flux de masse montant
-!     2. un flux de masse descendant
-!     3. un entrainement
-!     4. un detrainement
-!
-! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
-!    Introduction of an implicit computation of vertical advection in
-!    the environment of thermal plumes in thermcell_dq
-!    impl =     0 : explicit, 1 : implicit, -1 : old version
-!    controled by iflag_thermals =
-!       15, 16 run with impl=-1 : numerical convergence with NPv3
-!       17, 18 run with impl=1  : more stable
-!    15 and 17 correspond to the activation of the stratocumulus "bidouille"
-!
-!=======================================================================
-
-
-!-----------------------------------------------------------------------
-!   declarations:
-!   -------------
-
-
-!   arguments:
-!   ----------
-      integer, intent(in) :: itap,ngrid,nlay
-      real, intent(in) ::  ptimestep
-      real, intent(in), dimension(ngrid,nlay)    :: pt,pu,pv,po,pplay,pphi,zpspsk
-      real, intent(in), dimension(ngrid,nlay+1)  :: pplev
-      real, intent(out), dimension(ngrid,nlay)   :: pdtadj,pduadj,pdvadj,pdoadj,entr0,detr0
-      real, intent(out), dimension(ngrid,nlay)   :: ztla,zqla,zqta,zqsatth,zthl
-      real, intent(out), dimension(ngrid,nlay+1) :: fm0,zw2,fraca
-      real, intent(out), dimension(ngrid) :: zmax0,f0
-      real, intent(out), dimension(ngrid,nlay) :: ztva,ztv
-      logical, intent(in) :: debut
-
-      real, intent(out), dimension(ngrid) :: pcon
-      real, intent(out), dimension(ngrid,nlay) :: rhobarz,wth3
-      real, intent(out), dimension(ngrid) :: wmax_sec
-      integer,intent(out), dimension(ngrid) :: lalim
-      real, intent(out), dimension(ngrid,nlay+1) :: fm
-      real, intent(out), dimension(ngrid,nlay) :: alim_star
-      real, intent(out), dimension(ngrid) :: zmax
-
-!   local:
-!   ------
-
-
-      integer,save :: igout=1
-!$OMP THREADPRIVATE(igout)
-      integer,save :: lunout1=6
-!$OMP THREADPRIVATE(lunout1)
-      integer,save :: lev_out=10
-!$OMP THREADPRIVATE(lev_out)
-
-      real lambda, zf,zf2,var,vardiff,CHI
-      integer ig,k,l,ierr,ll
-      logical sorties
-      real, dimension(ngrid) :: linter,zmix, zmax_sec
-      integer,dimension(ngrid) :: lmax,lmin,lmix,lmix_bis,nivcon
-      real, dimension(ngrid,nlay) :: ztva_est
-      real, dimension(ngrid,nlay) :: deltaz,zlay,zh,zdthladj,zu,zv,zo,zl,zva,zua,zoa
-      real, dimension(ngrid,nlay) :: zta,zha,q2,wq,wthl,wthv,thetath2,wth2
-      real, dimension(ngrid,nlay) :: ratqscth,ratqsdiff,rho,masse
-      real, dimension(ngrid,nlay+1) :: zw_est,zlev
-      real, dimension(ngrid) :: wmax,wmax_tmp
-      real, dimension(ngrid,nlay+1) :: f_star
-      real, dimension(ngrid,nlay) :: entr,detr,entr_star,detr_star,alim_star_clos
-      real, dimension(ngrid,nlay) :: zqsat,csc
-      real, dimension(ngrid) :: zcon,zcon2,alim_star_tot,f
-
-      character (len=20) :: modname='thermcell_main'
-      character (len=80) :: abort_message
-
-
-#ifdef ISO
-      REAL xtpo(ntraciso,ngrid,nlay),xtpdoadj(ntraciso,ngrid,nlay)
-      REAL xtzo(ntraciso,ngrid,nlay)
-      REAL xtpdoadj_tmp(ngrid,nlay)
-      REAL xtpo_tmp(ngrid,nlay)
-      REAL xtzo_tmp(ngrid,nlay)
-      integer ixt
-#endif
-
-!
-
-!-----------------------------------------------------------------------
-!   initialisation:
-!   ---------------
-!
-    print*,'NEW THERMCELL cool'
-
-
-   fm=0. ; entr=0. ; detr=0.
-
-      if (prt_level.ge.1) print*,'thermcell_main V4'
-
-       sorties=.true.
-      IF(ngrid.NE.ngrid) THEN
-         PRINT*
-         PRINT*,'STOP dans convadj'
-         PRINT*,'ngrid    =',ngrid
-         PRINT*,'ngrid  =',ngrid
-      ENDIF
-!
-!     write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
-     do ig=1,ngrid
-         f0(ig)=max(f0(ig),1.e-2)
-         zmax0(ig)=max(zmax0(ig),40.)
-!IMmarche pas ?!       if (f0(ig)<1.e-2) f0(ig)=1.e-2
-     enddo
-
-      if (prt_level.ge.20) then
-       do ig=1,ngrid
-          print*,'th_main ig f0',ig,f0(ig)
-       enddo
-      endif
-!-----------------------------------------------------------------------
-! Calcul de T,q,ql a partir de Tl et qT dans l environnement
-!   --------------------------------------------------------------------
-!
-      CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
-     &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
-       
-      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
-
-!------------------------------------------------------------------------
-!                       --------------------
-!
-!
-!                       + + + + + + + + + + +
-!
-!
-!  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
-!  wh,wt,wo ...
-!
-!                       + + + + + + + + + + +  zh,zu,zv,zo,rho
-!
-!
-!                       --------------------   zlev(1)
-!                       \\\\\\\\\\\\\\\\\\\\
-!
-!
-
-!-----------------------------------------------------------------------
-!   Calcul des altitudes des couches
-!-----------------------------------------------------------------------
-
-      do l=2,nlay
-         zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
-      enddo
-      zlev(:,1)=0.
-      zlev(:,nlay+1)=(2.*pphi(:,nlay)-pphi(:,nlay-1))/RG
-      do l=1,nlay
-         zlay(:,l)=pphi(:,l)/RG
-      enddo
-      do l=1,nlay
-         deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
-      enddo
-
-!-----------------------------------------------------------------------
-!   Calcul des densites et masses
-!-----------------------------------------------------------------------
-
-      rho(:,:)=pplay(:,:)/(zpspsk(:,:)*RD*ztv(:,:))
-      if (prt_level.ge.10) write(lunout,*) 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
-      rhobarz(:,1)=rho(:,1)
-      do l=2,nlay
-         rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
-      enddo
-      do l=1,nlay
-         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
-      enddo
-      if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
-
-!------------------------------------------------------------------
-!
-!             /|\
-!    --------  |  F_k+1 -------   
-!                              ----> D_k
-!             /|\              <---- E_k , A_k
-!    --------  |  F_k --------- 
-!                              ----> D_k-1
-!                              <---- E_k-1 , A_k-1
-!
-!
-!
-!
-!
-!    ---------------------------
-!
-!    ----- F_lmax+1=0 ----------         \
-!            lmax     (zmax)              |
-!    ---------------------------          |
-!                                         |
-!    ---------------------------          |
-!                                         |
-!    ---------------------------          |
-!                                         |
-!    ---------------------------          |
-!                                         |
-!    ---------------------------          |
-!                                         |  E
-!    ---------------------------          |  D
-!                                         |
-!    ---------------------------          |
-!                                         |
-!    ---------------------------  \       |
-!            lalim                 |      |
-!    ---------------------------   |      |
-!                                  |      |
-!    ---------------------------   |      |
-!                                  | A    |
-!    ---------------------------   |      |
-!                                  |      |
-!    ---------------------------   |      |
-!    lmin  (=1 pour le moment)     |      |
-!    ----- F_lmin=0 ------------  /      /
-!
-!    ---------------------------
-!    //////////////////////////
-!
-!
-!=============================================================================
-!  Calculs initiaux ne faisant pas intervenir les changements de phase
-!=============================================================================
-
-!------------------------------------------------------------------
-!  1. alim_star est le profil vertical de l'alimentation a la base du
-!     panache thermique, calcule a partir de la flotabilite de l'air sec
-!  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
-!------------------------------------------------------------------
-!
-      entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
-      lmin=1
-
-!-----------------------------------------------------------------------------
-!  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
-!     panache sec conservatif (e=d=0) alimente selon alim_star 
-!     Il s'agit d'un calcul de type CAPE
-!     zmax_sec est utilise pour determiner la geometrie du thermique.
-!------------------------------------------------------------------------------
-!---------------------------------------------------------------------------------
-!calcul du melange et des variables dans le thermique
-!--------------------------------------------------------------------------------
-!
-      if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
-
-!=====================================================================
-! Old version of thermcell_plume in thermcell_plume_6A.F90
-! It includes both thermcell_plume_6A and thermcell_plume_5B corresponding
-! to the 5B and 6A versions used for CMIP5 and CMIP6.
-! The latest was previously named thermcellV1_plume.
-! The new thermcell_plume is a clean version (removing obsolete
-! options) of thermcell_plume_6A.
-! The 3 versions are controled by
-! flag_thermals_ed <= 9 thermcell_plume_6A
-!                  <= 19 thermcell_plume_5B
-!                  else thermcell_plume (default 20 for convergence with 6A)
-! Fredho
-!=====================================================================
-
-      if (iflag_thermals_ed<=9) then
-!         print*,'THERM NOUVELLE/NOUVELLE Arnaud'
-         CALL thermcell_plume_6A(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
-     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
-     &    lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
-     &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
-     &    ,lev_out,lunout1,igout)
-
-      elseif (iflag_thermals_ed<=19) then
-!        print*,'THERM RIO et al 2010, version d Arnaud'
-         CALL thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
-     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
-     &    lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
-     &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
-     &    ,lev_out,lunout1,igout)
-      else
-         CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
-     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
-     &    lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
-     &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
-     &    ,lev_out,lunout1,igout)
-      endif
-
-      if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
-
-      call test_ltherm(ngrid,nlay,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
-      call test_ltherm(ngrid,nlay,pplay,lmix ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
-
-      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
-      if (prt_level.ge.10) then
-         write(lunout1,*) 'Dans thermcell_main 2'
-         write(lunout1,*) 'lmin ',lmin(igout)
-         write(lunout1,*) 'lalim ',lalim(igout)
-         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
-         write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
-     &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
-      endif
-
-!-------------------------------------------------------------------------------
-! Calcul des caracteristiques du thermique:zmax,zmix,wmax
-!-------------------------------------------------------------------------------
-!
-      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,  &
-     &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
-! Attention, w2 est transforme en sa racine carree dans cette routine
-! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
-      wmax_tmp=0.
-      do  l=1,nlay
-         wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
-      enddo
-!     print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
-
-
-
-      call test_ltherm(ngrid,nlay,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
-      call test_ltherm(ngrid,nlay,pplay,lmin ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
-      call test_ltherm(ngrid,nlay,pplay,lmix ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
-      call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
-
-      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
-
-!-------------------------------------------------------------------------------
-! Fermeture,determination de f
-!-------------------------------------------------------------------------------
-!
-!
-!!      write(lunout,*)'THERM NOUVEAU XXXXX'
-      CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
-    &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
-
- 
-call test_ltherm(ngrid,nlay,pplay,lmin,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
-call test_ltherm(ngrid,nlay,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
-
-      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
-      if (prt_level.ge.10) then
-         write(lunout1,*) 'Dans thermcell_main 1b'
-         write(lunout1,*) 'lmin ',lmin(igout)
-         write(lunout1,*) 'lalim ',lalim(igout)
-         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
-         write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
-     &    ,l=1,lalim(igout)+4)
-      endif
-
-
-
-
-! Choix de la fonction d'alimentation utilisee pour la fermeture.
-! Apparemment sans importance
-      alim_star_clos(:,:)=alim_star(:,:)
-      alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
-!
-!CR Appel de la fermeture seche 
-      if (iflag_thermals_closure.eq.1) then
-
-      CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho,  &
-     &   zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! Appel avec les zmax et wmax tenant compte de la condensation
-! Semble moins bien marcher
-     else if (iflag_thermals_closure.eq.2) then
-
-     CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho,  &
-    &   zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
-
-     endif
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
-
-      if (tau_thermals>1.) then
-         lambda=exp(-ptimestep/tau_thermals)
-         f0=(1.-lambda)*f+lambda*f0
-      else
-         f0=f
-      endif
-
-! Test valable seulement en 1D mais pas genant
-      if (.not. (f0(1).ge.0.) ) then
-              abort_message = '.not. (f0(1).ge.0.)'
-              CALL abort_physic (modname,abort_message,1)
-      endif
-
-!-------------------------------------------------------------------------------
-!deduction des flux
-
-      CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
-     &       lalim,lmax,alim_star,  &
-     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
-     &       detr,zqla,lev_out,lunout1,igout)
-!IM 060508    &       detr,zqla,zmax,lev_out,lunout,igout)
-
-      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
-      call test_ltherm(ngrid,nlay,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
-      call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
-
-!------------------------------------------------------------------
-!   On ne prend pas directement les profils issus des calculs precedents
-!   mais on s'autorise genereusement une relaxation vers ceci avec
-!   une constante de temps tau_thermals (typiquement 1800s).
-!------------------------------------------------------------------
-
-      if (tau_thermals>1.) then
-         lambda=exp(-ptimestep/tau_thermals)
-         fm0=(1.-lambda)*fm+lambda*fm0
-         entr0=(1.-lambda)*entr+lambda*entr0
-         detr0=(1.-lambda)*detr+lambda*detr0
-      else
-         fm0=fm
-         entr0=entr
-         detr0=detr
-      endif
-
-!c------------------------------------------------------------------
-!   calcul du transport vertical
-!------------------------------------------------------------------
-
-      call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &
-     &                    zthl,zdthladj,zta,lev_out)
-      call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &
-     &                   po,pdoadj,zoa,lev_out)
-
-#ifdef ISO
-        ! C Risi: on utilise directement la même routine
-        do ixt=1,ntraciso
-          do ll=1,nlay
-            DO ig=1,ngrid
-                xtpo_tmp(ig,ll)=xtpo(ixt,ig,ll)
-                xtzo_tmp(ig,ll)=xtzo(ixt,ig,ll)
-            enddo
-          enddo
-          call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &
-     &                   xtpo_tmp,xtpdoadj_tmp,xtzo_tmp,lev_out)
-          do ll=1,nlay
-            DO ig=1,ngrid
-                xtpdoadj(ixt,ig,ll)=xtpdoadj_tmp(ig,ll)
-            enddo
-          enddo
-        enddo !do ixt=1,ntraciso
-#endif
-
-#ifdef ISO      
-#ifdef ISOVERIF
-      DO  ll=1,nlay
-        DO ig=1,ngrid
-          if (iso_eau.gt.0) then
-              call iso_verif_egalite(xtpo(iso_eau,ig,ll), &
-     &          po(ig,ll),'thermcell_main 594')
-              call iso_verif_egalite(xtpdoadj(iso_eau,ig,ll), &
-     &          pdoadj(ig,ll),'thermcell_main 596')
-          endif
-          if (iso_HDO.gt.0) then
-              call iso_verif_aberrant_encadre(xtpo(iso_hdo,ig,ll) &
-     &           /po(ig,ll),'thermcell_main 610')
-          endif
-        enddo
-      enddo !DO  ll=1,nlay 
-      write(*,*) 'thermcell_main 600 tmp: apres thermcell_dq'
-#endif      
-#endif
-
-
-
-!------------------------------------------------------------------
-! Calcul de la fraction de l'ascendance
-!------------------------------------------------------------------
-      do ig=1,ngrid
-         fraca(ig,1)=0.
-         fraca(ig,nlay+1)=0.
-      enddo
-      do l=2,nlay
-         do ig=1,ngrid
-            if (zw2(ig,l).gt.1.e-10) then
-            fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
-            else
-            fraca(ig,l)=0.
-            endif
-         enddo
-      enddo
-     
-!------------------------------------------------------------------
-!  calcul du transport vertical du moment horizontal
-!------------------------------------------------------------------
-
-!IM 090508  
-      if (dvdq == 0 ) then
-
-! Calcul du transport de V tenant compte d'echange par gradient
-! de pression horizontal avec l'environnement
-
-         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse  &
-!    &    ,fraca*dvdq,zmax &
-     &    ,fraca,zmax &
-     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
-
-      else
-
-! calcul purement conservatif pour le transport de V
-         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse  &
-     &    ,zu,pduadj,zua,lev_out)
-         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse  &
-     &    ,zv,pdvadj,zva,lev_out)
-
-      endif
-
-!     print*,'13 OK convect8'
-      do l=1,nlay
-         do ig=1,ngrid
-           pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)  
-         enddo
-      enddo
-
-      if (prt_level.ge.1) print*,'14 OK convect8'
-!------------------------------------------------------------------
-!   Calculs de diagnostiques pour les sorties
-!------------------------------------------------------------------
-!calcul de fraca pour les sorties
-      
-      if (sorties) then
-      if (prt_level.ge.1) print*,'14a OK convect8'
-! calcul du niveau de condensation
-! initialisation
-      do ig=1,ngrid
-         nivcon(ig)=0
-         zcon(ig)=0.
-      enddo 
-!nouveau calcul
-      do ig=1,ngrid
-      CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
-      pcon(ig)=pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI
-      enddo
-!IM   do k=1,nlay
-      do k=1,nlay-1
-         do ig=1,ngrid
-         if ((pcon(ig).le.pplay(ig,k))  &
-     &      .and.(pcon(ig).gt.pplay(ig,k+1))) then
-            zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
-         endif
-         enddo
-      enddo
-!IM
-      ierr=0
-      do ig=1,ngrid
-        if (pcon(ig).le.pplay(ig,nlay)) then 
-           zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
-           ierr=1
-        endif
-      enddo
-      if (ierr==1) then
-           abort_message = 'thermcellV0_main: les thermiques vont trop haut '
-           CALL abort_physic (modname,abort_message,1)
-      endif
-
-      if (prt_level.ge.1) print*,'14b OK convect8'
-      do k=nlay,1,-1
-         do ig=1,ngrid
-            if (zqla(ig,k).gt.1e-10) then
-               nivcon(ig)=k
-               zcon(ig)=zlev(ig,k)
-            endif
-         enddo
-      enddo
-      if (prt_level.ge.1) print*,'14c OK convect8'
-!calcul des moments
-!initialisation
-      do l=1,nlay
-         do ig=1,ngrid
-            q2(ig,l)=0.
-            wth2(ig,l)=0.
-            wth3(ig,l)=0.
-            ratqscth(ig,l)=0.
-            ratqsdiff(ig,l)=0.
-         enddo
-      enddo      
-      if (prt_level.ge.1) print*,'14d OK convect8'
-      if (prt_level.ge.10)write(lunout,*)                                &
-    &     'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
-      do l=1,nlay
-         do ig=1,ngrid
-            zf=fraca(ig,l)
-            zf2=zf/(1.-zf)
-!
-            thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2
-            if(zw2(ig,l).gt.1.e-10) then
-             wth2(ig,l)=zf2*(zw2(ig,l))**2
-            else
-             wth2(ig,l)=0.
-            endif
-            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
-     &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
-            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
-!test: on calcul q2/po=ratqsc
-            ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
-         enddo
-      enddo
-!calcul des flux: q, thetal et thetav
-      do l=1,nlay
-         do ig=1,ngrid
-      wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-po(ig,l)*1000.)
-      wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l))
-      wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l))
-         enddo
-      enddo
-
-!calcul du ratqscdiff
-      if (prt_level.ge.1) print*,'14e OK convect8'
-      var=0.
-      vardiff=0.
-      ratqsdiff(:,:)=0.
-
-      do l=1,nlay
-         do ig=1,ngrid
-            if (l<=lalim(ig)) then
-            var=var+alim_star(ig,l)*zqta(ig,l)*1000.
-            endif
-         enddo
-      enddo
-
-      if (prt_level.ge.1) print*,'14f OK convect8'
-
-      do l=1,nlay
-         do ig=1,ngrid
-            if (l<=lalim(ig)) then
-               zf=fraca(ig,l)
-               zf2=zf/(1.-zf)
-               vardiff=vardiff+alim_star(ig,l)*(zqta(ig,l)*1000.-var)**2
-            endif
-         enddo
-      enddo
-
-      if (prt_level.ge.1) print*,'14g OK convect8'
-         do l=1,nlay
-            do ig=1,ngrid
-               ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
-            enddo
-         enddo 
-      endif
-
-      if (prt_level.ge.1) print*,'thermcell_main FIN  OK'
-
-      return
-      end subroutine thermcell_main
-
-!=============================================================================
-!/////////////////////////////////////////////////////////////////////////////
-!=============================================================================
-      subroutine test_ltherm(ngrid,nlay,pplay,long,ztv,po,ztva, &  ! in
-    &            zqla,f_star,zw2,comment)                          ! in
-!=============================================================================
-      USE thermcell_ini_mod, ONLY: prt_level
-      IMPLICIT NONE
-
-      integer i, k, ngrid,nlay
-      real, intent(in), dimension(ngrid,nlay) :: pplay,ztv,po,ztva,zqla
-      real, intent(in), dimension(ngrid,nlay) :: f_star,zw2
-      integer, intent(in), dimension(ngrid) :: long
-      real seuil
-      character*21 comment
-      seuil=0.25
-
-      if (prt_level.ge.1) THEN
-       print*,'WARNING !!! TEST ',comment
-      endif
-      return
-
-!  test sur la hauteur des thermiques ...
-         do i=1,ngrid
-!IMtemp           if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
-           if (prt_level.ge.10) then
-               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
-               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
-               do k=1,nlay
-                  write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
-               enddo
-           endif
-         enddo
-
-
-      return
-      end
-
-! nrlmd le 10/04/2012   Transport de la TKE par le thermique moyen pour la fermeture en ALP 
-!                       On transporte pbl_tke pour donner therm_tke
-!                       Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
-
-!=======================================================================
-!///////////////////////////////////////////////////////////////////////
-!=======================================================================
-
-      subroutine thermcell_tke_transport( &
-     &     ngrid,nlay,ptimestep,fm0,entr0,rg,pplev,  &   ! in
-     &     therm_tke_max)                                ! out
-      USE thermcell_ini_mod, ONLY: prt_level
-      implicit none
-
-!=======================================================================
-!
-!   Calcul du transport verticale dans la couche limite en presence
-!   de "thermiques" explicitement representes
-!   calcul du dq/dt une fois qu'on connait les ascendances
-!
-!=======================================================================
-
-      integer ngrid,nlay
-
-      real, intent(in) :: ptimestep
-      real, intent(in), dimension(ngrid,nlay+1) :: fm0,pplev
-      real, intent(in), dimension(ngrid,nlay) :: entr0
-      real, intent(in) :: rg
-      real, intent(out), dimension(ngrid,nlay) :: therm_tke_max
-
-      real detr0(ngrid,nlay)
-      real masse0(ngrid,nlay)
-      real masse(ngrid,nlay),fm(ngrid,nlay+1)
-      real entr(ngrid,nlay)
-      real q(ngrid,nlay)
-      integer lev_out                           ! niveau pour les print
-
-      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
-      integer ig,k
-
-
-      lev_out=0
-
-
-      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
-
-!   calcul du detrainement
-      do k=1,nlay
-         detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
-         masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
-      enddo
-
-
-! Decalage vertical des entrainements et detrainements.
-      masse(:,1)=0.5*masse0(:,1)
-      entr(:,1)=0.5*entr0(:,1)
-      detr(:,1)=0.5*detr0(:,1)
-      fm(:,1)=0.
-      do k=1,nlay-1
-         masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
-         entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
-         detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
-         fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
-      enddo
-      fm(:,nlay+1)=0.
-
-
-   q(:,:)=therm_tke_max(:,:)
-!!! nrlmd le 16/09/2010
-      do ig=1,ngrid
-         qa(ig,1)=q(ig,1)
-      enddo
-!!!
-
-    if (1==1) then
-      do k=2,nlay
-         do ig=1,ngrid
-            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
-     &         1.e-5*masse(ig,k)) then
-         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
-     &         /(fm(ig,k+1)+detr(ig,k))
-            else
-               qa(ig,k)=q(ig,k)
-            endif
-            if (qa(ig,k).lt.0.) then
-!               print*,'qa<0!!!'
-            endif
-            if (q(ig,k).lt.0.) then
-!               print*,'q<0!!!'
-            endif
-         enddo
-      enddo
-
-! Calcul du flux subsident
-
-      do k=2,nlay
-         do ig=1,ngrid
-            wqd(ig,k)=fm(ig,k)*q(ig,k)
-            if (wqd(ig,k).lt.0.) then
-!               print*,'wqd<0!!!'
-            endif
-         enddo
-      enddo
-      do ig=1,ngrid
-         wqd(ig,1)=0.
-         wqd(ig,nlay+1)=0.
-      enddo
-
-! Calcul des tendances
-      do k=1,nlay
-         do ig=1,ngrid
-            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
-     &               -wqd(ig,k)+wqd(ig,k+1))  &
-     &               *ptimestep/masse(ig,k)
-         enddo
-      enddo
-
- endif
-
-   therm_tke_max(:,:)=q(:,:)
-
-      return
-!!! fin nrlmd le 10/04/2012
-     end
-
