Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 3156)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 3157)
@@ -343,4 +343,5 @@
       REAL zzlay(ngrid,nlayer)     ! altitude at the middle of the layers
       REAL zzlev(ngrid,nlayer+1)   ! altitude at layer boundaries
+      REAL gz(ngrid,nlayer)        ! variation of g with altitude from aeroid surface
 !      REAL latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
 
@@ -393,4 +394,6 @@
       INTEGER length
       PARAMETER (length=100)
+      REAL tlaymean ! temporary value of mean layer temperature for zzlay
+
 
 c     Variables for the total dust for diagnostics
@@ -836,4 +839,6 @@
       dwatercap(:,:)=0
 
+
+
       call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf,
      &       albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg)
@@ -869,25 +874,48 @@
       ps(:) = pplev(:,1)
 
+
+#ifndef MESOSCALE
+c-----------------------------------------------------------------------
+c    1.2.1 Compute mean mass, cp, and R
+c    concentrations outputs rnew(ngrid,nlayer), cpnew(ngrid,nlayer)
+c	, mmean(ngrid,nlayer) and Akknew(ngrid,nlayer)
+c    --------------------------------
+
+      if(photochem.or.callthermos) then
+         call concentrations(ngrid,nlayer,nq,
+     &                       zplay,pt,pdt,pq,pdq,ptimestep)
+      endif
+#endif
+
 c     Compute geopotential at interlayers
 c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 c     ponderation des altitudes au niveau des couches en dp/p
-
-      DO l=1,nlayer
-         DO ig=1,ngrid
-            zzlay(ig,l)=pphi(ig,l)/g
-         ENDDO
+cc ------------------------------------------
+	!Calculation zzlev & zzlay for first layer      
+      DO ig=1,ngrid
+      	 zzlay(ig,1)=-(log(pplay(ig,1)/ps(ig)))*rnew(ig,1)*pt(ig,1)/g                                      
+	 zzlev(ig,1)=0
+         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
+         gz(ig,1)=g
+      	
+        DO l=2,nlayer     
+         ! compute "mean" temperature of the layer
+          if(pt(ig,l) .eq. pt(ig,l-1)) then
+             tlaymean=pt(ig,l)
+          else
+             tlaymean=(pt(ig,l)- pt(ig,l-1))/log(pt(ig,l)/pt(ig,l-1))
+          endif
+                     
+          ! compute gravitational acceleration (at altitude zaeroid(nlayer-1))
+          gz(ig,l)=g*(rad**2)/(rad+zzlay(ig,l-1)+(phisfi(ig)/g))**2
+          
+          zzlay(ig,l)=zzlay(ig,l-1)- 
+     &   (log(pplay(ig,l)/pplay(ig,l-1))*rnew(ig,l)*tlaymean/gz(ig,l))
+                 
+          z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
+          z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
+          zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
+        ENDDO
       ENDDO
-      DO ig=1,ngrid
-         zzlev(ig,1)=0.
-         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
-      ENDDO
-      DO l=2,nlayer
-         DO ig=1,ngrid
-            z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
-            z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
-            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
-         ENDDO
-      ENDDO
-
 
 !     Potential temperature calculation not the same in physiq and dynamic
@@ -902,14 +930,4 @@
       ENDDO
 
-#ifndef MESOSCALE
-c-----------------------------------------------------------------------
-c    1.2.5 Compute mean mass, cp, and R
-c    --------------------------------
-
-      if(photochem.or.callthermos) then
-         call concentrations(ngrid,nlayer,nq,
-     &                       zplay,pt,pdt,pq,pdq,ptimestep)
-      endif
-#endif
 
       ! Compute vertical velocity (m/s) from vertical mass flux
@@ -1349,5 +1367,5 @@
      &                      pdqrds,wspeed,dsodust,dsords,dsotop,
      &                      tau_pref_scenario,tau_pref_gcm)
-
+     
 c      update the tendencies of both dust after vertical transport
          DO l=1,nlayer
@@ -1609,5 +1627,5 @@
      $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
      $ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
-
+       
          DO l=1,nlayer
            DO ig=1,ngrid
@@ -2061,4 +2079,5 @@
      &                tau,tauscaling)
 
+
 c Flux at the surface of co2 ice computed in co2cloud microtimestep
            IF (rdstorm) THEN
@@ -2282,12 +2301,37 @@
         ENDDO
         zplev(:,nlayer+1) = 0.
-        ! update layers altitude
-        DO l=2,nlayer
-          DO ig=1,ngrid
-           z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
-           z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
-           zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
+        
+ 	! Calculation of zzlay and zzlay with udpated pressure and temperature
+	DO ig=1,ngrid
+      	 zzlay(ig,1)=-(log(zplay(ig,1)/ps(ig)))*rnew(ig,1)*
+     &		(pt(ig,1)+pdt(ig,1)*ptimestep) /g                                      
+
+      	  DO l=2,nlayer
+      	         
+           ! compute "mean" temperature of the layer
+            if((pt(ig,l)+pdt(ig,l)*ptimestep) .eq.
+     &               (pt(ig,l-1)+pdt(ig,l-1)*ptimestep)) then
+               tlaymean= pt(ig,l)+pdt(ig,l)*ptimestep
+            else
+               tlaymean=((pt(ig,l)+pdt(ig,l)*ptimestep)- 
+     &		  (pt(ig,l-1)+pdt(ig,l-1)*ptimestep))/
+     &		  log((pt(ig,l)+pdt(ig,l)*ptimestep)/ 
+     &		  (pt(ig,l-1)+pdt(ig,l-1)*ptimestep))
+            endif
+                     
+            ! compute gravitational acceleration (at altitude zaeroid(nlayer-1))
+            gz(ig,l)=g*(rad**2)/(rad+zzlay(ig,l-1)+(phisfi(ig)/g))**2
+          
+
+            zzlay(ig,l)=zzlay(ig,l-1)- 
+     &     (log(zplay(ig,l)/zplay(ig,l-1))*rnew(ig,l)*tlaymean/gz(ig,l))
+              
+        
+        !   update layers altitude
+             z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
+             z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
+             zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
           ENDDO
-        ENDDO
+        ENDDO 
 #endif
       ENDIF  ! of IF (callcond)