Index: LMDZ4/trunk/libf/phylmd/physiq.F
===================================================================
--- LMDZ4/trunk/libf/phylmd/physiq.F	(revision 864)
+++ LMDZ4/trunk/libf/phylmd/physiq.F	(revision 878)
@@ -177,4 +177,21 @@
       REAL fm_therm(klon,klev+1)
       REAL entr_therm(klon,klev)
+      real,allocatable,save :: clwcon0th(:,:),rnebcon0th(:,:)
+c$OMP THREADPRIVATE(clwcon0th,rnebcon0th)
+
+      real weak_inversion(klon),dthmin(klon)
+      real seuil_inversion
+      save seuil_inversion
+c$OMP THREADPRIVATE(seuil_inversion)
+      integer iflag_ratqs
+      save iflag_ratqs
+c$OMP THREADPRIVATE(iflag_ratqs)
+
+      integer lmax_th(klon)
+      integer limbas(klon)
+      real ratqscth(klon,klev)
+      real ratqsdiff(klon,klev)
+      real zqsatth(klon,klev)
+
 c======================================================================
 c
@@ -929,4 +946,8 @@
       REAL fluxu(klon,klev, nbsrf)   ! flux turbulent de vitesse u
       REAL fluxv(klon,klev, nbsrf)   ! flux turbulent de vitesse v
+
+      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: pbl_tke ! turb kinetic energy
+c   !$OMP THREADPRIVATE(pbl_tke)
+
 c
       REAL zxfluxt(klon, klev)
@@ -1050,4 +1071,5 @@
 c$OMP THREADPRIVATE(d_u_con,d_v_con)
       REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev)
+      REAL d_t_ajsb(klon,klev), d_q_ajsb(klon,klev)
       REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
       REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
@@ -1471,4 +1493,5 @@
 #endif
       allocate( clwcon0(klon,klev),rnebcon0(klon,klev))
+      allocate( clwcon0th(klon,klev),rnebcon0th(klon,klev))
       allocate( tau_ae(klon,klev,2), piz_ae(klon,klev,2))
       allocate( cg_ae(klon,klev,2))
@@ -1507,4 +1530,6 @@
          CALL suphec ! initialiser constantes et parametres phys.
       ENDIF
+
+      print*,'CONVERGENCE PHYSIQUE THERM 1 '
 
 
@@ -1567,8 +1592,7 @@
 c
          call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
-     .                  ok_instan, ok_hf,
+     .                  ok_instan, ok_hf, seuil_inversion,
      .                  fact_cldcon, facttemps,ok_newmicro,
-cIM  .                  iflag_cldcon,ratqsbas,ratqshaut, if_ebil,
-     .                  iflag_cldcon,ratqsbas,ratqshaut,
+     .                  iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,
      .                  ok_ade, ok_aie, 
      .                  bl95_b0, bl95_b1,
@@ -1581,4 +1605,16 @@
          itap    = 0
          itaprad = 0
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! Un petit travail à faire ici.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+         if (iflag_pbl>1) then
+             PRINT*, "Using method MELLOR&YAMADA" 
+         endif
+          ! NB! pbl_tke could/should be read and written from (re)startphy.nc
+          ALLOCATE(pbl_tke(klon,klev+1,nbsrf))  
+          pbl_tke(:,:,:) = 1.e-8 
+
 
          CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
@@ -1589,5 +1625,12 @@
      .       radsol,clesphy0,
      .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
-     .       t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon)
+     .       t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon,
+     .       pbl_tke)
+
+
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
 
        DO i=1,klon
@@ -1620,5 +1663,6 @@
      .                        pdtphys
             abort_message='Pas physique n est pas correct '
-            call abort_gcm(modname,abort_message,1)
+!           call abort_gcm(modname,abort_message,1)
+            dtime=pdtphys
          ENDIF
          IF (nlon .NE. klon) THEN
@@ -1670,9 +1714,20 @@
          ENDIF
 
+           rugoro=0.
 c34EK
          IF (ok_orodr) THEN
-           DO i=1,klon
-             rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
-           ENDDO
+
+           rugoro=0.
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
+! justement quand ok_orodr = false.
+! ce rugoro est utilise par la couche limite et fait double emploi
+! avec les paramétrisations spécifiques de Francois Lott.
+!           DO i=1,klon
+!             rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
+!           ENDDO
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
            CALL SUGWD(klon,klev,paprs,pplay)
            DO i=1,klon
