Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/compile
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/compile	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/compile	(revision 226)
@@ -29,8 +29,8 @@
 nz=25
 ##################
-#tracers=2
-#nx=64
-#ny=48
-#nz=25
+tracers=2
+nx=64
+ny=48
+nz=25
 ##################
 
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/physiq.F
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/physiq.F	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/physiq.F	(revision 226)
@@ -1,10 +1,11 @@
-      SUBROUTINE physiq(ngrid,nlayer,nq,
-     $            firstcall,lastcall,
-     $            pday,ptime,ptimestep,
-     $            pplev,pplay,pphi,
-     $            pu,pv,pt,pq,
-     $            pw,
-     $            pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
-
+      SUBROUTINE physiq(
+     $            ngrid,nlayer,nq
+     $            ,firstcall,lastcall
+     $            ,pday,ptime,ptimestep
+     $            ,pplev,pplay,pphi
+     $            ,pu,pv,pt,pq
+     $            ,pw
+     $            ,pdu,pdv,pdt,pdq,pdpsrf,tracerdyn
+     $            )
 
       IMPLICIT NONE
@@ -187,5 +188,6 @@
                                      !   of the size distribution
 c     Albedo of deposited surface ice
-      REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45
+      !!REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45
+      REAL, PARAMETER :: alb_surfice = 0.45 !!TESTS_JB 
 
       SAVE day_ini, icount
@@ -278,4 +280,5 @@
       REAL ccn(ngridmx,nlayermx)   ! Cloud condensation nuclei
                                    !   (particules kg-1)
+      SAVE ccn  !! in case iradia != 1
       real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m)
       real qtot1,qtot2 ! total aerosol mass
@@ -298,4 +301,15 @@
 
       REAL time_phys
+
+c Variables from thermal
+
+      REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
+      REAL wmax_th(ngridmx)
+      REAL ,SAVE :: hfmax_th(ngridmx)
+      REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
+      REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
+      INTEGER lmax_th(ngridmx)
+      REAL dtke_th(ngridmx,nlayermx+1)
+      REAL dummycol(ngridmx)
 
 c=======================================================================
@@ -513,4 +527,5 @@
 c          Radiative transfer
 c          ------------------
+
            CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
      $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
@@ -549,6 +564,4 @@
            ENDIF
 
-
-
         ENDIF ! of if(mod(icount-1,iradia).eq.0)
 
@@ -565,5 +578,4 @@
                fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
            ENDDO
-
 
          DO l=1,nlayer
@@ -597,7 +609,6 @@
 c    4. Vertical diffusion (turbulent mixing):
 c    -----------------------------------------
-c
+
       IF (calldifv) THEN
-
 
          DO ig=1,ngrid
@@ -612,5 +623,5 @@
             enddo
          enddo
-         
+
 c        Calling vdif (Martian version WITH CO2 condensation)
          CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
@@ -633,7 +644,7 @@
          ENDDO
 
-         DO ig=1,ngrid
-            zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
-         ENDDO
+          DO ig=1,ngrid
+             zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
+          ENDDO
 
          if (tracer) then 
@@ -659,4 +670,44 @@
       ENDIF ! of IF (calldifv)
 
+c-----------------------------------------------------------------------
+c   TEST. Thermals :
+c HIGHLY EXPERIMENTAL, BEWARE !!
+c   -----------------------------
+ 
+      if(calltherm) then
+ 
+        call calltherm_interface(firstcall,
+     $ long,lati,zzlev,zzlay,
+     $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
+     $ pplay,pplev,pphi,zpopsk,
+     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
+     $ dtke_th,hfmax_th,wmax_th)
+ 
+         DO l=1,nlayer
+           DO ig=1,ngrid
+              pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
+              pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
+              pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
+              q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
+           ENDDO
+        ENDDO
+ 
+        DO ig=1,ngrid
+          q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
+        ENDDO      
+    
+        if (tracer) then
+        DO iq=1,nq
+         DO l=1,nlayer
+           DO ig=1,ngrid
+             pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
+           ENDDO
+         ENDDO
+        ENDDO
+        endif
+
+        else   !of if calltherm
+        lmax_th(:)=0
+        end if
 
 c-----------------------------------------------------------------------
@@ -676,9 +727,10 @@
 
          CALL convadj(ngrid,nlayer,nq,ptimestep,
-     $                pplay,pplev,zpopsk,
+     $                pplay,pplev,zpopsk,lmax_th,
      $                pu,pv,zh,pq,
      $                pdu,pdv,zdh,pdq,
      $                zduadj,zdvadj,zdhadj,
      $                zdqadj)
+
 
          DO l=1,nlayer
@@ -713,5 +765,5 @@
      $              co2ice,albedo,emis,
      $              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
-     $	            fluxsurf_sw,zls) 
+     $              fluxsurf_sw,zls) 
 
          DO l=1,nlayer
@@ -1021,7 +1073,7 @@
       ENDDO
 
