Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/calfis_p.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/calfis_p.F	(revision 1247)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/calfis_p.F	(revision 1250)
@@ -158,6 +158,4 @@
       REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
 c
-c      REAL,ALLOCATABLE,SAVE :: pvervel(:,:)
-c
       REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
       REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
@@ -175,5 +173,4 @@
       REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
       REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
-c      REAL,ALLOCATABLE,SAVE :: pvervel_omp(:,:) 
       REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
       REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
@@ -252,5 +249,4 @@
       ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
       ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
-c      ALLOCATE(pvervel(klon,llm))
       ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm))
       ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
@@ -278,4 +274,6 @@
 
 c$OMP MASTER             
+!CDIR ON_ADB(index_i)
+!CDIR ON_ADB(index_j) 
       do ig0=1,klon
         i=index_i(ig0)
@@ -300,4 +298,6 @@
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l = 1, llmp1
+!CDIR ON_ADB(index_i)
+!CDIR ON_ADB(index_j) 
         do ig0=1,klon
           i=index_i(ig0)
@@ -314,5 +314,6 @@
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
-
+!CDIR ON_ADB(index_i)
+!CDIR ON_ADB(index_j) 
         do ig0=1,klon
           i=index_i(ig0)
@@ -321,5 +322,4 @@
           zplay(ig0,l)   = preff * pksurcp ** unskap
           ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
-c          pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
         enddo
 
@@ -335,4 +335,6 @@
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
          DO l=1,llm
+!CDIR ON_ADB(index_i)
+!CDIR ON_ADB(index_j) 
            do ig0=1,klon
              i=index_i(ig0)
@@ -344,20 +346,4 @@
       ENDDO
 
-c   convergence dynamique pour les traceurs "EAU"
-! Earth-specific treatment of first 2 tracers (water)
-      if (planet_type=="earth") then
-       DO iq=1,2
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-         DO l=1,llm
-           do ig0=1,klon
-             i=index_i(ig0)
-             j=index_j(ig0)
-c             pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
-           enddo
-         ENDDO
-c$OMP END DO NOWAIT	 
-       ENDDO
-      endif ! of if (planet_type=="earth")
-
 
 c   Geopotentiel calcule par rapport a la surface locale:
@@ -378,19 +364,4 @@
 c$OMP END DO NOWAIT
       
-c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
-c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux 
-c      de masse est calclue dans advtrac_p.F  
-c
-cc$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
-c      DO l=1,llm
-c        do ig0=1,klon
-c           i=index_i(ig0)
-c           j=index_j(ig0)
-c           pvervel(ig0,l) = pw(i,j,l)*g* unsaire(i,j)
-c        enddo
-c	if (is_north_pole) pvervel(1,l)=pw(1,1,l)*g /apoln
-c	if (is_south_pole) pvervel(klon,l)=pw(1,jjp1,l)*g/apols
-c      ENDDO
-cc$OMP END DO NOWAIT
 
 c
@@ -406,4 +377,7 @@
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
       DO l=1,llm
+!CDIR ON_ADB(index_i)
+!CDIR ON_ADB(index_j) 
+!CDIR SPARSE
         do ig0=kstart,kend
           i=index_i(ig0)
@@ -412,19 +386,19 @@
             zufi(ig0,l)= 0.5 *(  pucov(iim,j,l)/cu(iim,j)
      $                         + pucov(1,j,l)/cu(1,j) )
-c            pcvgu(ig0,l)= 0.5*(  pducov(iim,j,l)/cu(iim,j) 
-c     $                         + pducov(1,j,l)/cu(1,j) )
           else
             zufi(ig0,l)= 0.5*(  pucov(i-1,j,l)/cu(i-1,j) 
      $                       + pucov(i,j,l)/cu(i,j) )