@@ -2088,5 +2143,5 @@
      d     wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
      -     dsens,     devap,     zxsnow,
-     -     zxfluxt,   zxfluxq,   q2m,     fluxq )
+     -     zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
 c
 c
@@ -2260,4 +2315,8 @@
           DO i = 1, klon
             ema_pct(i)  = paprs(i,itop_con(i))
+            if (itop_con(i).gt.klev-3) then
+               print*,'La convection monte trop haut '
+               print*,'itop_con(,',i,',)=',itop_con(i)
+            endif
           ENDDO
           DO i = 1, klon
@@ -2343,8 +2402,16 @@
 c===================================================================
 c
+       call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
+     s ,seuil_inversion,weak_inversion,dthmin)
+
+
+
+      d_t_ajsb(:,:)=0.
+      d_q_ajsb(:,:)=0.
       d_t_ajs(:,:)=0.
       d_u_ajs(:,:)=0.
       d_v_ajs(:,:)=0.
       d_q_ajs(:,:)=0.
+      clwcon0th(:,:)=0.
       fm_therm(:,:)=0.
       entr_therm(:,:)=0.
@@ -2357,22 +2424,68 @@
 c  ====
          IF(prt_level>9)WRITE(lunout,*)'pas de convection'
-      else if(iflag_thermals.eq.0) then
-
-c  Ajustement sec
-c  ==============
-         IF(prt_level>9)WRITE(lunout,*)'ajsec'
-         CALL ajsec(paprs, pplay, t_seri,q_seri, d_t_ajs, d_q_ajs)
-         t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
-         q_seri(:,:) = q_seri(:,:) + d_q_ajs(:,:)
+
+
       else
+
 c  Thermiques
 c  ==========
          IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
      s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
+         print*,'JUSTE AVANT , iflag_thermals='
+     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
+
+
+         if (iflag_thermals.gt.1) then
          call calltherm(pdtphys
-     s      ,pplay,paprs,pphi
-     s      ,u_seri,v_seri,t_seri,q_seri
+     s      ,pplay,paprs,pphi,weak_inversion
+     s      ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
      s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
-     s      ,fm_therm,entr_therm)
+     s      ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth,
+     s       ratqsdiff,zqsatth)
+         endif
+
+!        call calltherm(pdtphys
+!    s      ,pplay,paprs,pphi
+!    s      ,u_seri,v_seri,t_seri,q_seri
+!    s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
+!    s      ,fm_therm,entr_therm)
+
+c  Ajustement sec
+c  ==============
+
+! Dans le cas où on active les thermiques, on fait partir l'ajustement
+! a partir du sommet des thermiques.
+! Dans le cas contraire, on demarre au niveau 1.
+
+         if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
+
+         if(iflag_thermals.eq.0) then
+            IF(prt_level>9)WRITE(lunout,*)'ajsec'
+            limbas(:)=1
+         else
+            limbas(:)=lmax_th(:)
+         endif
+  
+! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
+! pour des test de convergence numerique.
+! Le nouveau ajsec est a priori mieux, meme pour le cas 
+! iflag_thermals = 0 (l'ancienne version peut faire des tendances
+! non nulles numeriquement pour des mailles non concernees.
+
+         if (iflag_thermals.eq.0) then
+            CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
+     s      , d_t_ajsb, d_q_ajsb)
+         else
+            CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
+     s      , d_t_ajsb, d_q_ajsb)
+         endif
+
+         t_seri(:,:) = t_seri(:,:) + d_t_ajsb(:,:)
+         q_seri(:,:) = q_seri(:,:) + d_q_ajsb(:,:)
+         d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
+         d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
+
+         endif
+
       endif
 c
@@ -2406,34 +2519,99 @@
          enddo
          enddo
