Index: trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90	(revision 4081)
+++ trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90	(revision 4082)
@@ -407,4 +407,7 @@
       real zdpsrfmr(ngrid)        ! Pressure tendency for mass_redistribution routine (Pa/s).
 
+      real gz(ngrid,nlayer)       ! Variation of g with altitude from aeroid surface
+      real tlaymean               ! Temporary value of mean layer temperature for zzlay
+
       ! Local variables for MICROPHYSICS:
       ! ---------------------------------
@@ -776,26 +779,46 @@
       endif
 
-      glat(:) = g !AF24: removed oblateness
-
       ! Compute geopotential between layers.
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
-      do l=1,nlayer
-         zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
-      enddo
-
-      zzlev(1:ngrid,1)=0.
-
-      do l=2,nlayer
-         do ig=1,ngrid
+      glat(:) = g !AF24: removed oblateness
+
+      !zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
+      !do l=1,nlayer
+      !   zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
+      !enddo
+      !zzlev(1:ngrid,1)=0.
+      !do l=2,nlayer
+      !   do ig=1,ngrid
+      !      z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
+      !      z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
+      !      zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
+      !   enddo
+      !enddo
+      !Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
+      !zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
+      
+      ! Calculation zzlev & zzlay with g variable (from Mars PCM)
+      do ig = 1, ngrid
+         ! First layer
+         zzlay(ig,1) = -(log(pplay(ig,1)/pplev(ig,1)))*r*pt(ig,1)/g
+         zzlev(ig,1) = 0
+         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))*&
+                        r*tlaymean/gz(ig,l))
             z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
             z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
             zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
          enddo
+         zzlev(ig,nlayer+1) = 2*zzlev(ig,nlayer)-zzlev(ig,nlayer-1)
       enddo
-
-      !Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
-
-      zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
 
       ! Compute potential temperature
@@ -838,7 +861,8 @@
       endif
 
-      !  Compute variations of g with latitude (to do).
+      !  Compute variations of g with altitude (to do).
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       gzlat(:,:) = g
+      !gzlat(:,:) = gz(:,:)
 
       ! Initialize microphysical diagnostics.
@@ -1143,9 +1167,9 @@
               print*, 'Read temp file from ',file_path
               ! levs have been put in km
-              DO ig=1,ngrid
-                CALL interp_line(levdat,tempdat,Nfine,zzlay(ig,:)/1000.,tmean(ig,:),nlayer)
-              enddo
             endif
             if (fixed_temp_prof) then
+               do ig=1,ngrid
+                call interp_line(levdat,tempdat,Nfine,zzlay(ig,:)/1000.,tmean(ig,:),nlayer)
+               enddo
                DO ig=1,ngrid
                  dtrad(ig,1:nlayer)=(tmean(ig,1:nlayer)-(pt(ig,1:nlayer)+    &