-c            pcvgu(ig0,l)= 0.5*(  pducov(i-1,j,l)/cu(i-1,j) 
-c     $                        + pducov(i,j,l)/cu(i,j) )
           endif
         enddo
       ENDDO
 c$OMP END DO NOWAIT
+
 c   46.champ v:
 c   -----------
+
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
+!CDIR ON_ADB(index_i)
+!CDIR ON_ADB(index_j) 
         DO ig0=kstart,kend
           i=index_i(ig0)
@@ -433,6 +407,4 @@
      $                       + pvcov(i,j,l)/cv(i,j) )
     
-c          pcvgv(ig0+i,l)= 0.5 * (  pdvcov(i,j-1,l)/cv(i,j-1) 
-c     $                           + pdvcov(i,j,l)/cv(i,j) )
          ENDDO
       ENDDO
@@ -449,21 +421,15 @@
 
            z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
-c           z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
            DO i=2,iim
               z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
-c              z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
            ENDDO
   
            DO i=1,iim
               zcos(i)   = COS(rlonv(i))*z1(i)
-c              zcosbis(i)= COS(rlonv(i))*z1bis(i)
               zsin(i)   = SIN(rlonv(i))*z1(i)
-c              zsinbis(i)= SIN(rlonv(i))*z1bis(i)
            ENDDO
   
            zufi(1,l)  = SSUM(iim,zcos,1)/pi
-c           pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
            zvfi(1,l)  = SSUM(iim,zsin,1)/pi
-c           pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
   
         ENDDO
@@ -482,22 +448,15 @@
   
          z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
-c         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
            DO i=2,iim
-           z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
-c           z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
+             z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
 	   ENDDO
   
            DO i=1,iim
               zcos(i)    = COS(rlonv(i))*z1(i)
-c              zcosbis(i) = COS(rlonv(i))*z1bis(i)
               zsin(i)    = SIN(rlonv(i))*z1(i)
-c              zsinbis(i) = SIN(rlonv(i))*z1bis(i)
 	   ENDDO
   
            zufi(klon,l)  = SSUM(iim,zcos,1)/pi
-c           pcvgu(klon,l) = SSUM(iim,zcosbis,1)/pi
            zvfi(klon,l)  = SSUM(iim,zsin,1)/pi
-c           pcvgv(klon,l) = SSUM(iim,zsinbis,1)/pi
-
         ENDDO
 c$OMP END DO NOWAIT       
@@ -506,6 +465,4 @@
 
       IF (is_sequential) THEN
-        if (planet_type=="earth") then
-#ifdef CPP_EARTH
 c
 cIM calcul PV a teta=350, 380, 405K
@@ -514,6 +471,4 @@
      $           ntetaSTD,rtetaSTD,PVteta)
 c
-#endif
-        endif 
       ENDIF
 
@@ -525,15 +480,4 @@
 c   ---------------------
 
-cc$OMP  PARALLEL DEFAULT(NONE) 
-cc$OMP+ PRIVATE(i,l,offset,iq)
-cc$OMP+ SHARED(klon_omp_nb,nqtot,klon_omp_begin,
-cc$OMP+        debut,lafin,rdayvrai,heure,dtphys,zplev,zplay,
-cc$OMP+        zphi,zphis,presnivs,clesphy0,zufi,zvfi,ztfi,
-cc$OMP+        zqfi,pvervel,zdufi,zdvfi,zdtfi,zdqfi,zdpsrf)
-
-c PRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
-c c$OMP+                 presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
-c c$OMP+                 zqfi_omp,pvervel_omp,zdufi_omp,zdvfi_omp,
-c c$OMP+                 zdtfi_omp,zdqfi_omp,zdpsrf_omp)
 
 c$OMP BARRIER
@@ -550,5 +494,4 @@
         allocate(ztfi_omp(klon,llm))
         allocate(zqfi_omp(klon,llm,nqtot))
-c        allocate(pvervel_omp(klon,llm))
         allocate(zdufi_omp(klon,llm))
         allocate(zdvfi_omp(klon,llm))
