Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 882)
+++ trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 883)
@@ -128,4 +128,5 @@
 #include "tracer.h"
 #include "nlteparams.h"
+#include "comvert.h"
 
 #include "chimiedata.h"
@@ -306,5 +307,7 @@
 
       real co2col(ngridmx)        ! CO2 column
-      REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
+      ! pplev and pplay are dynamical inputs and must not be modified in the physics.
+      ! instead, use zplay and zplev :
+      REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx) 
       REAL zstress(ngrid), cd
       real hco2(nqmx),tmean, zlocal(nlayermx)
@@ -511,4 +514,11 @@
       end if
 
+c     Initialize pressure levels
+c     ~~~~~~~~~~~~~~~~~~
+      zplev(:,:) = pplev(:,:)
+      zplay(:,:) = pplay(:,:)
+      ps(:) = pplev(:,1)
+
+
 c     Compute geopotential at interlayers
 c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -526,6 +536,6 @@
       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))
+            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
@@ -539,5 +549,5 @@
       DO l=1,nlayer
          DO ig=1,ngrid 
-            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
+            zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp
             zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
          ENDDO
@@ -550,5 +560,5 @@
 
       if(photochem.or.callthermos) then
-         call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
+         call concentrations(zplay,pt,pdt,pq,pdq,ptimestep)
       endif
 #endif
@@ -581,5 +591,5 @@
            IF(callnlte) then
               if(nltemodel.eq.0.or.nltemodel.eq.1) then
-                 CALL nltecool(ngrid,nlayer,nq,pplay,pt,pq,zdtnlte)
+                 CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte)
               else if(nltemodel.eq.2) then
                 co2vmr_gcm(1:ngrid,1:nlayer)=
@@ -596,5 +606,5 @@
      &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
                  
-                 CALL nlte_tcool(ngrid,nlayer,pplay*9.869e-6,
+                 CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
      $                pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm, 
      $                ovmr_gcm,  zdtnlte )
@@ -609,5 +619,5 @@
 c          Find number of layers for LTE radiation calculations
            IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
-     &          CALL nlthermeq(ngrid,nlayer,pplev,pplay)
+     &          CALL nlthermeq(ngrid,nlayer,zplev,zplay)
 
 c          Note: Dustopacity.F has been transferred to callradite.F
@@ -622,5 +632,5 @@
 
            CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
-     $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
+     $     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
      $     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
      $     tauref,tau,aerosol,tauscaling,taucloudtes,rdust,rice,
@@ -635,5 +645,5 @@
      &             " tau(",f4.0," Pa) =",f9.6)')
      &            odpref,tauref(igout),
-     &            odpref,tau(igout,1)*odpref/pplev(igout,1)
+     &            odpref,tau(igout,1)*odpref/zplev(igout,1)
 c          ---------------------------------------------------------
 c          Call slope parameterization for direct and scattered flux
@@ -680,5 +690,5 @@
            zdtnirco2(:,:)=0
            if (callnirco2) then