+      else if (iflag_cldcon.eq.4) then
+         ptconvth(:,:)=.false.
+         ratqsc(:,:)=0.
+         print*,'avant clouds_gno thermique'
+         call clouds_gno
+     s   (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
+         print*,' CLOUDS_GNO OK'
       endif
 
 c   ratqs stables
 c   -------------
-      do k=1,klev
-cIM RAJOUT boucle do=i
-       do i=1, klon
-cIM         ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
-cIM     s   min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.) 
-         ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
-     s   min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.) 
-cIM      print*,' IMratqs STABLE i, k',i,k,ratqss(i,k)
-       enddo 
-      enddo
+
+      if (iflag_ratqs.eq.0) then
+
+! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
+         do k=1,klev
+            do i=1, klon
+               ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
+     s         min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.) 
+            enddo 
+         enddo
+
+! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de 
+! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
+! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
+! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
+! Il s'agit de differents tests dans la phase de reglage du modele
+! avec thermiques.
+
+      else if (iflag_ratqs.eq.1) then
+
+         do k=1,klev
+            do i=1, klon
+               if (pplay(i,k).ge.60000.) then
+                  ratqss(i,k)=ratqsbas
+               else if ((pplay(i,k).ge.30000.).and.
+     s            (pplay(i,k).lt.60000.)) then
+                  ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
+     s            (60000.-pplay(i,k))/(60000.-30000.)
+               else
+                  ratqss(i,k)=ratqshaut
+               endif
+            enddo
+         enddo
+
+      else
+
+         do k=1,klev
+            do i=1, klon
+               if (pplay(i,k).ge.60000.) then
+                  ratqss(i,k)=ratqsbas
+     s            *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
+               else if ((pplay(i,k).ge.30000.).and.
+     s             (pplay(i,k).lt.60000.)) then
+                    ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
+     s              (60000.-pplay(i,k))/(60000.-30000.)
+               else
+                    ratqss(i,k)=ratqshaut
+               endif
+            enddo
+         enddo
+      endif
+
+
 
 
 c  ratqs final
 c  -----------
-      if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
-c   les ratqs sont une conbinaison de ratqss et ratqsc
-c   ratqs final
-c   1e4 (en gros 3 heures), en dur pour le moment, est le temps de
-c   relaxation des ratqs
-c         facttemps=exp(-pdtphys/1.e4)
-         facteur=exp(-pdtphys*facttemps)
-         ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
-         ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
-c         print*,'calcul des ratqs fini'
+
+      if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
+     s    .or.iflag_cldcon.eq.4) then
+
+! On ajoute une constante au ratqsc*2 pour tenir compte de 
+! fluctuations turbulentes de petite echelle
+
+         do k=1,klev
+            do i=1,klon
+               if ((fm_therm(i,k).gt.1.e-10)) then
+                  ratqsc(i,k)=sqrt(ratqsc(i,k)**2+0.05**2)
+               endif
+            enddo
+         enddo
+
+!   les ratqs sont une conbinaison de ratqss et ratqsc
+!   1800s, en dur pour le moment, est le temps de
+!   relaxation des ratqs
+
+         facteur=exp(-pdtphys/1800.)
+
+         print*,'WARNING ratqs a revoir '
+         ratqs(:,:)=ratqsc(:,:)*(1.-facteur)+ratqs(:,:)*facteur
+         ratqs(:,:)=sqrt(ratqs(:,:)**2+ratqss(:,:)**2)
+
       else
-c   on ne prend que le ratqs stable pour fisrtilp
+!   on ne prend que le ratqs stable pour fisrtilp
          ratqs(:,:)=ratqss(:,:)
       endif
@@ -2505,5 +2683,5 @@
 cIM cf FH
 c     IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
-       IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
+      IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
        snow_tiedtke=0.
 c     print*,'avant calcul de la pseudo precip '
@@ -2542,9 +2720,8 @@
       ENDDO
 
-      ELSE IF (iflag_cldcon.eq.3) THEN
+      ELSE IF (iflag_cldcon.ge.3) THEN
 c  On prend pour les nuages convectifs le max du calcul de la
 c  convection et du calcul du pas de temps precedent diminue d'un facteur
 c  facttemps
-c      facttemps=pdtphys/1.e4
       facteur = pdtphys *facttemps
       do k=1,klev
@@ -3314,5 +3491,6 @@
      .      radsol,
      .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
-     .      t_ancien, q_ancien, rnebcon, ratqs, clwcon)
+     .      t_ancien, q_ancien, rnebcon, ratqs, clwcon,
+     .      pbl_tke)
       ENDIF
       