-c    Compute surface stress : (NB: z0 is a common in planete.h)
+c    Compute surface stress : (NB: z0 is a common in surfdat.h)
 c     DO ig=1,ngrid
-c        cd = (0.4/log(zzlay(ig,1)/z0))**2
+c        cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
 c        zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2)
 c     ENDDO
@@ -1029,6 +1081,6 @@
 c     Sum of fluxes in solar spectral bands (for output only)
       DO ig=1,ngrid
-	     fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
-	     fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
+             fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
+             fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
       ENDDO
 c ******* TEST ******************************************************
@@ -1070,7 +1122,7 @@
 
       IF (ngrid.NE.1) THEN
-         print*,'Ls =',zls*180./pi,
-     &   ' tauref(700 Pa,lat=0) =',tauref(ngrid/2),
-     &   ' tau(Viking1) =',tau(ig_vl1,1)
+         print*,'Ls =',zls*180./pi
+     &   ,' tauref(700 Pa,lat=0) =',tauref(ngrid/2)
+     &   ,' tau(Viking1) =',tau(ig_vl1,1)
 
 
@@ -1142,5 +1194,4 @@
 c        which can later be used to make the statistic files of the run:
 c        "stats")          only possible in 3D runs !
-
          
          IF (callstats) THEN
@@ -1257,6 +1308,8 @@
 c        WRITEDIAGFI can ALSO be called from any other subroutines
 c        for any variables !!
-        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
+         call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
      &                  emis)
+!         call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
+!         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
          call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
      &                  tsurf)
@@ -1264,4 +1317,6 @@
         call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
      &                  co2ice)
+c         call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
+c     &                  "K",2,zt(1,7))
 c         call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
 c     &                  fluxsurf_lw)
@@ -1273,5 +1328,4 @@
 c     &                  fluxtop_sw_tot)
          call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
-c        call WRITEDIAGFI(ngrid,"tau","tau"," ",2,tau)
         call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
         call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
@@ -1279,5 +1333,5 @@
 c         call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
 c        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
-c        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
+!        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
 c        call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
 c        call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
@@ -1287,4 +1341,7 @@
 c        call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
 c    &                   'w.m-2',3,zdtlw)
+c        CALL WRITEDIAGFI(ngridmx,'tauTESap',
+c     &                         'tau abs 825 cm-1',
+c     &                         '',2,tauTES)
 
 !!!!!!!!!!!!!!!!!!!!!!!!SOIL SOIL SOIL
@@ -1324,40 +1381,40 @@
 
             !!!! waterice = q01, voir readmeteo.F90
-            call WRITEDIAGFI(ngridmx,'q01',noms(iq),
+            call WRITEDIAGFI(ngridmx,'q01',noms(igcm_h2o_ice),
      &                      'kg/kg',3,
      &                       zq(1:ngridmx,1:nlayermx,igcm_h2o_ice))
             !!!! watervapor = q02, voir readmeteo.F90
-            call WRITEDIAGFI(ngridmx,'q02',noms(iq),
+            call WRITEDIAGFI(ngridmx,'q02',noms(igcm_h2o_vap),
      &                      'kg/kg',3,
      &                       zq(1:ngridmx,1:nlayermx,igcm_h2o_vap))
-
-c            call WRITEDIAGFI(ngridmx,'qsurf'//str2,noms(iq),
-c     &                     'kg.m-2',2,qsurf(1,iq))
-
-
-c             CALL WRITEDIAGFI(ngridmx,'mtot',
-c     &                       'total mass of water vapor',
-c     &                       'kg/m2',2,mtot)
-c             CALL WRITEDIAGFI(ngridmx,'icetot',
-c     &                       'total mass of water ice',
-c     &                       'kg/m2',2,icetot)
-cc            vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
-cc    &                *mugaz/mmol(igcm_h2o_ice)
-cc            call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
-cc    &                       'mol/mol',3,vmr)
-c             CALL WRITEDIAGFI(ngridmx,'reffice',
-c     &                       'Mean reff',
-c     &                       'm',2,rave)
-cc            call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
-cc    &                       'm',3,rice)
-cc            If activice is true, tauTES is computed in aeropacity.F;
-c             if (.not.activice) then
-c               CALL WRITEDIAGFI(ngridmx,'tauTESap',
-c     &                         'tau abs 825 cm-1',
-c     &                         '',2,tauTES)
-c             endif
-c             call WRITEDIAGFI(ngridmx,'h2o_ice_s',
-c     &                       'surface h2o_ice',
-c     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
+            !!!! surface waterice qsurf02 (voir readmeteo)
+            call WRITEDIAGFI(ngridmx,'qsurf02','surface tracer',
+     &                      'kg.m-2',2,
+     &                       qsurf(1:ngridmx,igcm_h2o_ice))
+
+             CALL WRITEDIAGFI(ngridmx,'mtot',
+     &                       'total mass of water vapor',
+     &                       'kg/m2',2,mtot)
+             CALL WRITEDIAGFI(ngridmx,'icetot',
+     &                       'total mass of water ice',
+     &                       'kg/m2',2,icetot)
+c            vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
+c    &                *mugaz/mmol(igcm_h2o_ice)
+c            call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
+c    &                       'mol/mol',3,vmr)
+             CALL WRITEDIAGFI(ngridmx,'reffice',
+     &                       'Mean reff',
+     &                       'm',2,rave)
+c            call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
+c    &                       'm',3,rice)
+c            If activice is true, tauTES is computed in aeropacity.F;
+             if (.not.activice) then
+               CALL WRITEDIAGFI(ngridmx,'tauTESap',
+     &                         'tau abs 825 cm-1',
+     &                         '',2,tauTES)
+             endif
+             call WRITEDIAGFI(ngridmx,'h2o_ice_s',
+     &                       'surface h2o_ice',
+     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
            endif !(water)
 
@@ -1401,5 +1458,5 @@
              call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
      &                        'kg/kg',3,pq(1,1,igcm_dust_mass))