@@ -617,10 +560,4 @@
       enddo
 	
-c      do l=1,llm
-c        do i=1,klon
-c	  pvervel_omp(i,l)=pvervel(offset+i,l)
-c	enddo 
-c      enddo
-	
       do l=1,llm
         do i=1,klon
@@ -660,5 +597,4 @@
       
 c$OMP BARRIER
-cym      call WriteField_phy_p('zdtfi_omp',zdtfi_omp(:,:),llm)
       
       if (planet_type=="earth") then
@@ -681,5 +617,4 @@
      .             ztfi_omp,
      .             zqfi_omp,
-c     .             pvervel_omp,
 c#ifdef INCA
      .             flxwfi_omp,
@@ -695,7 +630,4 @@
 #endif
       endif !of if (planet_type=="earth")
-
-cym      call WriteField_phy_p('zdtfi_omp',zdtfi_omp(:,:),llm)
-
 c$OMP BARRIER
 
@@ -753,10 +685,4 @@
 	enddo 
       enddo
-	
-c      do l=1,llm
-c        do i=1,klon
-c	  pvervel(offset+i,l)=pvervel_omp(i,l)
-c	enddo 
-c      enddo
 	
       do l=1,llm
@@ -791,5 +717,4 @@
       
 
-cc$OMP END PARALLEL
       klon=klon_mpi
 500   CONTINUE
@@ -797,5 +722,4 @@
 
 c$OMP MASTER
-cym      call WriteField_phy('zdtfi',zdtfi(:,:),llm)
       call stop_timer(timer_physic)
 c$OMP END MASTER
@@ -913,5 +837,7 @@
       DO l=1,llm
 
-!!cdir NODEP
+!CDIR ON_ADB(index_i)
+!CDIR ON_ADB(index_j) 
+!cdir NODEP
         do ig0=kstart,kend
           i=index_i(ig0)
@@ -975,11 +901,12 @@
 
 C
-
+!cdir NODEP
       DO iq=1,nqtot
          iiq=niadv(iq)
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
          DO l=1,llm
-
-!!cdir NODEP           
+!CDIR ON_ADB(index_i)
+!CDIR ON_ADB(index_j) 
+!cdir NODEP           
 	     DO ig0=kstart,kend
               i=index_i(ig0)
@@ -1009,5 +936,7 @@
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
-!!cdir NODEP
+!CDIR ON_ADB(index_i)
+!CDIR ON_ADB(index_j) 
+!cdir NODEP
          do ig0=kstart,kend
            i=index_i(ig0)
@@ -1052,5 +981,7 @@
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
       DO l=1,llm
-!!cdir NODEP
+!CDIR ON_ADB(index_i)
+!CDIR ON_ADB(index_j) 
+!cdir NODEP
         do ig0=kstart,kend
            i=index_i(ig0)
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/cray.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/cray.F	(revision 1247)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/cray.F	(revision 1250)
@@ -13,11 +13,17 @@
       real sx((n-1)*incx+1),sy((n-1)*incy+1)
 c
+      if (incx.eq.1.and.incy.eq.1) then
+      do 10 i=1,n
+         sy(i)=sx(i)
+10    continue
+      else
       iy=1
       ix=1
-      do 10 i=1,n
+      do 11 i=1,n
          sy(iy)=sx(ix)
          ix=ix+incx
          iy=iy+incy
-10    continue
+11    continue
+      endif
 c
       return
@@ -32,9 +38,15 @@
 c
       ssum=0.
+      if (incx.eq.1) then
+      do 10 i=1,n
+         ssum=ssum+sx(i)
+10    continue
+      else
       ix=1
-      do 10 i=1,n
+      do 11 i=1,n
          ssum=ssum+sx(ix)
          ix=ix+incx