-              call nirco2abs (ngrid,nlayer,pplay,dist_sol,nq,pq,
+              call nirco2abs (ngrid,nlayer,zplay,dist_sol,nq,pq,
      .                       mu0,fract,declin, zdtnirco2)
            endif
@@ -696,5 +706,5 @@
 c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            IF(callnlte) THEN
-              CALL blendrad(ngrid, nlayer, pplay,
+              CALL blendrad(ngrid, nlayer, zplay,
      &             zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
            ELSE
@@ -742,5 +752,5 @@
 
         CALL calldrag_noro(ngrid,nlayer,ptimestep,
-     &                 pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
+     &                 zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
  
         DO l=1,nlayer
@@ -799,5 +809,5 @@
          IF (tke_heat_flux .ne. 0.) THEN
              zz1(:)=(pt(:,1)+pdt(:,1)*ptimestep)*(r/g)*
-     &                 (-alog(pplay(:,1)/pplev(:,1)))
+     &                 (-alog(zplay(:,1)/zplev(:,1)))
              pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
          ENDIF
@@ -806,5 +816,5 @@
          CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
      $        ptimestep,capcal,lwrite,
-     $        pplay,pplev,zzlay,zzlev,z0,
+     $        zplay,zplev,zzlay,zzlev,z0,
      $        pu,pv,zh,pq,tsurf,emis,qsurf,
      $        zdum1,zdum2,zdh,pdq,zflubid,
@@ -879,5 +889,5 @@
      $ zzlev,zzlay,
      $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
-     $ pplay,pplev,pphi,zpopsk,
+     $ zplay,zplev,pphi,zpopsk,
      $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
      $ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
@@ -931,5 +941,5 @@
 
          CALL convadj(ngrid,nlayer,nq,ptimestep,
-     $                pplay,pplev,zpopsk,lmax_th,
+     $                zplay,zplev,zpopsk,lmax_th,
      $                pu,pv,zh,pq,
      $                pdu,pdv,zdh,pdq,
@@ -969,8 +979,11 @@
          co2ice = co2ice * 10000.
       ENDIF
+      
+      
+      pdpsrf(:) = 0
 
       IF (callcond) THEN
          CALL newcondens(ngrid,nlayer,nq,ptimestep,
-     $              capcal,pplay,pplev,tsurf,pt,
+     $              capcal,zplay,zplev,tsurf,pt,
      $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
      $              co2ice,albedo,emis,
@@ -999,5 +1012,30 @@
          ENDIF ! of IF (tracer)
 
+        ! update surface pressure
+        DO ig=1,ngrid
+          ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
+        ENDDO
+        
+        ! update pressure levels
+        DO l=1,nlayer
+         DO ig=1,ngrid
+          zplay(ig,l) = aps(l) + bps(l)*ps(ig)
+          zplev(ig,l) = ap(l)  + bp(l)*ps(ig)
+         ENDDO
+        ENDDO
+        zplev(:,l) = 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)
+          ENDDO
+        ENDDO
+      
       ENDIF  ! of IF (callcond)
+      
+
 
 c-----------------------------------------------------------------------
@@ -1016,5 +1054,5 @@
 
            call watercloud(ngrid,nlayer,ptimestep,
-     &                pplev,pplay,pdpsrf,zzlay, pt,pdt,
+     &                zplev,zplay,pdpsrf,zzlay, pt,pdt,
      &                pq,pdq,zdqcloud,zdtcloud,
      &                nq,tau,tauscaling,rdust,rice,nuice,
@@ -1037,5 +1075,30 @@
 
 ! increment dust and ccn masses and numbers
+! We need to check that we have Nccn & Ndust > 0
+! This is due to single precision rounding problems
            if (microphys) then
+             pdq(1:ngrid,1:nlayer,igcm_ccn_mass) = 
+     &         pdq(1:ngrid,1:nlayer,igcm_ccn_mass) + 
+     &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
+             pdq(1:ngrid,1:nlayer,igcm_ccn_number) = 
+     &         pdq(1:ngrid,1:nlayer,igcm_ccn_number) + 
+     &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
+             where (pq(:,:,igcm_ccn_mass) + 
+     &       ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
+               pdq(:,:,igcm_ccn_mass) = 
+     &         - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
+               pdq(:,:,igcm_ccn_number) = 
+     &         - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
+             end where
+             where (pq(:,:,igcm_ccn_number) + 
+     &       ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
+               pdq(:,:,igcm_ccn_mass) = 
+     &         - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
+               pdq(:,:,igcm_ccn_number) = 
+     &         - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
+             end where
+           endif
+
+           if (scavenging) then
              pdq(1:ngrid,1:nlayer,igcm_dust_mass) = 
      &         pdq(1:ngrid,1:nlayer,igcm_dust_mass) + 
@@ -1044,46 +1107,19 @@
      &         pdq(1:ngrid,1:nlayer,igcm_dust_number) + 
      &         zdqcloud(1:ngrid,1:nlayer,igcm_dust_number)
-             pdq(1:ngrid,1:nlayer,igcm_ccn_mass) = 
-     &         pdq(1:ngrid,1:nlayer,igcm_ccn_mass) + 
-     &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
-             pdq(1:ngrid,1:nlayer,igcm_ccn_number) = 
-     &         pdq(1:ngrid,1:nlayer,igcm_ccn_number) + 
-     &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
-           endif
-
-! We need to check that we have Nccn & Ndust > 0
-! This is due to single precision rounding problems
-         if (scavenging) then
-
-           where (pq(:,:,igcm_ccn_mass) + 
-     &      ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
-              pdq(:,:,igcm_ccn_mass) = 
-     &        - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
-              pdq(:,:,igcm_ccn_number) = 
-     &        - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
-           end where
-           where (pq(:,:,igcm_ccn_number) + 
-     &      ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
-              pdq(:,:,igcm_ccn_mass) = 
-     &        - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
-              pdq(:,:,igcm_ccn_number) = 
-     &        - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
-           end where
-           where (pq(:,:,igcm_dust_mass) + 
-     &      ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
-              pdq(:,:,igcm_dust_mass) = 
-     &        - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
-              pdq(:,:,igcm_dust_number) = 
-     &        - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
-           end where
-           where (pq(:,:,igcm_dust_number) + 
-     &      ptimestep*pdq(:,:,igcm_dust_number) < 0.)
-              pdq(:,:,igcm_dust_mass) = 
-     &        - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
-              pdq(:,:,igcm_dust_number) = 
-     &        - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
-           end where
-           
-         endif ! of if scavenging
+             where (pq(:,:,igcm_dust_mass) + 
+     &       ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
+               pdq(:,:,igcm_dust_mass) = 
+     &         - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
+               pdq(:,:,igcm_dust_number) = 
+     &         - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
+             end where
+             where (pq(:,:,igcm_dust_number) + 
+     &       ptimestep*pdq(:,:,igcm_dust_number) < 0.)
+               pdq(:,:,igcm_dust_mass) = 
+     &         - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
+               pdq(:,:,igcm_dust_number) = 
+     &         - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
+             end where
+           endif ! of if scavenging
                       
 
@@ -1097,5 +1133,5 @@
 c        ----------
          IF(callddevil) then 
-           call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
+           call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
      &                zdqdev,zdqsdev)
  
@@ -1127,5 +1163,5 @@
 
            call callsedim(ngrid,nlayer, ptimestep,
-     &                pplev,zzlev, zzlay, pt, rdust, rice,
+     &                zplev,zzlev, zzlay, pt, rdust, rice,
      &                rsedcloud,rhocloud,
      &                pq, pdq, zdqsed, zdqssed,nq, 
@@ -1162,10 +1198,10 @@
 
 !           dust and ice surface area
-            call surfacearea(ngrid, nlayer, ptimestep, pplay, zzlay, 
+            call surfacearea(ngrid, nlayer, ptimestep, zplay, zzlay, 
      $                       pt, pq, pdq, nq, 
      $                       rdust, rice, tau, tauscaling, 
      $                       surfdust, surfice)
 !           call photochemistry
-            call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
+            call calchim(ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
      $                   zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
      $                   zdqcloud,zdqscloud,tauref,co2ice,
@@ -1233,5 +1269,5 @@
 
       if (callthermos) then
-        call thermosphere(pplev,pplay,dist_sol,
+        call thermosphere(zplev,zplay,dist_sol,
      $     mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
      &     pt,pq,pu,pv,pdt,pdq,
@@ -1275,7 +1311,7 @@
          if (caps.and.(obliquit.lt.27.)) then
            ! NB: Updated surface pressure, at grid point 'ngrid', is
-           !     ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
+           !     ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
            tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095*
-     &                     (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
+     &                     (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
          endif
 #endif
@@ -1342,17 +1378,4 @@
             zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
           ENDDO
-        ENDDO
-      ENDDO
-
-      ! surface pressure
-      DO ig=1,ngrid
-          ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
-      ENDDO
-
-      ! pressure
-      DO l=1,nlayer
-        DO ig=1,ngrid
-             zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
-             zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
         ENDDO
       ENDDO
@@ -1413,5 +1436,5 @@
          WRITE(*,'(2f10.5,2f15.5)')
      s   tsurf(igout),zdtsurf(igout)*ptimestep,
-     s   pplev(igout,1),pdpsrf(igout)*ptimestep
+     s   zplev(igout,1),pdpsrf(igout)*ptimestep
          WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
          WRITE(*,'(i4,6f10.5)') (l,
@@ -1545,8 +1568,8 @@
                  mtot(ig) = mtot(ig) + 
      &                      zq(ig,l,igcm_h2o_vap) * 
-     &                      (pplev(ig,l) - pplev(ig,l+1)) / g
+     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
                  icetot(ig) = icetot(ig) + 
      &                        zq(ig,l,igcm_h2o_ice) * 
-     &                        (pplev(ig,l) - pplev(ig,l+1)) / g
+     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
 c                Computing abs optical depth at 825 cm-1 in each
 c                  layer to simulate NEW TES retrieval
@@ -1556,5 +1579,5 @@
                  opTES(ig,l)= 0.75 * Qabsice * 
      &             zq(ig,l,igcm_h2o_ice) *
-     &             (pplev(ig,l) - pplev(ig,l+1)) / g
+     &             (zplev(ig,l) - zplev(ig,l+1)) / g
      &             / (rho_ice * rice(ig,l) * (1.+nuice_ref))
                  tauTES(ig)=tauTES(ig)+ opTES(ig,l) 
@@ -1563,5 +1586,5 @@
 c               if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
              enddo
-             call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
+             call watersat(ngridmx*nlayermx,zt,zplay,zqsat)
              satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
 
@@ -1574,8 +1597,8 @@
                     Nccntot(ig) = Nccntot(ig) + 
      &              zq(ig,l,igcm_ccn_number)*tauscaling(ig)
-     &              *(pplev(ig,l) - pplev(ig,l+1)) / g
+     &              *(zplev(ig,l) - zplev(ig,l+1)) / g
                     Mccntot(ig) = Mccntot(ig) + 
      &              zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
-     &              *(pplev(ig,l) - pplev(ig,l+1)) / g
+     &              *(zplev(ig,l) - zplev(ig,l+1)) / g
 cccc Column integrated effective ice radius 
 cccc is weighted by total ice surface area (BETTER than total ice mass) 
@@ -1583,5 +1606,5 @@
      &                      tauscaling(ig) *
      &                      zq(ig,l,igcm_ccn_number) *
-     &                      (pplev(ig,l) - pplev(ig,l+1)) / g * 
+     &                      (zplev(ig,l) - zplev(ig,l+1)) / g * 
      &                      rice(ig,l) * rice(ig,l)*  (1.+nuice_ref)
                  enddo
@@ -1596,5 +1619,5 @@
                  rave(ig) = rave(ig) + 
      &                      zq(ig,l,igcm_h2o_ice) *
-     &                      (pplev(ig,l) - pplev(ig,l+1)) / g * 
+     &                      (zplev(ig,l) - zplev(ig,l+1)) / g * 
      &                      rice(ig,l) * (1.+nuice_ref)
                  enddo
@@ -1641,5 +1664,5 @@
      &                "m.s-1",3,pw)
         call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
-        call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
+        call wstats(ngrid,"pressure","Pressure","Pa",3,zplay)
           call wstats(ngrid,"q2",
      &                "Boundary layer eddy kinetic energy",
@@ -1795,5 +1818,5 @@
                       do ig=1,ngrid                                  
                          colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
-     $                                  *(pplev(ig,l)-pplev(ig,l+1)) 
+     $                                  *(zplev(ig,l)-zplev(ig,l+1)) 
      $                                  *6.022e22/(mmol(iq)*g)       
                       end do                                        
@@ -1848,5 +1871,5 @@
         dustot(ig) = dustot(ig) +
      &               zq(ig,l,igcm_dust_mass)
-     &               *  (pplev(ig,l) - pplev(ig,l+1)) / g
+     &               *  (zplev(ig,l) - zplev(ig,l+1)) / g
        enddo
       enddo
@@ -1945,5 +1968,5 @@
            do ig=1,ngrid
              co2col(ig)=co2col(ig)+
-     &                  zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
+     &                  zq(ig,l,igcm_co2)*(zplev(ig,l)-zplev(ig,l+1))/g
            enddo
          enddo
@@ -2217,5 +2240,5 @@
            do ig=1,ngrid
              co2col(ig)=co2col(ig)+
-     &                  zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
+     &                  zq(ig,l,1)*(zplev(ig,l)-zplev(ig,l+1))/g
          enddo
          enddo
@@ -2285,5 +2308,5 @@
                opTES(1,l)= 0.75 * Qabsice *
      &             zq(1,l,igcm_h2o_ice) *
-     &             (pplev(1,l) - pplev(1,l+1)) / g
+     &             (zplev(1,l) - zplev(1,l+1)) / g
      &             / (rho_ice * rice(1,l) * (1.+nuice_ref))
                tauTES=tauTES+ opTES(1,l)
@@ -2305,7 +2328,7 @@
            do l=1,nlayer
              mtot = mtot +  zq(1,l,igcm_h2o_vap) 
-     &                 * (pplev(1,l) - pplev(1,l+1)) / g
+     &                 * (zplev(1,l) - zplev(1,l+1)) / g
              icetot = icetot +  zq(1,l,igcm_h2o_ice) 
-     &                 * (pplev(1,l) - pplev(1,l+1)) / g
+     &                 * (zplev(1,l) - zplev(1,l+1)) / g
            end do
            h2otot = h2otot+mtot+icetot
@@ -2329,5 +2352,5 @@
              rave = rave + tauscaling(1) *
      &              zq(1,l,igcm_ccn_number) *
-     &              (pplev(1,l) - pplev(1,l+1)) / g * 
+     &              (zplev(1,l) - zplev(1,l+1)) / g * 
      &              rice(1,l) * rice(1,l)*  (1.+nuice_ref)
              enddo
@@ -2335,13 +2358,13 @@
 
               Nccntot= 0
-              call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
+              call watersat(ngridmx*nlayermx,zt,zplay,zqsat)
                do l=1,nlayermx
                    Nccntot = Nccntot + 
      &              zq(1,l,igcm_ccn_number)*tauscaling(1)
-     &              *(pplev(1,l) - pplev(1,l+1)) / g
+     &              *(zplev(1,l) - zplev(1,l+1)) / g
                    satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l)
                    satu(1,l) = (max(satu(1,l),float(1))-1)
 !     &                      * zq(1,l,igcm_h2o_vap) *
-!     &                      (pplev(1,l) - pplev(1,l+1)) / g
+!     &                      (zplev(1,l) - zplev(1,l+1)) / g
                enddo
                call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1,
@@ -2364,5 +2387,5 @@
              do l=1,nlayer
                rave = rave + zq(1,l,igcm_h2o_ice)
-     &              * (pplev(1,l) - pplev(1,l+1)) / g
+     &              * (zplev(1,l) - zplev(1,l+1)) / g
      &              *  rice(1,l) * (1.+nuice_ref)
              enddo
@@ -2395,5 +2418,5 @@
 
 
-         zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
+         zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g
 
          do l=2,nlayer-1
@@ -2402,8 +2425,8 @@
      &        tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
               zlocal(l)= zlocal(l-1)
-     &        -log(pplay(1,l)/pplay(1,l-1))*rnew(1,l)*tmean/g
+     &        -log(zplay(1,l)/zplay(1,l-1))*rnew(1,l)*tmean/g
          enddo
          zlocal(nlayer)= zlocal(nlayer-1)-
-     &                   log(pplay(1,nlayer)/pplay(1,nlayer-1))*
+     &                   log(zplay(1,nlayer)/zplay(1,nlayer-1))*
      &                   rnew(1,nlayer)*tmean/g
 