-            call WRITEDIAGFI(ngridmx,'dustN','Dust number',
+             call WRITEDIAGFI(ngridmx,'dustN','Dust number',
      &                        'part/kg',3,pq(1,1,igcm_dust_number))
            else
@@ -1419,4 +1476,40 @@
 
 c        ----------------------------------------------------------
+c        Outputs of thermals
+c        ----------------------------------------------------------
+         if (calltherm) then
+
+!        call WRITEDIAGFI(ngrid,'dtke',
+!     &              'tendance tke thermiques','m**2/s**2',
+!     &                         3,dtke_th)
+!        call WRITEDIAGFI(ngrid,'d_u_ajs',
+!     &              'tendance u thermiques','m/s',
+!     &                         3,pdu_th*ptimestep)
+!        call WRITEDIAGFI(ngrid,'d_v_ajs',
+!     &              'tendance v thermiques','m/s',
+!     &                         3,pdv_th*ptimestep)
+!        if (tracer) then
+!        if (nq .eq. 2) then
+!        call WRITEDIAGFI(ngrid,'deltaq_th',
+!     &              'delta q thermiques','kg/kg',
+!     &                         3,ptimestep*pdq_th(:,:,2))
+!        endif
+!        endif
+
+        lmax_th_out(:)=real(lmax_th(:))
+
+        call WRITEDIAGFI(ngridmx,'lmax_th',
+     &              'hauteur du thermique','K',
+     &                         2,lmax_th_out)
+        call WRITEDIAGFI(ngridmx,'hfmax_th',
+     &              'maximum TH heat flux','K.m/s',
+     &                         2,hfmax_th)
+        call WRITEDIAGFI(ngridmx,'wmax_th',
+     &              'maximum TH vertical velocity','m/s',
+     &                         2,wmax_th)
+
+         endif
+
+c        ----------------------------------------------------------
 c        Output in netcdf file "diagsoil.nc" for subterranean
 c          variables (output every "ecritphy", as for writediagfi)
@@ -1452,4 +1545,51 @@
 c        CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
 
+! THERMALS STUFF 1D
+
+         if(calltherm) then
+
+        lmax_th_out(:)=real(lmax_th(:))
+
+        call WRITEDIAGFI(ngridmx,'lmax_th',
+     &              'hauteur du thermique','point',
+     &                         0,lmax_th_out)
+        call WRITEDIAGFI(ngridmx,'hfmax_th',
+     &              'maximum TH heat flux','K.m/s',
+     &                         0,hfmax_th)
+        call WRITEDIAGFI(ngridmx,'wmax_th',
+     &              'maximum TH vertical velocity','m/s',
+     &                         0,wmax_th)
+
+
+         co2col(:)=0.
+         dummycol(:)=0.
+         if (tracer) then
+         do l=1,nlayermx
+           do ig=1,ngrid
+             co2col(ig)=co2col(ig)+
+     &                  zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
+         if (nqmx .gt. 1) then
+             iq=2 ! to avoid out-of-bounds spotting by picky compilers
+                  ! when gcm is compiled with only one tracer
+             dummycol(ig)=dummycol(ig)+
+     &                  zq(ig,l,iq)*(pplev(ig,l)-pplev(ig,l+1))/g
+         endif
+         enddo
+         enddo
+
+         end if
+         call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass'          &
+     &                                      ,'kg/m-2',0,co2col)
+         call WRITEDIAGFI(ngrid,'dummycol','integrated dummy mass'      &
+     &                                      ,'kg/m-2',0,dummycol)
+         endif
+         call WRITEDIAGFI(ngrid,'w','vertical velocity'                 &
+     &                              ,'m/s',1,pw)
+         call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
+         call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
+     &                  tsurf)
+
+         call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
+         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
 ! or output in diagfi.nc (for testphys1d)
          call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