-10    continue
+11    continue
+      endif
 c
       return
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/groupeun_p.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/groupeun_p.F	(revision 1247)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/groupeun_p.F	(revision 1250)
@@ -1,4 +1,5 @@
       SUBROUTINE groupeun_p(jjmax,llmax,jjb,jje,q)
       USE parallel
+      USE Write_Field_p
       IMPLICIT NONE
 
@@ -17,5 +18,5 @@
       REAL airecs,qs
 
-      INTEGER i,j,l,ig,j1,j2,i0,jd
+      INTEGER i,j,l,ig,ig2,j1,j2,i0,jd
 
 c--------------------------------------------------------------------c 
@@ -37,4 +38,7 @@
       LOGICAL, SAVE :: first = .TRUE.
 !$OMP THREADPRIVATE(first)
+      INTEGER,SAVE :: i_index(iim,ngroup)
+      INTEGER      :: offset
+      REAL         :: qsum(iim/ngroup)
 
       IF (first) THEN
@@ -54,32 +58,67 @@
             j_start  = MAX(jjb, j1-jd)
             j_finish = MIN(jje, j2-jd)
-            DO j=j_start, j_finish
-               DO i0=1,iim,2**(ngroup-ig+1)
-                  qn=0.
-                  DO i=i0,i0+2**(ngroup-ig+1)-1
-                     qn=qn+q(i,j,l)
-                  ENDDO
-                  DO i=i0,i0+2**(ngroup-ig+1)-1
-                     q(i,j,l)=qn*airen_tab(i,j,jd)
-                  ENDDO
+            DO ig2=1,ngroup-ig+1
+              offset=2**(ig2-1)
+              DO j=j_start, j_finish
+!CDIR NODEP
+!CDIR ON_ADB(q)
+                 DO i0=1,iim,2**ig2
+                   q(i0,j,l)=q(i0,j,l)+q(i0+offset,j,l) 
+                 ENDDO
+              ENDDO
+            ENDDO
+            
+            DO j=j_start, j_finish
+!CDIR NODEP
+!CDIR ON_ADB(q)
+               DO i=1,iim
+                 q(i,j,l)=q(i-MOD(i-1,2**(ngroup-ig+1)),j,l)
+               ENDDO
+            ENDDO
+
+            DO j=j_start, j_finish
+!CDIR ON_ADB(airen_tab)
+!CDIR ON_ADB(q)
+               DO i=1,iim
+                 q(i,j,l)=q(i,j,l)*airen_tab(i,j,jd)
                ENDDO
                q(iip1,j,l)=q(1,j,l)
             ENDDO
-        
+       
 !c     Concerne le pole sud
             j_start  = MAX(1+jjp1-jje-jd, j1-jd)
             j_finish = MIN(1+jjp1-jjb-jd, j2-jd)
-            DO j=j_start, j_finish
-               DO i0=1,iim,2**(ngroup-ig+1)
-                  qs=0.
-                  DO i=i0,i0+2**(ngroup-ig+1)-1
-                     qs=qs+q(i,jjp1-j+1-jd,l)
-                  ENDDO
-                  DO i=i0,i0+2**(ngroup-ig+1)-1
-                     q(i,jjp1-j+1-jd,l)=qs*aires_tab(i,jjp1-j+1,jd)
-                  ENDDO
+            DO ig2=1,ngroup-ig+1
+              offset=2**(ig2-1)
+              DO j=j_start, j_finish
+!CDIR NODEP
+!CDIR ON_ADB(q)
+                 DO i0=1,iim,2**ig2
+                   q(i0,jjp1-j+1-jd,l)= q(i0,jjp1-j+1-jd,l)
+     &                                 +q(i0+offset,jjp1-j+1-jd,l) 
+                 ENDDO
+              ENDDO
+            ENDDO
+
+
+            DO j=j_start, j_finish
+!CDIR NODEP
+!CDIR ON_ADB(q)
+               DO i=1,iim
+                 q(i,jjp1-j+1-jd,l)=q(i-MOD(i-1,2**(ngroup-ig+1)),
+     &                                jjp1-j+1-jd,l)
+               ENDDO
+            ENDDO
+
+            DO j=j_start, j_finish
+!CDIR ON_ADB(aires_tab)
+!CDIR ON_ADB(q)
+               DO i=1,iim
+                 q(i,jjp1-j+1-jd,l)=q(i,jjp1-j+1-jd,l)*  
+     &                              aires_tab(i,jjp1-j+1,jd)
                ENDDO
                q(iip1,jjp1-j+1-jd,l)=q(1,jjp1-j+1-jd,l)
             ENDDO
+
         
             j1=j2+1
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/leapfrog_p.F	(revision 1247)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/leapfrog_p.F	(revision 1250)
@@ -327,5 +327,5 @@
 c$OMP MASTER
       ItCount=ItCount+1
-      if (MOD(ItCount,1)==1) then
+      if (MOD(ItCount,12)==0) then
         debug=.true.
       else
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/pres2lev.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/pres2lev.F	(revision 1247)
+++ 	(revision )
@@ -1,84 +1,0 @@
-! $Id$
-!
-c******************************************************
-      SUBROUTINE   pres2lev(varo,varn,lmo,lmn,po,pn,
-     %                      ni,nj,ok_invertp)
-c
-c interpolation lineaire pour passer
-c a une nouvelle discretisation verticale pour
-c les variables de GCM
-c Francois Forget (01/1995)
-c MOdif remy roca 12/97 pour passer de pres2sig
-c Modif F.Codron 07/08 po en 3D
-c**********************************************************
-
-      IMPLICIT NONE
-
-c   Declarations:
-c ==============
-c
-c  ARGUMENTS
-c  """""""""
-       LOGICAL, INTENT(IN) :: ok_invertp
-       INTEGER, INTENT(IN) :: lmo ! dimensions ancienne couches
-       INTEGER, INTENT(IN) :: lmn ! dimensions nouvelle couches
-       INTEGER lmomx ! dimensions ancienne couches
-       INTEGER lmnmx ! dimensions nouvelle couches
-
-       parameter(lmomx=10000,lmnmx=10000)
-
-        real, INTENT(IN) :: po(ni,nj,lmo) ! niveau de pression ancienne grille
-        real, INTENT(IN) :: pn(ni,nj,lmn) ! niveau de pression nouvelle grille
-
-       INTEGER, INTENT(IN) :: ni,nj ! nombre de point horizontale
-
-       REAL, INTENT(IN)  :: varo(ni,nj,lmo) ! var dans l'ancienne grille
-       REAL, INTENT(OUT) :: varn(ni,nj,lmn) ! var dans la nouvelle grille
-
-       real zvaro(lmomx),zpo(lmomx)
-
-c Autres variables
-c """"""""""""""""
-       INTEGER n, ln ,lo, i, j, Nhoriz
-       REAL coef
-
-c run
-c ====
-        do i=1,ni
-        do j=1,nj
-! Inversion de l'ordre des niveaux verticaux
-          IF (ok_invertp) THEN
-           do lo=1,lmo
-              zpo(lo)=po(i,j,lmo+1-lo)
-              zvaro(lo)=varo(i,j,lmo+1-lo)
-           enddo
-          ELSE
-           do lo=1,lmo
-              zpo(lo)=po(i,j,lo)
-              zvaro(lo)=varo(i,j,lo)
-           enddo
-          ENDIF 
-
-           do ln=1,lmn
-              if (pn(i,j,ln).ge.zpo(1))then
-                 varn(i,j,ln) =  zvaro(1)
-              else if (pn(i,j,ln).le.zpo(lmo)) then
-                 varn(i,j,ln) =  zvaro(lmo)
-              else
-                 do lo=1,lmo-1 
-                    if ( (pn(i,j,ln).le.zpo(lo)).and.
-     &                 (pn(i,j,ln).gt.zpo(lo+1)) )then
-                       coef=(pn(i,j,ln)-zpo(lo))
-     &                 /(zpo(lo+1)-zpo(lo))
-                       varn(i,j,ln)=zvaro(lo)
-     &                 +coef*(zvaro(lo+1)-zvaro(lo))
-c       print*,'pn(',ln,')=',pn(i,j,ln),varn(i,j,ln)
-                    end if
-                 enddo           
-              endif
-           enddo
-
-        enddo
-        enddo
-      return
-      end    
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/pres2lev.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/pres2lev.F90	(revision 1250)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/pres2lev.F90	(revision 1250)
@@ -0,0 +1,74 @@
+! $Id: pres2lev.F 1179 2009-06-11 14:18:47Z jghattas $
+!
+!******************************************************
+SUBROUTINE pres2lev(varo,varn,lmo,lmn,po,pn,ni,nj,ok_invertp)
+!
+! interpolation lineaire pour passer
+! a une nouvelle discretisation verticale pour
+! les variables de GCM
+! Francois Forget (01/1995)
+! MOdif remy roca 12/97 pour passer de pres2sig
+! Modif F.Codron 07/08 po en 3D
+!**********************************************************
+
+  IMPLICIT NONE
+
+!   Declarations:
+! ==============
+!
+!  ARGUMENTS
+!  """""""""
+  LOGICAL, INTENT(IN) :: ok_invertp
+  INTEGER, INTENT(IN) :: lmo ! dimensions ancienne couches
+  INTEGER, INTENT(IN) :: lmn ! dimensions nouvelle couches
+  
+  REAL, INTENT(IN) :: po(ni*nj,lmo) ! niveau de pression ancienne grille
+  REAL, INTENT(IN) :: pn(ni*nj,lmn) ! niveau de pression nouvelle grille
+
+  INTEGER, INTENT(IN) :: ni,nj ! nombre de point horizontal
+
+  REAL, INTENT(IN)  :: varo(ni*nj,lmo) ! var dans l'ancienne grille
+  REAL, INTENT(OUT) :: varn(ni*nj,lmn) ! var dans la nouvelle grille
+
+  REAL :: zvaro(ni*nj,lmo),zpo(ni*nj,lmn)
+
+! Autres variables
+! """"""""""""""""
+  INTEGER ::  ln ,lo, k
+  REAL    :: coef
+
+
+! Inversion de l'ordre des niveaux verticaux
+  IF (ok_invertp) THEN
+    DO lo=1,lmo
+      DO k=1,ni*nj
+        zpo(k,lo)=po(k,lmo+1-lo)
+        zvaro(k,lo)=varo(k,lmo+1-lo)
+      ENDDO
+    ENDDO
+  ELSE
+    DO lo=1,lmo
+      DO k=1,ni*nj
+        zpo(k,lo)=po(k,lo)
+        zvaro(k,lo)=varo(k,lo)
+      ENDDO
+    ENDDO
+  ENDIF 
+
+  DO ln=1,lmn
+    DO lo=1,lmo-1
+      DO k=1,ni*nj
+        IF (pn(k,ln) >= zpo(k,1) ) THEN
+          varn(k,ln) = varo(k,1)
+        ELSE IF (pn(k,ln) <= zpo(k,lmo)) THEN
+          varn(k,ln) = zvaro(k,lmo)
+        ELSE IF ( pn(k,ln) <= zpo(k,lo) .AND. pn(k,ln) > zpo(k,lo+1) ) THEN
+          coef = (pn(k,ln)-zpo(k,lo)) / (zpo(k,lo+1)-zpo(k,lo))
+          varn(k,ln) = zvaro(k,lo) + coef*(zvaro(k,lo+1)-zvaro(k,lo))
+        ENDIF
+         
+      ENDDO  
+    ENDDO
+  ENDDO                
+
+END SUBROUTINE pres2lev    
