Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/dyn3dmem/gcm.F
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/dyn3dmem/gcm.F	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/dyn3dmem/gcm.F	(revision 3825)
@@ -13,7 +13,4 @@
       USE parallel_lmdz
       USE infotrac
-#ifdef CPP_PHYS
-      USE mod_interface_dyn_phys, ONLY: init_interface_dyn_phys
-#endif
       USE mod_hallo
       USE Bands
@@ -28,10 +25,4 @@
 #endif
 
-#ifdef CPP_PHYS
-!      USE mod_grid_phy_lmdz
-!      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
-!      USE dimphy
-!      USE comgeomphy
-#endif
       IMPLICIT NONE
 
@@ -160,12 +151,12 @@
       call Read_Distrib
 
-#ifdef CPP_PHYS
-        CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys,
-     &                      COMM_LMDZ)
+!#ifdef CPP_PHYS
+!        CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys,
+!     &                      COMM_LMDZ)
 !#endif
 !      CALL set_bands
 !#ifdef CPP_PHYS
-      CALL Init_interface_dyn_phys
-#endif
+!      CALL Init_interface_dyn_phys
+!#endif
       CALL barrier
 
@@ -429,5 +420,5 @@
 ! Physics:
 #ifdef CPP_PHYS
-         CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys,
+         CALL iniphysiq(iim,jjm, distrib_phys(mpi_rank), comm_lmdz, llm, daysec, day_ini, dtphys/nsplit_phys,
      &                rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,
      &                iflag_phys)
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/iniphysiq.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/iniphysiq.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/iniphysiq.F90	(revision 3825)
@@ -3,22 +3,17 @@
 
 
-SUBROUTINE iniphysiq(ii,jj,nlayer,punjours, pdayref,ptimestep,         &
+SUBROUTINE iniphysiq(ii, jj, nbp, communicator, nlayer,punjours, pdayref,ptimestep,         &
                      rlatu,rlatv,rlonu,rlonv,aire,cu,cv,               &
                      prad,pg,pr,pcpp,iflag_phys)
+  
   USE dimphy, ONLY: klev ! number of atmospheric levels
-  USE mod_grid_phy_lmdz, ONLY: klon_glo ! number of atmospheric columns
-                                        ! (on full grid)
+  USE mod_grid_phy_lmdz, ONLY: klon_glo,  & ! number of atmospheric columns (on full grid)
+                               regular_lonlat  ! regular longitude-latitude grid type
   USE mod_phys_lmdz_para, ONLY: klon_omp, & ! number of columns (on local omp grid)
                                 klon_omp_begin, & ! start index of local omp subgrid
                                 klon_omp_end, & ! end index of local omp subgrid
                                 klon_mpi_begin ! start indes of columns (on local mpi grid)
-  USE comgeomphy, ONLY: initcomgeomphy, &
-                        initcomgeomphy_vert, &
-                        initcomgeomphy_horiz,&
-                        airephy, & ! physics grid area (m2)
-                        cuphy, & ! cu coeff. (u_covariant = cu * u)
-                        cvphy, & ! cv coeff. (v_covariant = cv * v)
-                        rlond, & ! longitudes
-                        rlatd ! latitudes
+  USE geometry_mod, ONLY : init_geometry
+  USE vertical_layers_mod, ONLY : init_vertical_layers
   USE misc_mod, ONLY: debug
   USE infotrac, ONLY: nqtot,nqo,nbtr,tname,ttext,type_trac,&
@@ -29,4 +24,7 @@
   USE infotrac_phy, ONLY: init_infotrac_phy
   USE phyaqua_mod, ONLY: iniaqua
+  USE physics_distribution_mod, ONLY : init_physics_distribution
+  USE regular_lonlat_mod, ONLY : init_regular_lonlat, east, west, north, south, north_east, north_west, south_west, south_east
+  USE mod_interface_dyn_phys, ONLY :  init_interface_dyn_phys
   IMPLICIT NONE
 
@@ -51,4 +49,6 @@
   INTEGER, INTENT (IN) :: ii ! number of atmospheric columns along longitudes
   INTEGER, INTENT (IN) :: jj ! number of atompsheric columns along latitudes
+  INTEGER, INTENT (IN) :: nbp ! number of physics points (local)
+  INTEGER, INTENT (IN) :: communicator ! mpi communicator
   REAL, INTENT (IN) :: rlatu(jj+1) ! latitudes of the physics grid
   REAL, INTENT (IN) :: rlatv(jj) ! latitude boundaries of the physics grid
@@ -63,11 +63,20 @@
 
   INTEGER :: ibegin, iend, offset
-  INTEGER :: i,j
+  INTEGER :: i,j,k
   CHARACTER (LEN=20) :: modname = 'iniphysiq'
   CHARACTER (LEN=80) :: abort_message
   REAL :: total_area_phy, total_area_dyn
 
+  REAL,ALLOCATABLE :: boundslon_reg(:,:)
+  REAL,ALLOCATABLE :: boundslat_reg(:,:)
 
   ! global array, on full physics grid:
+  REAL,ALLOCATABLE :: latfi_glo(:)
+  REAL,ALLOCATABLE :: lonfi_glo(:)
+  REAL,ALLOCATABLE :: cufi_glo(:)
+  REAL,ALLOCATABLE :: cvfi_glo(:)
+  REAL,ALLOCATABLE :: airefi_glo(:)
+  REAL,ALLOCATABLE :: boundslonfi_glo(:,:)
+  REAL,ALLOCATABLE :: boundslatfi_glo(:,:)
   REAL,ALLOCATABLE :: latfi(:)
   REAL,ALLOCATABLE :: lonfi(:)
@@ -75,49 +84,100 @@
   REAL,ALLOCATABLE :: cvfi(:)
   REAL,ALLOCATABLE :: airefi(:)
-
-  IF (nlayer/=klev) THEN
-    WRITE (lunout, *) 'STOP in ', trim(modname)
-    WRITE (lunout, *) 'Problem with dimensions :'
-    WRITE (lunout, *) 'nlayer     = ', nlayer
-    WRITE (lunout, *) 'klev   = ', klev
-    abort_message = ''
-    CALL abort_gcm(modname, abort_message, 1)
-  END IF
-
-  !call init_phys_lmdz(ii,jj+1,llm,1,(/(jj-1)*ii+2/))
+  REAL,ALLOCATABLE :: boundslonfi(:,:)
+  REAL,ALLOCATABLE :: boundslatfi(:,:)
+
+  CALL init_physics_distribution(regular_lonlat, 4, nbp, ii, jj+1, nlayer, communicator)
+  CALL init_interface_dyn_phys
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! init regular longitude-latitude grid
+
+  ALLOCATE(boundslon_reg(ii,2))
+  ALLOCATE(boundslat_reg(jj+1,2))
+  
+  DO i=1,ii
+   boundslon_reg(i,east)=rlonu(i) 
+   boundslon_reg(i,west)=rlonu(i+1) 
+  ENDDO
+
+  boundslat_reg(1,north)= PI/2 
+  boundslat_reg(1,south)= rlatv(1)
+  DO j=2,jj
+   boundslat_reg(i,north)=rlatv(j-1) 
+   boundslat_reg(i,south)=rlatv(j) 
+  ENDDO
+  boundslat_reg(jj+1,north)= rlatv(jj) 
+  boundslat_reg(jj+1,south)= -PI/2
+  
+  CALL init_regular_lonlat(ii,jj+1, rlonv(1:ii), rlatu, boundslon_reg, boundslat_reg)
+  
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   
   ! Generate global arrays on full physics grid
+  ALLOCATE(latfi_glo(klon_glo),lonfi_glo(klon_glo),cufi_glo(klon_glo),cvfi_glo(klon_glo))
+  ALLOCATE(airefi_glo(klon_glo))
+  ALLOCATE(boundslonfi_glo(klon_glo,4))
+  ALLOCATE(boundslatfi_glo(klon_glo,4))
+ 
   ALLOCATE(latfi(klon_glo),lonfi(klon_glo),cufi(klon_glo),cvfi(klon_glo))
   ALLOCATE(airefi(klon_glo))
+  ALLOCATE(boundslonfi(klon_glo,4))
+  ALLOCATE(boundslatfi(klon_glo,4))
 
   IF (klon_glo>1) THEN ! general case
     ! North pole
-    latfi(1)=rlatu(1)
-    lonfi(1)=0.
-    cufi(1) = cu(1)
-    cvfi(1) = cv(1)
+    latfi_glo(1)=rlatu(1)
+    lonfi_glo(1)=0.
+    cufi_glo(1) = cu(1)
+    cvfi_glo(1) = cv(1)
+    boundslonfi_glo(1,north_east)=0
+    boundslatfi_glo(1,north_east)=PI/2
+    boundslonfi_glo(1,north_west)=2*PI
+    boundslatfi_glo(1,north_west)=PI/2
+    boundslonfi_glo(1,south_west)=2*PI
+    boundslatfi_glo(1,south_west)=rlatv(1)
+    boundslonfi_glo(1,south_east)=0
+    boundslatfi_glo(1,south_east)=rlatv(1)
     DO j=2,jj
       DO i=1,ii
-        latfi((j-2)*ii+1+i)= rlatu(j)
-        lonfi((j-2)*ii+1+i)= rlonv(i)
-        cufi((j-2)*ii+1+i) = cu((j-1)*ii+1+i)
-        cvfi((j-2)*ii+1+i) = cv((j-1)*ii+1+i)
+        k=(j-2)*ii+1+i
+        latfi_glo(k)= rlatu(j)
+        lonfi_glo(k)= rlonv(i)
+        cufi_glo(k) = cu((j-1)*ii+1+i)
+        cvfi_glo(k) = cv((j-1)*ii+1+i)
+        boundslonfi_glo(k,north_east)=rlonu(i)
+        boundslatfi_glo(k,north_east)=rlatv(j-1)
+        boundslonfi_glo(k,north_west)=rlonu(i+1)
+        boundslatfi_glo(k,north_west)=rlatv(j-1)
+        boundslonfi_glo(k,south_west)=rlonu(i+1)
+        boundslatfi_glo(k,south_west)=rlatv(j)
+        boundslonfi_glo(k,south_east)=rlonu(i)
+        boundslatfi_glo(k,south_east)=rlatv(j)
       ENDDO
     ENDDO
     ! South pole
-    latfi(klon_glo)= rlatu(jj+1)
-    lonfi(klon_glo)= 0.
-    cufi(klon_glo) = cu((ii+1)*jj+1)
-    cvfi(klon_glo) = cv((ii+1)*jj-ii)
+    latfi_glo(klon_glo)= rlatu(jj+1)
+    lonfi_glo(klon_glo)= 0.
+    cufi_glo(klon_glo) = cu((ii+1)*jj+1)
+    cvfi_glo(klon_glo) = cv((ii+1)*jj-ii)
+    boundslonfi_glo(klon_glo,north_east)=0
+    boundslatfi_glo(klon_glo,north_east)=rlatv(jj)
+    boundslonfi_glo(klon_glo,north_west)=2*PI
+    boundslatfi_glo(klon_glo,north_west)=rlatv(jj)
+    boundslonfi_glo(klon_glo,south_west)=2*PI
+    boundslatfi_glo(klon_glo,south_west)=-PI/2
+    boundslonfi_glo(klon_glo,south_east)=rlonu(0)
+    boundslatfi_glo(klon_glo,south_east)=-Pi/2
 
     ! build airefi(), mesh area on physics grid
-    CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,aire,airefi)
+    CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,aire,airefi_glo)
     ! Poles are single points on physics grid
-    airefi(1)=sum(aire(1:ii,1))
-    airefi(klon_glo)=sum(aire(1:ii,jj+1))
+    airefi_glo(1)=sum(aire(1:ii,1))
+    airefi_glo(klon_glo)=sum(aire(1:ii,jj+1))
 
     ! Sanity check: do total planet area match between physics and dynamics?
     total_area_dyn=sum(aire(1:ii,1:jj+1))
-    total_area_phy=sum(airefi(1:klon_glo))
+    total_area_phy=sum(airefi_glo(1:klon_glo))
     IF (total_area_dyn/=total_area_phy) THEN
       WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
@@ -132,32 +192,38 @@
   ELSE ! klon_glo==1, running the 1D model
     ! just copy over input values
-    latfi(1)=rlatu(1)
-    lonfi(1)=rlonv(1)
-    cufi(1)=cu(1)
-    cvfi(1)=cv(1)
-    airefi(1)=aire(1,1)
+    latfi_glo(1)=rlatu(1)
+    lonfi_glo(1)=rlonv(1)
+    cufi_glo(1)=cu(1)
+    cvfi_glo(1)=cv(1)
+    airefi_glo(1)=aire(1,1)
+    boundslonfi_glo(1,north_east)=rlonu(1)
+    boundslatfi_glo(1,north_east)=PI/2
+    boundslonfi_glo(1,north_west)=rlonu(2)
+    boundslatfi_glo(1,north_west)=PI/2
+    boundslonfi_glo(1,south_west)=rlonu(2)
+    boundslatfi_glo(1,south_west)=rlatv(1)
+    boundslonfi_glo(1,south_east)=rlonu(1)
+    boundslatfi_glo(1,south_east)=rlatv(1)
   ENDIF ! of IF (klon_glo>1)
 
 !$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/)
   ! Now generate local lon/lat/cu/cv/area arrays
-  CALL initcomgeomphy(klon_omp)
+
 
   offset = klon_mpi_begin - 1
-  airephy(1:klon_omp) = airefi(offset+klon_omp_begin:offset+klon_omp_end)
-  cuphy(1:klon_omp) = cufi(offset+klon_omp_begin:offset+klon_omp_end)
-  cvphy(1:klon_omp) = cvfi(offset+klon_omp_begin:offset+klon_omp_end)
-  rlond(1:klon_omp) = lonfi(offset+klon_omp_begin:offset+klon_omp_end)
-  rlatd(1:klon_omp) = latfi(offset+klon_omp_begin:offset+klon_omp_end)
+  airefi(1:klon_omp) = airefi_glo(offset+klon_omp_begin:offset+klon_omp_end)
+  cufi(1:klon_omp) = cufi_glo(offset+klon_omp_begin:offset+klon_omp_end)
+  cvfi(1:klon_omp) = cvfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
+  lonfi(1:klon_omp) = lonfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
+  latfi(1:klon_omp) = latfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
+  boundslonfi(1:klon_omp,:) = boundslonfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
+  boundslatfi(1:klon_omp,:) = boundslatfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
 
   ! copy over global grid longitudes and latitudes
-  CALL initcomgeomphy_horiz(iim,jjm,rlonu,rlonv,rlatu,rlatv)
+  CALL init_geometry(lonfi, latfi, boundslonfi, boundslatfi, airefi, cufi, cvfi)
+
   
   ! copy over preff , ap(), bp(), etc 
-  CALL initcomgeomphy_vert(nlayer,preff,ap,bp,presnivs,pseudoalt)
-
-!    ! suphel => initialize some physical constants (orbital parameters,
-!    !           geoid, gravity, thermodynamical constants, etc.) in the
-!    !           physics
-!  CALL suphel
+  CALL init_vertical_layers(nlayer,preff,ap,bp,presnivs,pseudoalt)
 
   ! Initialize tracer names, numbers, etc. for physics
@@ -179,5 +245,5 @@
 !!$OMP PARALLEL
   IF (iflag_phys>=100) THEN
-    CALL iniaqua(klon_omp, rlatd, rlond, iflag_phys)
+    CALL iniaqua(klon_omp, iflag_phys)
   END IF
 !$OMP END PARALLEL
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/add_phys_tend.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/add_phys_tend.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/add_phys_tend.F90	(revision 3825)
@@ -18,4 +18,5 @@
 use phys_state_var_mod
 USE inifis_mod, ONLY: prt_level
+USE geometry_mod, ONLY : lon_degrees, lat_degrees
 IMPLICIT none
   include "YOMCST.h"
@@ -106,5 +107,5 @@
          i=jadrs(j)
          if(prt_level.ge.debug_level) THEN
-          print*,'PLANTAGE POUR LE POINT i rlon rlat =',i,rlon(i),rlat(i),text
+          print*,'PLANTAGE POUR LE POINT i lon lat =',i,lon_degrees(i),lat_degrees(i),text
           print*,'l    T     dT       Q     dQ    '
           DO k = 1, klev
@@ -124,5 +125,5 @@
         i=jqadrs(j)
           if(prt_level.ge.debug_level) THEN
-           print*,'WARNING  : EAU POUR LE POINT i rlon rlat =',i,rlon(i),rlat(i),text
+           print*,'WARNING  : EAU POUR LE POINT i lon lat =',i,lon_degrees(i),lat_degrees(i),text
            print*,'l    T     dT       Q     dQ    '
            DO k = 1, klev
@@ -205,5 +206,5 @@
          k=kadrs(j)
          if(prt_level.ge.debug_level) THEN
-          print*,'PLANTAGE2 POUR LE POINT i itap rlon rlat txt jbad zdt t',i,itap,rlon(i),rlat(i),text,jbad, &
+          print*,'PLANTAGE2 POUR LE POINT i itap lon lat txt jbad zdt t',i,itap,lon_degrees(i),lat_degrees(i),text,jbad, &
        &        zdt(i,k),t_seri(i,k)-zdt(i,k)
 !!!       if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN 
@@ -222,5 +223,5 @@
          k=kqadrs(j)
          if(prt_level.ge.debug_level) THEN
-          print*,'WARNING  : EAU2 POUR LE POINT i itap rlon rlat txt jqbad zdq q zdql ql',i,itap,rlon(i),rlat(i),text,jqbad,&
+          print*,'WARNING  : EAU2 POUR LE POINT i itap lon lat txt jqbad zdq q zdql ql',i,itap,lon_degrees(i), lat_degrees(i), text,jqbad,&
        &        zdq(i,k), q_seri(i,k)-zdq(i,k), zdql(i,k), ql_seri(i,k)-zdql(i,k)
 !!!       if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN 
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/calcul_divers.h
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/calcul_divers.h	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/calcul_divers.h	(revision 3825)
@@ -8,5 +8,5 @@
          DO i=1, klon
             IF(pctsrf(i,is_ter).GT.0.) THEN
-               paire_ter(i)=airephy(i)*pctsrf(i,is_ter)
+               paire_ter(i)=cell_area(i)*pctsrf(i,is_ter)
             ENDIF 
          ENDDO
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/carbon_cycle_mod.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/carbon_cycle_mod.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/carbon_cycle_mod.F90	(revision 3825)
@@ -86,5 +86,5 @@
 
     USE dimphy
-    USE comgeomphy
+    USE geometry_mod, ONLY : cell_area 
     USE mod_phys_lmdz_transfert_para
     USE infotrac_phy
@@ -289,5 +289,5 @@
 ! 5) Calculate total area of the earth surface
 ! --------------------------------------------
-    CALL reduce_sum(SUM(airephy),airetot)
+    CALL reduce_sum(SUM(cell_area),airetot)
     CALL bcast(airetot)
 
@@ -308,7 +308,7 @@
     USE phys_cal_mod, ONLY : mth_cur, mth_len
     USE phys_cal_mod, ONLY : day_cur
-    USE comgeomphy
     USE indice_sol_mod
     USE inifis_mod, ONLY: lunout
+    USE geometry_mod, ONLY : cell_area
 
     IMPLICIT NONE
@@ -443,5 +443,5 @@
 
           ! Calculate a global mean value of delta CO2 flux
-          fco2_tmp(1:klon) = fco2_tmp(1:klon) * airephy(1:klon)
+          fco2_tmp(1:klon) = fco2_tmp(1:klon) * cell_area(1:klon)
           CALL reduce_sum(SUM(fco2_tmp),sumtmp)
           CALL bcast(sumtmp)
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/comgeomphy.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/comgeomphy.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/comgeomphy.F90	(revision 3825)
@@ -1,24 +1,24 @@
 module comgeomphy
-   real,save,allocatable :: airephy(:)
-   real,save,allocatable :: cuphy(:)
-   real,save,allocatable :: cvphy(:)
-   real,save,allocatable :: rlatd(:)
-   real,save,allocatable :: rlond(:)
-!$OMP THREADPRIVATE(airephy,cuphy,cvphy,rlatd,rlond)
+!   real,save,allocatable :: airephy(:) => transfered in geometry_mod (area)
+!   real,save,allocatable :: cuphy(:) => transfered in regular_lonlat_mod : dlon_reg
+!   real,save,allocatable :: cvphy(:) => transfered in regular_lonlat_mod : dlat_reg
+!   real,save,allocatable :: rlatd(:) => transfered in geometry_mod (longitude) 
+!   real,save,allocatable :: rlond(:) => transfered in geometry_mod (latitude) 
+!!$OMP THREADPRIVATE(airephy,rlatd,rlond)
    ! general (over global grid) latitudes and longitudes:
    ! useful for outputs
-   real,save,allocatable :: rlonv(:) ! longitudes (rad)
-   real,save,allocatable :: rlonu(:) ! longitude boundaries (rad)
-   real,save,allocatable :: rlatu(:) ! latitudes (rad)
-   real,save,allocatable :: rlatv(:) ! latitude boundaries (rad)
-!$OMP THREADPRIVATE(rlonv,rlonu,rlatu,rlatv)
+!   real,save,allocatable :: rlonv(:) ! longitudes (rad)           => transfered in regular_lonlat_mod 
+!   real,save,allocatable :: rlonu(:) ! longitude boundaries (rad) => transfered in regular_lonlat_mod 
+!   real,save,allocatable :: rlatu(:) ! latitudes (rad)            => transfered in regular_lonlat_mod
+!   real,save,allocatable :: rlatv(:) ! latitude boundaries (rad)  => transfered in regular_lonlat_mod 
+!!$OMP THREADPRIVATE(rlonv,rlonu,rlatu,rlatv)
 
    ! vertical grid
-   real,save :: preff  ! reference surface pressure (Pa)
-   real,save,allocatable :: ap(:)
-   real,save,allocatable :: bp(:)
-   real,save,allocatable :: presnivs(:)
-   real,save,allocatable :: pseudoalt(:)
-!$OMP THREADPRIVATE(preff,ap,bp,presnivs,pseudoalt)
+!   real,save :: preff  ! reference surface pressure (Pa)
+!   real,save,allocatable :: ap(:)        => transfered in vertical_mod
+!   real,save,allocatable :: bp(:)        => transfered in vertical_mod
+!   real,save,allocatable :: presnivs(:)  => transfered in vertical_mod
+!   real,save,allocatable :: pseudoalt(:) => transfered in vertical_mod
+!!$OMP THREADPRIVATE(preff,ap,bp,presnivs,pseudoalt)
 
 contains
@@ -28,32 +28,16 @@
     integer,intent(in) :: klon_omp ! number of columns (on local omp grid)
  
-    allocate(airephy(klon_omp))
-    allocate(cuphy(klon_omp))
-    allocate(cvphy(klon_omp))
-    allocate(rlatd(klon_omp))
-    allocate(rlond(klon_omp))
+!    allocate(airephy(klon_omp))
+!    allocate(rlatd(klon_omp))
+!    allocate(rlond(klon_omp))
 
   end subroutine InitComgeomphy
   
-  subroutine initcomgeomphy_horiz(iim,jjm,rlonu_dyn,rlonv_dyn,&
-                                  rlatu_dyn,rlatv_dyn)
+  subroutine initcomgeomphy_horiz(iim,jjm)
     IMPLICIT NONE
     integer,intent(in) :: iim
     integer,intent(in) :: jjm
-    real,intent(in) :: rlonu_dyn(iim+1) ! dyn grid boundaries (rad)
-    real,intent(in) :: rlonv_dyn(iim+1) ! dyn grid longitudes (rad)
-    real,intent(in) :: rlatu_dyn(jjm+1) ! dyn grid latitudes (rad)
-    real,intent(in) :: rlatv_dyn(jjm)  ! dyn grid boundaries (rad)
     
-    allocate(rlonu(iim+1))
-    allocate(rlonv(iim+1))
-    allocate(rlatu(jjm+1))
-    allocate(rlatv(jjm))
-    
-    rlonu(:)=rlonu_dyn(:)
-    rlonv(:)=rlonv_dyn(:)
-    rlatu(:)=rlatu_dyn(:)
-    rlatv(:)=rlatv_dyn(:)
-    
+   
   end subroutine initcomgeomphy_horiz
   
@@ -68,14 +52,14 @@
     real,intent(in) :: pseudoalt_dyn(nlayer) ! pseudo-altitude of atm. layers (km)
   
-    allocate(ap(nlayer+1))
-    allocate(bp(nlayer+1))
-    allocate(presnivs(nlayer))
-    allocate(pseudoalt(nlayer))
-  
-    preff=preff_dyn
-    ap(:)=ap_dyn(:)
-    bp(:)=bp_dyn(:)
-    presnivs(:)=presnivs_dyn(:)
-    pseudoalt(:)=pseudoalt_dyn(:)
+!    allocate(ap(nlayer+1))
+!    allocate(bp(nlayer+1))
+!    allocate(presnivs(nlayer))
+!    allocate(pseudoalt(nlayer))
+!  
+!    preff=preff_dyn
+!    ap(:)=ap_dyn(:)
+!    bp(:)=bp_dyn(:)
+!    presnivs(:)=presnivs_dyn(:)
+!    pseudoalt(:)=pseudoalt_dyn(:)
 
   end subroutine initcomgeomphy_vert
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/cpl_mod.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/cpl_mod.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/cpl_mod.F90	(revision 3825)
@@ -294,5 +294,5 @@
 ! are stored in this module.
     USE surface_data
-    USE phys_state_var_mod, ONLY : rlon, rlat
+    USE geometry_mod, ONLY : lon_degrees, lat_degrees
     USE carbon_cycle_mod, ONLY : carbon_cycle_cpl
     USE indice_sol_mod
@@ -365,6 +365,6 @@
 
 ! Transform the longitudes and latitudes on 2D arrays
-          CALL gather_omp(rlon,rlon_mpi)
-          CALL gather_omp(rlat,rlat_mpi)
+          CALL gather_omp(lon_degrees,rlon_mpi)
+          CALL gather_omp(lat_degrees,rlat_mpi)
 !$OMP MASTER
           CALL Grid1DTo2D_mpi(rlon_mpi,tmp_lon)
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/geometry_mod.f90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/geometry_mod.f90	(revision 3825)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/geometry_mod.f90	(revision 3825)
@@ -0,0 +1,70 @@
+MODULE geometry_mod
+
+
+  REAL,SAVE,ALLOCATABLE :: longitude(:) ! longitude of the cell (rad)
+!$OMP THREADPRIVATE(longitude)
+
+  REAL,SAVE,ALLOCATABLE :: latitude(:)! latitude of the cell (rad)
+!$OMP THREADPRIVATE(latitude)
+
+  REAL,SAVE,ALLOCATABLE :: lon_degrees(:) ! longitude of the cell (degree)
+!$OMP THREADPRIVATE(lon_degrees)
+
+  REAL,SAVE,ALLOCATABLE :: lat_degrees(:)! latitude of the cell (degree)
+!$OMP THREADPRIVATE(lat_degrees)
+
+  REAL,SAVE,ALLOCATABLE :: boundslon(:,:)  ! boundaries of the cell (rad)
+!$OMP THREADPRIVATE(boundslon)
+
+  REAL,SAVE,ALLOCATABLE :: boundslat(:,:) ! boundaries of the cell (rad)
+!$OMP THREADPRIVATE(boundslat)
+
+  REAL,SAVE,ALLOCATABLE :: dx(:)      ! resolution of longitude cell (valid only for 2D grid)
+!$OMP THREADPRIVATE(dx)
+  
+  REAL,SAVE,ALLOCATABLE :: dy(:)      ! resolution of latitude cell (valid only for 2D grid)
+!$OMP THREADPRIVATE(dy)
+
+  REAL,SAVE,ALLOCATABLE :: cell_area(:)      ! area of the cell
+!$OMP THREADPRIVATE(cell_area)
+
+
+CONTAINS
+
+  SUBROUTINE init_geometry(longitude_, latitude_, boundslon_, boundslat_, cell_area_, dx_, dy_)
+  USE dimphy, ONLY: klon
+  USE mod_grid_phy_lmdz, ONLY: nvertex
+  USE nrtype, ONLY : PI
+  IMPLICIT NONE
+    REAL :: longitude_(klon)
+    REAL :: latitude_(klon)
+    REAL :: boundslon_(klon,nvertex)
+    REAL :: boundslat_(klon,nvertex)
+    REAL :: cell_area_(klon)
+    REAL,OPTIONAL :: dx_(klon)
+    REAL,OPTIONAL :: dy_(klon)
+    
+    ALLOCATE(longitude(klon))
+    ALLOCATE(latitude(klon))
+    ALLOCATE(lon_degrees(klon))
+    ALLOCATE(lat_degrees(klon))
+    ALLOCATE(boundslon(klon,nvertex))
+    ALLOCATE(boundslat(klon,nvertex))
+    ALLOCATE(cell_area(klon))
+    IF (PRESENT(dx_)) ALLOCATE(dx(klon))
+    IF (PRESENT(dy_))ALLOCATE(dy(klon))
+    
+    longitude(:) = longitude_(:)
+    latitude(:) = latitude_(:)
+    lon_degrees(:) = longitude(:)*180./PI
+    lat_degrees(:) = latitude(:)*180./PI
+    boundslon(:,:) = boundslon_(:,:)
+    boundslat(:,:) = boundslat_(:,:)
+    cell_area(:) = cell_area_(:)
+    IF (PRESENT(dx_)) dx(:) = dx_(:)
+    IF (PRESENT(dy_)) dy(:) = dy_(:)
+    
+  END SUBROUTINE init_geometry
+
+
+END MODULE geometry_mod
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/hgardfou.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/hgardfou.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/hgardfou.F90	(revision 3825)
@@ -6,4 +6,5 @@
   USE indice_sol_mod
   USE inifis_mod, ONLY: lunout
+  USE geometry_mod, ONLY : lon_degrees, lat_degrees 
   IMPLICIT NONE
   ! ======================================================================
@@ -30,5 +31,5 @@
     firstcall = .FALSE.
     ! DO i = 1, klon
-    ! WRITE(lunout,*)'i=',i,'rlon=',rlon(i),'rlat=',rlat(i)
+    ! WRITE(lunout,*)'i=',i,'lon=',lon_degrees(i),'lat=',lat_degrees(i)
     ! ENDDO
 
@@ -55,5 +56,5 @@
       DO i = 1, jbad
         WRITE (lunout, *) 'i,k,temperature,lon,lat,pourc ter,lic,oce,sic =', &
-          jadrs(i), k, zt(jadrs(i)), rlon(jadrs(i)), rlat(jadrs(i)), &
+          jadrs(i), k, zt(jadrs(i)), lon_degrees(jadrs(i)), lat_degrees(jadrs(i)), &
           (pctsrf(jadrs(i),nsrf), nsrf=1, nbsrf)
       END DO
@@ -75,5 +76,5 @@
       DO i = 1, jbad
         WRITE (lunout, *) 'i,k,temperature,lon,lat,pourc ter,lic,oce,sic =', &
-          jadrs(i), k, zt(jadrs(i)), rlon(jadrs(i)), rlat(jadrs(i)), &
+          jadrs(i), k, zt(jadrs(i)), lon_degrees(jadrs(i)), lat_degrees(jadrs(i)), &
           (pctsrf(jadrs(i),nsrf), nsrf=1, nbsrf)
       END DO
@@ -101,5 +102,5 @@
         WRITE (lunout, *) &
           'i,nsrf,temperature,lon,lat,pourc ter,lic,oce,sic =', jadrs(i), &
-          nsrf, zt(jadrs(i)), rlon(jadrs(i)), rlat(jadrs(i)), &
+          nsrf, zt(jadrs(i)), lon_degrees(jadrs(i)), lat_degrees(jadrs(i)), &
           pctsrf(jadrs(i), nsrf)
       END DO
@@ -122,5 +123,5 @@
         WRITE (lunout, *) &
           'i,nsrf,temperature,lon,lat,pourc ter,lic,oce,sic =', jadrs(i), &
-          nsrf, zt(jadrs(i)), rlon(jadrs(i)), rlat(jadrs(i)), &
+          nsrf, zt(jadrs(i)), lon_degrees(jadrs(i)), lat_degrees(jadrs(i)), &
           pctsrf(jadrs(i), nsrf)
       END DO
Index: namico_lmdz/aquaplanet/LMDZ5/libf/phylmd/ini_coord_REGDYN.h
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/ini_coord_REGDYN.h	(revision 3824)
+++ 	(revision )
@@ -1,81 +1,0 @@
-c
-c $Header$
-c
-       nsrf=3
-       DO nreg=1, nbregdyn
-       DO i=1, klon
-
-c       IF (debut) THEN
-         IF(rlon(i).LT.0.) THEN
-           rlonPOS(i)=rlon(i)+360.
-         ELSE
-           rlonPOS(i)=rlon(i)  
-         ENDIF
-c       ENDIF
-
-        pct_ocean(i,nreg)=0
-
-c test si c'est 1 point d'ocean
-        IF(pctsrf(i,nsrf).EQ.1.) THEN
-
-         IF(nreg.EQ.1) THEN
-
-c TROP
-          IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
-           pct_ocean(i,nreg)=1
-          ENDIF
-
-c PACIFIQUE NORD
-          ELSEIF(nreg.EQ.2) THEN
-           IF(rlat(i).GE.40.AND.rlat(i).LE.60.) THEN
-            IF(rlonPOS(i).GE.160..AND.rlonPOS(i).LE.235.) THEN 
-             pct_ocean(i,nreg)=1
-            ENDIF
-           ENDIF
-c CALIFORNIE ST-CU
-         ELSEIF(nreg.EQ.3) THEN
-          IF(rlonPOS(i).GE.220..AND.rlonPOS(i).LE.250.) THEN
-           IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
-            pct_ocean(i,nreg)=1
-           ENDIF
-          ENDIF
-c HAWAI
-        ELSEIF(nreg.EQ.4) THEN 
-         IF(rlonPOS(i).GE.180..AND.rlonPOS(i).LE.220.) THEN
-          IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
-           pct_ocean(i,nreg)=1
-          ENDIF
-         ENDIF
-c WARM POOL
-        ELSEIF(nreg.EQ.5) THEN 
-         IF(rlonPOS(i).GE.70..AND.rlonPOS(i).LE.150.) THEN
-          IF(rlat(i).GE.-5.AND.rlat(i).LE.20.) THEN
-           pct_ocean(i,nreg)=1
-          ENDIF
-         ENDIF
-        ENDIF !nbregdyn
-c TROP
-c        IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
-c         pct_ocean(i)=.TRUE.
-c         WRITE(*,*) 'pct_ocean =',i, rlon(i), rlat(i)
-c          ENDIF !lon
-c         ENDIF !lat
-
-        ENDIF !pctsrf
-       ENDDO !klon
-       ENDDO !nbregdyn
-cIM 190504      ENDIF !ok_regdyn
- 
-cIM somme de toutes les nhistoW BEG
-      IF (debut) THEN
-      DO nreg = 1, nbregdyn
-       DO k = 1, kmaxm1
-        DO l = 1, lmaxm1
-         DO iw = 1, iwmax
-          nhistoWt(k,l,iw,nreg)=0.
-         ENDDO !iw
-        ENDDO !l
-       ENDDO !k
-      ENDDO !nreg
-      ENDIF !(debut) THEN
-cIM 190504 BEG
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/ini_histday_seri.h
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/ini_histday_seri.h	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/ini_histday_seri.h	(revision 3825)
@@ -13,13 +13,13 @@
          CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
 !
-         CALL grid1dTo2d_glo(rlon,zx_lon)
+         CALL grid1dTo2d_glo(lon_degrees,zx_lon)
          DO i = 1, nbp_lon
-            zx_lon(i,1) = rlon(i+1)
-            zx_lon(i,jjmp1) = rlon(i+1)
+            zx_lon(i,1) = lon_degrees(i+1)
+            zx_lon(i,jjmp1) = lon_degrees(i+1)
          ENDDO
          DO ll=1,klev
             znivsig(ll)=REAL(ll)
          ENDDO
-         CALL grid1dTo2d_glo(rlat,zx_lat)
+         CALL grid1dTo2d_glo(lat_degrees,zx_lat)
 !
          imin_debut=1 
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/ini_paramLMDZ_phy.h
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/ini_paramLMDZ_phy.h	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/ini_paramLMDZ_phy.h	(revision 3825)
@@ -1,7 +1,7 @@
 !IM    Implemente en modes sequentiel et parallele
 
-       CALL gather(rlat,rlat_glo)
+       CALL gather(lat_degrees,rlat_glo)
        CALL bcast(rlat_glo)
-       CALL gather(rlon,rlon_glo)
+       CALL gather(lon_degrees,rlon_glo)
        CALL bcast(rlon_glo)
 
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/init_be.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/init_be.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/init_be.F90	(revision 3825)
@@ -5,7 +5,7 @@
 
   USE dimphy
-  USE comgeomphy
   USE infotrac_phy, ONLY : nbtr
   USE indice_sol_mod
+  USE geometry_mod, ONLY : longitude, latitude
     
   IMPLICIT NONE 
@@ -114,8 +114,8 @@
 
   DO i = 1,klon
-     qcos=sin(glt)*sin(rlatd(i))
+     qcos=sin(glt)*sin(latitude(i))
 !!jyg
 !!     qcos=qcos+cos(glt)*cos(rlatd(i))*cos(rlond(i)+glg)
-     qcos=qcos+cos(glt)*cos(rlatd(i))*cos(rlond(i)-glg)
+     qcos=qcos+cos(glt)*cos(latitude(i))*cos(longitude(i)-glg)
 !!jyg end
      IF ( qcos .LT. -1.) qcos = -1.
Index: namico_lmdz/aquaplanet/LMDZ5/libf/phylmd/init_phys_lmdz.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/init_phys_lmdz.F90	(revision 3824)
+++ 	(revision )
@@ -1,37 +1,0 @@
-!
-!$Header$
-!
-SUBROUTINE Init_Phys_lmdz(iim,jjp1,llm,nb_proc,distrib,communicator)
-  USE mod_phys_lmdz_para, ONLY: Init_phys_lmdz_para, klon_omp
-  USE mod_grid_phy_lmdz, ONLY: Init_grid_phy_lmdz, nbp_lev
-  USE dimphy, ONLY : Init_dimphy
-  USE infotrac_phy, ONLY : type_trac
-#ifdef REPROBUS
-  USE CHEM_REP, ONLY : Init_chem_rep_phys
-#endif
-
-  IMPLICIT NONE
-  
-    INTEGER,INTENT(in) :: iim
-    INTEGER,INTENT(in) :: jjp1
-    INTEGER,INTENT(in) :: llm
-    INTEGER,INTENT(in) :: nb_proc
-    INTEGER,INTENT(in) :: distrib(0:nb_proc-1)
-    INTEGER,INTENT(in) :: communicator
-
-
-    CALL Init_grid_phy_lmdz(iim,jjp1,llm)
-    CALL Init_phys_lmdz_para(iim,jjp1,nb_proc,distrib,communicator)
-!$OMP PARALLEL
-    CALL Init_dimphy(klon_omp,nbp_lev)
-
-! Initialization of Reprobus
-    IF (type_trac == 'repr') THEN
-#ifdef REPROBUS
-       CALL Init_chem_rep_phys(klon_omp,nbp_lev)
-#endif
-    END IF
-
-!$OMP END PARALLEL
- 
-END SUBROUTINE Init_Phys_lmdz  
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/lsc_scav.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/lsc_scav.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/lsc_scav.F90	(revision 3825)
@@ -10,5 +10,4 @@
   USE traclmdz_mod
   USE infotrac_phy,ONLY : nbtr
-  USE comgeomphy
   USE iophy
   IMPLICIT NONE 
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/mod_grid_phy_lmdz.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/mod_grid_phy_lmdz.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/mod_grid_phy_lmdz.F90	(revision 3825)
@@ -1,5 +1,2 @@
-!
-!$Id $
-!
 MODULE mod_grid_phy_lmdz
 
@@ -8,8 +5,14 @@
              grid2dTo1d_glo_igen, grid2dTo1d_glo_rgen, grid2dTo1d_glo_lgen
   
+  INTEGER,PARAMETER :: unstructured=0
+  INTEGER,PARAMETER :: regular_lonlat=1
+
+  INTEGER,SAVE :: grid_type
+  INTEGER,SAVE :: nvertex
   INTEGER,SAVE :: nbp_lon  ! == iim
   INTEGER,SAVE :: nbp_lat  ! == jjmp1
   INTEGER,SAVE :: nbp_lev  ! == llm
   INTEGER,SAVE :: klon_glo
+
 
   INTERFACE grid1dTo2d_glo
@@ -32,16 +35,20 @@
 
 
-  SUBROUTINE Init_grid_phy_lmdz(iim,jjp1,llm)
+  SUBROUTINE init_grid_phy_lmdz(grid_type_,nvertex_,nbp_lon_,nbp_lat_,nbp_lev_)
   IMPLICIT NONE
-  INTEGER, INTENT(in) :: iim
-  INTEGER, INTENT(in) :: jjp1
-  INTEGER, INTENT(in) :: llm
-  
-    nbp_lon=iim
-    nbp_lat=jjp1
-    nbp_lev=llm
-    klon_glo=(iim*jjp1)-2*(iim-1)
-  
-  END SUBROUTINE Init_grid_phy_lmdz
+  INTEGER,INTENT(IN)  :: grid_type_ 
+  INTEGER,INTENT(IN)  :: nvertex_ 
+  INTEGER, INTENT(IN) :: nbp_lon_
+  INTEGER, INTENT(IN) :: nbp_lat_
+  INTEGER, INTENT(IN) :: nbp_lev_
+  
+    grid_type = grid_type_
+    nvertex   = nvertex_
+    nbp_lon   = nbp_lon_
+    nbp_lat   = nbp_lat_
+    nbp_lev   = nbp_lev_
+    klon_glo  = (nbp_lon*nbp_lat)-2*(nbp_lon-1)
+  
+  END SUBROUTINE init_grid_phy_lmdz
   
   
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/mod_phys_lmdz_mpi_data.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/mod_phys_lmdz_mpi_data.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/mod_phys_lmdz_mpi_data.F90	(revision 3825)
@@ -1,5 +1,2 @@
-!
-!$Header$
-!
 MODULE mod_phys_lmdz_mpi_data
   
@@ -46,15 +43,15 @@
 CONTAINS
   
-  SUBROUTINE Init_phys_lmdz_mpi_data(iim,jjp1,nb_proc,distrib,communicator)
+  SUBROUTINE init_phys_lmdz_mpi_data(nbp, nbp_lon, nbp_lat, communicator)
   IMPLICIT NONE
 #ifdef CPP_MPI
     INCLUDE 'mpif.h'
 #endif
-    INTEGER,INTENT(in) :: iim
-    INTEGER,INTENT(in) :: jjp1
-    INTEGER,INTENT(in) :: nb_proc
-    INTEGER,INTENT(in) :: distrib(0:nb_proc-1)
-    INTEGER,INTENT(in) :: communicator
-    
+    INTEGER,INTENT(IN) :: nbp
+    INTEGER,INTENT(IN) :: nbp_lon
+    INTEGER,INTENT(IN) :: nbp_lat
+    INTEGER,INTENT(IN) :: communicator
+    
+    INTEGER,ALLOCATABLE :: distrib(:)
     INTEGER :: ierr
     INTEGER :: klon_glo
@@ -67,8 +64,8 @@
 #endif
     
-    if (iim.eq.1) then
+    if (nbp_lon.eq.1) then
        klon_glo=1
     else
-       klon_glo=iim*(jjp1-2)+2
+       klon_glo=nbp_lon*(nbp_lat-2)+2
     endif
     
@@ -85,5 +82,16 @@
       mpi_rank=0
     ENDIF
-    
+
+    ALLOCATE(distrib(0:mpi_size-1))
+
+    IF (is_using_mpi) THEN    
+#ifdef CPP_MPI
+    CALL MPI_ALLGATHER(nbp,1,MPI_REAL8,distrib,1,MPI_REAL_LMDZ,COMM_LMDZ_PHY)
+#endif
+    ELSE
+     distrib(:)=nbp
+    ENDIF
+    
+        
     IF (mpi_rank == 0) THEN
       mpi_master = 0
@@ -119,5 +127,5 @@
   
       
-    klon_mpi_para_nb(0:mpi_size-1)=distrib(0:nb_proc-1)
+    klon_mpi_para_nb(0:mpi_size-1)=distrib(0:mpi_size-1)
 
     DO i=0,mpi_size-1
@@ -136,14 +144,14 @@
         ij_para_begin(i) = 1
       ELSE
-        ij_para_begin(i) = klon_mpi_para_begin(i)+iim-1
+        ij_para_begin(i) = klon_mpi_para_begin(i)+nbp_lon-1
       ENDIF
 
-      jj_para_begin(i) = (ij_para_begin(i)-1)/iim + 1
-      ii_para_begin(i) = MOD(ij_para_begin(i)-1,iim) + 1
+      jj_para_begin(i) = (ij_para_begin(i)-1)/nbp_lon + 1
+      ii_para_begin(i) = MOD(ij_para_begin(i)-1,nbp_lon) + 1
 
       
-      ij_para_end(i) = klon_mpi_para_end(i)+iim-1
-      jj_para_end(i) = (ij_para_end(i)-1)/iim + 1
-      ii_para_end(i) = MOD(ij_para_end(i)-1,iim) + 1
+      ij_para_end(i) = klon_mpi_para_end(i)+nbp_lon-1
+      jj_para_end(i) = (ij_para_end(i)-1)/nbp_lon + 1
+      ii_para_end(i) = MOD(ij_para_end(i)-1,nbp_lon) + 1
 
 
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/mod_phys_lmdz_para.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/mod_phys_lmdz_para.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/mod_phys_lmdz_para.F90	(revision 3825)
@@ -16,13 +16,12 @@
 CONTAINS
 
-  SUBROUTINE Init_phys_lmdz_para(iim,jjp1,nb_proc,distrib,communicator)
+  SUBROUTINE Init_phys_lmdz_para(nbp, nbp_lon, nbp_lat, communicator)
   IMPLICIT NONE
-    INTEGER,INTENT(in) :: iim
-    INTEGER,INTENT(in) :: jjp1
-    INTEGER,INTENT(in) :: nb_proc
-    INTEGER,INTENT(in) :: distrib(0:nb_proc-1)
-    INTEGER,INTENT(in) :: communicator
+    INTEGER,INTENT(IN) :: nbp
+    INTEGER,INTENT(IN) :: nbp_lon
+    INTEGER,INTENT(IN) :: nbp_lat
+    INTEGER,INTENT(IN) :: communicator
 
-    CALL Init_phys_lmdz_mpi_data(iim,jjp1,nb_proc,distrib,communicator)
+    CALL Init_phys_lmdz_mpi_data(nbp, nbp_lon, nbp_lat, communicator)
 !$OMP PARALLEL
     CALL Init_phys_lmdz_omp_data(klon_mpi)
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phyaqua_mod.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phyaqua_mod.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phyaqua_mod.F90	(revision 3825)
@@ -5,5 +5,5 @@
 CONTAINS
 
-  SUBROUTINE iniaqua(nlon, latfi, lonfi, iflag_phys)
+  SUBROUTINE iniaqua(nlon, iflag_phys)
 
     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -19,30 +19,21 @@
     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-    USE comgeomphy, ONLY: rlatd, rlond
     USE dimphy, ONLY: klon
+    USE geometry_mod, ONLY : latitude
     USE surface_data, ONLY: type_ocean, ok_veget
     USE pbl_surface_mod, ONLY: pbl_surface_init
     USE fonte_neige_mod, ONLY: fonte_neige_init
     USE phys_state_var_mod
-    !USE control_phy_mod, ONLY: dayref, nday, iphysiq
     USE inifis_mod, ONLY: dayref, nday, iphysiq, dtphys, &
                           daysec, day_ini,day_end
     USE indice_sol_mod
-    !USE temps_phy_mod
-    !USE comconst_phy_mod
     USE nrtype, ONLY: pi
     USE ioipsl
     IMPLICIT NONE
 
-    ! #include "dimphy.h"
-    ! #include "YOMCST.h"
-!    include "comconst.h"
     include "clesphys.h"
     include "dimsoil.h"
-!    include "temps.h"
 
     INTEGER, INTENT (IN) :: nlon, iflag_phys
-    ! IM ajout latfi, lonfi
-    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
 
     INTEGER type_profil, type_aqua
@@ -63,6 +54,4 @@
     INTEGER it, unit, i, k, itap
 
-    REAL airefi, zcufi, zcvfi
-
     REAL rugos, albedo
     REAL tsurf
@@ -70,9 +59,5 @@
     REAL qsol_f
     REAL rugsrel(nlon)
-    ! real zmea(nlon),zstd(nlon),zsig(nlon)
-    ! real zgam(nlon),zthe(nlon),zpic(nlon),zval(nlon)
-    ! real rlon(nlon),rlat(nlon)
     LOGICAL alb_ocean
-    ! integer demih_pas
 
     CHARACTER *80 ans, file_forctl, file_fordat, file_start
@@ -123,8 +108,6 @@
     INTEGER l, ierr, aslun
 
-    REAL longitude, latitude
     REAL paire
 
-    DATA latitude, longitude/48., 0./
 
     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -154,10 +137,4 @@
     time = 0.
 
-    ! IM ajout latfi, lonfi
-    rlatd = latfi
-    rlond = lonfi
-    rlat = rlatd*180./pi
-    rlon = rlond*180./pi
-
     ! -----------------------------------------------------------------------
     ! initialisations de la physique
@@ -166,7 +143,5 @@
     day_ini = dayref
     day_end = day_ini + nday
-    airefi = 1.
-    zcufi = 1.
-    zcvfi = 1.
+
     !$OMP MASTER
     nbapp_rad_omp = 24
@@ -266,5 +241,5 @@
     END DO
     ! IM calcul profil sst
-    CALL profil_sst(nlon, rlatd, type_profil, phy_sst)
+    CALL profil_sst(nlon, latitude, type_profil, phy_sst)
 
     CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phyetat0.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phyetat0.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phyetat0.F90	(revision 3825)
@@ -9,4 +9,5 @@
   USE pbl_surface_mod,  ONLY : pbl_surface_init
   USE surface_data,     ONLY : type_ocean, version_ocean
+  USE geometry_mod,         ONLY : lon_degrees, lat_degrees
   USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, dtime, &
        qsol, fevap, z0m, z0h, agesno, &
@@ -14,5 +15,5 @@
        falb_dir, falb_dif, &
        ftsol, pbl_tke, pctsrf, q_ancien, radpas, radsol, rain_fall, ratqs, &
-       rlat, rlon, rnebcon, rugoro, sig1, snow_fall, solaire_etat0, sollw, sollwdown, &
+       rnebcon, rugoro, sig1, snow_fall, solaire_etat0, sollw, sollwdown, &
        solsw, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
        wake_deltat, wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
@@ -135,11 +136,4 @@
   clesphy0(8)=tab_cntrl( 12 )
 
-  ! Lecture des latitudes (coordonnees):
-
-  CALL get_field("latitude", rlat)
-
-  ! Lecture des longitudes (coordonnees):
-
-  CALL get_field("longitude", rlon)
 
   ! Lecture du masque terre mer
@@ -845,8 +839,8 @@
   ! Initialize module ocean_cpl_mod for the case of coupled ocean
   IF ( type_ocean == 'couple' ) THEN
-     CALL ocean_cpl_init(dtime, rlon, rlat)
-  ENDIF
-
-  CALL init_iophy_new(rlat, rlon)
+     CALL ocean_cpl_init(dtime, lon_degrees, lat_degrees)
+  ENDIF
+
+  CALL init_iophy_new(lon_degrees, lat_degrees)
 
   ! Initilialize module fonte_neige_mod      
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phyredem.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phyredem.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phyredem.F90	(revision 3825)
@@ -9,4 +9,5 @@
   USE pbl_surface_mod,  ONLY : pbl_surface_final
   USE phys_state_var_mod
+  USE geometry_mod, ONLY : lon_degrees, lat_degrees
   USE iostart
   USE traclmdz_mod, ONLY : traclmdz_to_restart
@@ -93,7 +94,7 @@
 
   CALL put_field("longitude", &
-       "Longitudes de la grille physique", rlon)
-
-  CALL put_field("latitude", "Latitudes de la grille physique", rlat)
+       "Longitudes de la grille physique", lon_degrees)
+
+  CALL put_field("latitude", "Latitudes de la grille physique", lat_degrees)
 
   ! PB ajout du masque terre/mer
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phys_output_mod.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phys_output_mod.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phys_output_mod.F90	(revision 3825)
@@ -47,5 +47,5 @@
     USE mod_grid_phy_lmdz, only: klon_glo
     !USE temps_phy_mod
-    USE comgeomphy, ONLY: ap,bp,presnivs,preff
+    USE vertical_layers_mod,ONLY: ap,bp,presnivs,preff 
     !Use comvert_phy_mod
     USE inifis_mod, ONLY: day_ini, itau_phy, start_time, annee_ref, day_ref
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phys_output_write_mod.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phys_output_write_mod.F90	(revision 3825)
@@ -235,5 +235,5 @@
     USE indice_sol_mod, only: nbsrf
     USE infotrac_phy, only: nqtot, nqo, type_trac
-    USE comgeomphy, only: airephy
+    USE geometry_mod, ONLY : cell_area
     USE surface_data, only: type_ocean, version_ocean, ok_veget, ok_snow
 !    USE aero_mod, only: naero_spc
@@ -320,5 +320,5 @@
 !!! Champs 1D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        CALL histwrite_phy(o_phis, pphis)
-       CALL histwrite_phy(o_aire, airephy)
+       CALL histwrite_phy(o_aire, cell_area)
 
        IF (vars_defined) THEN
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phys_state_var_mod.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phys_state_var_mod.F90	(revision 3825)
@@ -20,6 +20,7 @@
 !$OMP THREADPRIVATE(dtime, solaire_etat0)
 
-      REAL, ALLOCATABLE, SAVE :: rlat(:), rlon(:), pctsrf(:,:)
-!$OMP THREADPRIVATE(rlat, rlon, pctsrf)
+      REAL, ALLOCATABLE, SAVE ::  pctsrf(:,:)
+!$OMP THREADPRIVATE(pctsrf)
+
       REAL, ALLOCATABLE, SAVE :: ftsol(:,:)
 !$OMP THREADPRIVATE(ftsol)
@@ -418,5 +419,4 @@
 
 include "clesphys.h"
-      ALLOCATE(rlat(klon), rlon(klon))
       ALLOCATE(pctsrf(klon,nbsrf))
       ALLOCATE(ftsol(klon,nbsrf))
@@ -593,5 +593,5 @@
 include "clesphys.h"
 
-      deallocate(rlat, rlon, pctsrf, ftsol, falb1, falb2)
+      deallocate(pctsrf, ftsol, falb1, falb2)
       deallocate(qsol,fevap,z0m,z0h,agesno)
       deallocate(rain_fall, snow_fall, solsw, sollw, radsol)
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/physics_distribution_mod.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/physics_distribution_mod.F90	(revision 3825)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/physics_distribution_mod.F90	(revision 3825)
@@ -0,0 +1,47 @@
+MODULE physics_distribution_mod
+
+
+CONTAINS
+
+  SUBROUTINE init_physics_distribution(grid_type, nvertex, nbp, nbp_lon, nbp_lat, nbp_lev, communicator)
+  USE mod_phys_lmdz_para, ONLY: init_phys_lmdz_para, klon_omp
+  USE mod_grid_phy_lmdz, ONLY: init_grid_phy_lmdz
+  USE dimphy, ONLY : Init_dimphy
+  USE infotrac_phy, ONLY : type_trac
+#ifdef REPROBUS
+  USE CHEM_REP, ONLY : Init_chem_rep_phys
+#endif
+
+  IMPLICIT NONE
+    INTEGER,INTENT(IN) :: grid_type 
+    INTEGER,INTENT(IN) :: nvertex 
+    INTEGER,INTENT(IN) :: nbp           
+    INTEGER,INTENT(IN) :: nbp_lon
+    INTEGER,INTENT(IN) :: nbp_lat
+    INTEGER,INTENT(IN) :: nbp_lev
+    INTEGER,INTENT(IN) :: communicator
+
+
+    CALL init_grid_phy_lmdz(grid_type,nvertex, nbp_lon,nbp_lat,nbp_lev)
+    CALL init_phys_lmdz_para(nbp,nbp_lon, nbp_lat, communicator)
+!$OMP PARALLEL
+    CALL init_dimphy(klon_omp,nbp_lev)
+
+! Initialization of Reprobus
+    IF (type_trac == 'repr') THEN
+#ifdef REPROBUS
+       CALL Init_chem_rep_phys(klon_omp,nbp_lev)
+#endif
+    END IF
+
+!$OMP END PARALLEL
+ 
+END SUBROUTINE init_physics_distribution  
+
+
+
+
+
+
+
+END MODULE physics_distribution_mod
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/physiq.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/physiq.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/physiq.F90	(revision 3825)
@@ -7,10 +7,9 @@
      u,v,t,qx, &
      flxmass_w, &
-     d_u, d_v, d_t, d_qx, d_ps &
-     , dudyn)
+     d_u, d_v, d_t, d_qx, d_ps)
 
   USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
        histwrite, ju2ymds, ymds2ju, getin
-  USE comgeomphy, ONLY: airephy
+  USE geometry_mod, ONLY: cell_area
   USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, year_cur, &
        mth_cur, phys_cal_update
@@ -33,4 +32,5 @@
   USE phys_output_ctrlout_mod
   USE iophy
+  USE geometry_mod, ONLY : lon_degrees, lat_degrees
   use open_climoz_m, only: open_climoz ! ozone climatology from a file
   use regr_pr_av_m, only: regr_pr_av
@@ -681,9 +681,5 @@
 !  PARAMETER(imp1jmp1=(iim+1)*jjmp1) ! => (imp1jmp1=(nbp_lon+1)*(nbp_lat-1/(nbp_lat-1)))
   !ym A voir plus tard
-  REAL zx_tmp((nbp_lon+1)*(nbp_lat-1/(nbp_lat-1)))
-  REAL airedyn(nbp_lon+1,nbp_lat-1/(nbp_lat-1))
-  REAL padyn(nbp_lon+1,nbp_lat-1/(nbp_lat-1),klev+1)
-  REAL dudyn(nbp_lon+1,nbp_lat-1/(nbp_lat-1),klev)
-  REAL rlatdyn(nbp_lon+1,nbp_lat-1/(nbp_lat-1))
+!  REAL zx_tmp((nbp_lon+1)*(nbp_lat-1/(nbp_lat-1)))
   !IM 190504 END
   LOGICAL ok_msk
@@ -1264,5 +1260,5 @@
      ok_sync_omp=.false.
      CALL getin('ok_sync',ok_sync_omp)
-     call phys_output_open(rlon,rlat,nCFMIP,tabijGCM, &
+     call phys_output_open(lon_degrees,lat_degrees,nCFMIP,tabijGCM, &
           iGCM,jGCM,lonGCM,latGCM, &
           jjmp1,nlevSTD,clevSTD,rlevSTD, dtime,ok_veget, &
@@ -1318,7 +1314,7 @@
              rg, &
              ra, &
-             airephy, &
-             rlat, &
-             rlon, &
+             cell_area, &
+             lat_degrees, &
+             lon_degrees, &
              presnivs, &
              calday, &
@@ -1526,5 +1522,5 @@
   IF (ip_ebil_phy.ge.1) THEN 
      ztit='after dynamic'
-     CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
+     CALL diagetpq(cell_area,ztit,ip_ebil_phy,1,1,dtime &
           , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
           , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
@@ -1533,5 +1529,5 @@
      !     est egale a la variation de la physique au pas de temps precedent.
      !     Donc la somme de ces 2 variations devrait etre nulle.
-     call diagphy(airephy,ztit,ip_ebil_phy &
+     call diagphy(cell_area,ztit,ip_ebil_phy &
           , zero_v, zero_v, zero_v, zero_v, zero_v &
           , zero_v, zero_v, zero_v, ztsol &
@@ -1675,5 +1671,5 @@
         zzz=real(90) ! could be revisited
      ENDIF
-     wo(:,:,1)=ozonecm(rlat, paprs,read_climoz,rjour=zzz)
+     wo(:,:,1)=ozonecm(lat_degrees, paprs,read_climoz,rjour=zzz)
   ENDIF
   !
@@ -1730,8 +1726,8 @@
   IF (ip_ebil_phy.ge.2) THEN 
      ztit='after reevap'
-     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime &
+     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,1,dtime &
           , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
           , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-     call diagphy(airephy,ztit,ip_ebil_phy &
+     call diagphy(cell_area,ztit,ip_ebil_phy &
           , zero_v, zero_v, zero_v, zero_v, zero_v &
           , zero_v, zero_v, zero_v, ztsol &
@@ -1777,12 +1773,12 @@
   ! non nul aux poles.
   IF (abs(solarlong0-1000.)<1.e-4) then
-     call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract)
+     call zenang_an(cycle_diurne,jH_cur,lat_degrees,lon_degrees,rmu0,fract)
   ELSE
      !  Avec ou sans cycle diurne
      IF (cycle_diurne) THEN
         zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
-        CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
+        CALL zenang(zlongi,jH_cur,zdtime,lat_degrees,lon_degrees,rmu0,fract)
      ELSE
-        CALL angle(zlongi, rlat, fract, rmu0)
+        CALL angle(zlongi, lat_degrees, fract, rmu0)
      ENDIF
   ENDIF
@@ -1854,5 +1850,5 @@
           dtime,     date0,     itap,    days_elapsed+1, &
           debut,     lafin, &
-          rlon,      rlat,      rugoro,  rmu0,      &
+          lon_degrees,      lat_degrees,      rugoro,  rmu0,      &
           zsig,      sollwdown, pphi,    cldt,      &
           rain_fall, snow_fall, solsw,   sollw,     &
@@ -1974,8 +1970,8 @@
      IF (ip_ebil_phy.ge.2) THEN 
         ztit='after surface_main'
-        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+        CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
              , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
              , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-        call diagphy(airephy,ztit,ip_ebil_phy &
+        call diagphy(cell_area,ztit,ip_ebil_phy &
              , zero_v, zero_v, zero_v, zero_v, sens &
              , evap  , zero_v, zero_v, ztsol &
@@ -2024,5 +2020,5 @@
   ENDDO
   IF (check) THEN
-     za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
+     za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
      WRITE(lunout,*) "avantcon=", za
   ENDIF
@@ -2044,5 +2040,5 @@
   DO k = 1, klev
      DO i = 1, klon
-        omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
+        omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
      END DO
   END DO
@@ -2377,8 +2373,8 @@
   IF (ip_ebil_phy.ge.2) THEN 
      ztit='after convect'
-     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
           , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
           , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-     call diagphy(airephy,ztit,ip_ebil_phy &
+     call diagphy(cell_area,ztit,ip_ebil_phy &
           , zero_v, zero_v, zero_v, zero_v, zero_v &
           , zero_v, rain_con, snow_con, ztsol &
@@ -2388,12 +2384,12 @@
   !
   IF (check) THEN
-     za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
+     za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
      WRITE(lunout,*)"aprescon=", za
      zx_t = 0.0
      za = 0.0
      DO i = 1, klon
-        za = za + airephy(i)/REAL(klon)
+        za = za + cell_area(i)/REAL(klon)
         zx_t = zx_t + (rain_con(i)+ &
-             snow_con(i))*airephy(i)/REAL(klon)
+             snow_con(i))*cell_area(i)/REAL(klon)
      ENDDO
      zx_t = zx_t/za*dtime
@@ -2542,8 +2538,8 @@
   IF (ip_ebil_phy.ge.2) THEN 
      ztit='after wake'
-     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
           , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
           , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-     call diagphy(airephy,ztit,ip_ebil_phy &
+     call diagphy(cell_area,ztit,ip_ebil_phy &
           , zero_v, zero_v, zero_v, zero_v, zero_v &
           , zero_v, zero_v, zero_v, ztsol &
@@ -2635,5 +2631,5 @@
              ,ztv,zpspsk,ztla,zthl &
              !cc nrlmd le 10/04/2012
-             ,pbl_tke_input,pctsrf,omega,airephy &
+             ,pbl_tke_input,pctsrf,omega,cell_area &
              ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
              ,n2,s2,ale_bl_stat &
@@ -2858,8 +2854,8 @@
   IF (ip_ebil_phy.ge.2) THEN 
      ztit='after dry_adjust'
-     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
           , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
           , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-     call diagphy(airephy,ztit,ip_ebil_phy &
+     call diagphy(cell_area,ztit,ip_ebil_phy &
           , zero_v, zero_v, zero_v, zero_v, zero_v &
           , zero_v, zero_v, zero_v, ztsol &
@@ -2911,12 +2907,12 @@
   ENDDO
   IF (check) THEN
-     za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
+     za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
      WRITE(lunout,*)"apresilp=", za
      zx_t = 0.0
      za = 0.0
      DO i = 1, klon
-        za = za + airephy(i)/REAL(klon)
+        za = za + cell_area(i)/REAL(klon)
         zx_t = zx_t + (rain_lsc(i) &
-             + snow_lsc(i))*airephy(i)/REAL(klon)
+             + snow_lsc(i))*cell_area(i)/REAL(klon)
      ENDDO
      zx_t = zx_t/za*dtime
@@ -2926,8 +2922,8 @@
   IF (ip_ebil_phy.ge.2) THEN 
      ztit='after fisrt'
-     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
           , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
           , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-     call diagphy(airephy,ztit,ip_ebil_phy &
+     call diagphy(cell_area,ztit,ip_ebil_phy &
           , zero_v, zero_v, zero_v, zero_v, zero_v &
           , zero_v, rain_lsc, snow_lsc, ztsol &
@@ -3170,8 +3166,8 @@
   IF (ip_ebil_phy.ge.2) THEN 
      ztit="after diagcld"
-     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
           , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
           , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-     call diagphy(airephy,ztit,ip_ebil_phy &
+     call diagphy(cell_area,ztit,ip_ebil_phy &
           , zero_v, zero_v, zero_v, zero_v, zero_v &
           , zero_v, zero_v, zero_v, ztsol &
@@ -3239,5 +3235,5 @@
         CALL AEROSOL_METEO_CALC( &
              calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
-             prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
+             prfl,psfl,pctsrf,cell_area,lat_degrees,lon_degrees,u10m,v10m)
      END IF
 
@@ -3248,7 +3244,7 @@
           jH_cur, &
           pctsrf(1,1), &
-          rlat, &
-          rlon, &
-          airephy, &
+          lat_degrees, &
+          lon_degrees, &
+          cell_area, &
           paprs, &
           pplay, &
@@ -3360,6 +3356,6 @@
         DO i=1,klon
            !
-           if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND. &
-                rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
+           if (lon_degrees(i).ge.lon1_beta.AND.lat_degrees(i).le.lon2_beta.AND. &
+                lat_degrees(i).le.lat1_beta.AND.lat_degrees(i).ge.lat2_beta) THEN
               if (pplay(i,k).GE.pfree) THEN
                  beta(i,k) = beta_pbl
@@ -3582,8 +3578,8 @@
   IF (ip_ebil_phy.ge.2) THEN 
      ztit='after rad'
-     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
           , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
           , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-     call diagphy(airephy,ztit,ip_ebil_phy &
+     call diagphy(cell_area,ztit,ip_ebil_phy &
           , topsw, toplw, solsw, sollw, zero_v &
           , zero_v, zero_v, zero_v, ztsol &
@@ -3681,5 +3677,5 @@
 
         CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
-             rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
+             lat_degrees,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
              igwd,idx,itest, &
              t_seri, u_seri, v_seri, &
@@ -3689,5 +3685,5 @@
      ELSE
         CALL lift_noro(klon,klev,dtime,paprs,pplay, &
-             rlat,zmea,zstd,zpic, &
+             lat_degrees,zmea,zstd,zpic, &
              itest, &
              t_seri, u_seri, v_seri, &
@@ -3707,5 +3703,5 @@
 
      CALL hines_gwd(klon,klev,dtime,paprs,pplay, &
-          rlat,t_seri,u_seri,v_seri, &
+          lat_degrees,t_seri,u_seri,v_seri, &
           zustrhi,zvstrhi, &
           d_t_hin, d_u_hin, d_v_hin)
@@ -3751,5 +3747,5 @@
      CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
           ra,rg,romega, &
-          rlat,rlon,pphis, &
+          lat_degrees,lon_degrees,pphis, &
           zustrdr,zustrli,zustrph, &
           zvstrdr,zvstrli,zvstrph, &
@@ -3761,8 +3757,8 @@
   IF (ip_ebil_phy.ge.2) THEN 
      ztit='after orography'
-     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
           , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
           , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-     call diagphy(airephy,ztit,ip_ebil_phy &
+     call diagphy(cell_area,ztit,ip_ebil_phy &
           , zero_v, zero_v, zero_v, zero_v, zero_v &
           , zero_v, zero_v, zero_v, ztsol &
@@ -3798,5 +3794,5 @@
              ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
              ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, &
-             klon,klev,rlon,rlat,presnivs,overlap, &
+             klon,klev,lon_degrees,lat_degrees,presnivs,overlap, &
              fract,ref_liq,ref_ice, &
              pctsrf(:,is_ter)+pctsrf(:,is_lic), &
@@ -3842,5 +3838,5 @@
        zustar,   zu10m,     zv10m, &
        wstar(:,is_ave),    ale_bl,         ale_wake, &
-       rlat,     rlon, &
+       lat_degrees,     lon_degrees, &
        frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
        presnivs, pphis,     pphi,     albsol1, &
@@ -3863,10 +3859,10 @@
           print*,'Attention on met a 0 les thermiques pour phystoke'
      call phystokenc ( &
-          nlon,klev,pdtphys,rlon,rlat, &
+          nlon,klev,pdtphys,lon_degrees,lat_degrees, &
           t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
           fm_therm,entr_therm, &
           cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
           frac_impa, frac_nucl, &
-          pphis,airephy,dtime,itap, &
+          pphis,cell_area,dtime,itap, &
           qx(:,:,ivap),da,phi,mp,upwd,dnwd)
 
@@ -3908,5 +3904,5 @@
   IF (ip_ebil_phy.ge.1) THEN 
      ztit='after physic'
-     CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
+     CALL diagetpq(cell_area,ztit,ip_ebil_phy,1,1,dtime &
           , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
           , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
@@ -3916,5 +3912,5 @@
      !     Donc la somme de ces 2 variations devrait etre nulle.
 
-     call diagphy(airephy,ztit,ip_ebil_phy &
+     call diagphy(cell_area,ztit,ip_ebil_phy &
           , topsw, toplw, solsw, sollw, sens &
           , evap, rain_fall, snow_fall, ztsol &
@@ -3964,5 +3960,5 @@
           paprs, &
           q_seri, &
-          airephy, &
+          cell_area, &
           pphi, &
           pphis, &
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/print_debug_phys.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/print_debug_phys.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/print_debug_phys.F90	(revision 3825)
@@ -4,4 +4,5 @@
 use phys_local_var_mod
 use phys_state_var_mod
+use geometry_mod, ONLY : lon_degrees, lat_degrees
 IMPLICIT NONE
 integer i,debug_lev
@@ -14,5 +15,5 @@
 print*,'l    u, v, T, q, ql'
 DO k = 1, klev
-   write(*,'(i3,2f8.4,3f14.4,2e14.2)') k,rlon(i),rlat(i),u_seri(i,k),v_seri(i,k),t_seri(i,k),q_seri(i,k),ql_seri(i,k)
+   write(*,'(i3,2f8.4,3f14.4,2e14.2)') k,lon_degrees(i),lat_degrees(i),u_seri(i,k),v_seri(i,k),t_seri(i,k),q_seri(i,k),ql_seri(i,k)
 ENDDO
 
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/radlwsw_m.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/radlwsw_m.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/radlwsw_m.F90	(revision 3825)
@@ -769,5 +769,5 @@
       ENDDO
 !----- Fin des mises a zero des tableaux output de RECMWF -------------------              
-!        GEMU(1:klon)=sin(rlatd(1:klon))
+!        GEMU(1:klon)=sin(latitude(1:klon))
 ! On met les donnees dans l'ordre des niveaux arpege
          paprs_i(:,1)=paprs(:,klev+1)
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/regr_lat_time_climoz_m.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/regr_lat_time_climoz_m.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/regr_lat_time_climoz_m.F90	(revision 3825)
@@ -73,7 +73,5 @@
     use netcdf, only: nf90_get_att, nf90_get_var, nf90_noerr, nf90_nowrite
     use assert_m, only: assert
-    !use comconst_phy_mod, only : pi
-    !use comgeom2_phy_mod, only : rlatv
-    use comgeomphy, only: rlatv
+    use regular_lonlat_mod, only : boundslat_reg, south
     use nrtype, only: pi
     integer, intent(in):: read_climoz ! read ozone climatology
@@ -293,5 +291,5 @@
             "Found 12 months in ozone climatologies, assuming periodicity..."
        o3_regr_lat(nbp_lat:1:-1, :, 1:12, :) = regr1_step_av(o3_in, &
-            xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(nbp_lat-1:1:-1), pi / 2/)))
+            xs=sin(lat_in_edg), xt=sin((/- pi / 2, boundslat_reg(nbp_lat-1:1:-1,south), pi / 2/)))
        ! (invert order of indices in "o3_regr_lat" because "rlatu" is
        ! in descending order)
@@ -304,5 +302,5 @@
        print *, "Using 14 months in ozone climatologies..."
        o3_regr_lat(nbp_lat:1:-1, :, :, :) = regr1_step_av(o3_in, &
-            xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(nbp_lat-1:1:-1), pi / 2/)))
+            xs=sin(lat_in_edg), xt=sin((/- pi / 2, boundslat_reg(nbp_lat-1:1:-1,south), pi / 2/)))
        ! (invert order of indices in "o3_regr_lat" because "rlatu" is
        ! in descending order)
@@ -334,8 +332,6 @@
          nf95_put_att, nf95_enddef, nf95_copy_att, nf95_put_var
     use netcdf, only: nf90_clobber, nf90_float, nf90_global
-    !use comconst_phy_mod, only : pi
-    !use comgeom2_phy_mod, only : rlatu
     use nrtype, only: pi
-    use comgeomphy, only: rlatu
+    use regular_lonlat_mod, only : lat_reg
     integer, intent(in):: ncid_in, n_plev
     integer, intent(out):: ncid_out, varid_plev, varid_time
@@ -421,5 +417,5 @@
 
     ! Write one of the coordinate variables:
-    call nf95_put_var(ncid_out, varid_rlatu, rlatu(nbp_lat:1:-1) / pi * 180.)
+    call nf95_put_var(ncid_out, varid_rlatu, lat_reg(nbp_lat:1:-1) / pi * 180.)
     ! (convert from rad to degrees and sort in ascending order)
 
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/regr_lat_time_coefoz_m.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/regr_lat_time_coefoz_m.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/regr_lat_time_coefoz_m.F90	(revision 3825)
@@ -49,5 +49,5 @@
     !use comconst_phy_mod, only : pi
     use nrtype, only: pi
-    use comgeomphy, only: rlatv
+    use regular_lonlat_mod, only: boundslat_reg, south
     ! Variables local to the procedure:
 
@@ -211,5 +211,5 @@
        ! equivalent to weighting by cosine of latitude:
        v_regr_lat(nbp_lat:1:-1, :, 1:12) = regr1_step_av(o3_par_in, &
-            xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(nbp_lat-1:1:-1), pi / 2/)))
+            xs=sin(lat_in_edg), xt=sin((/- pi / 2, boundslat_reg(nbp_lat-1:1:-1,south), pi / 2/)))
        ! (invert order of indices in "v_regr_lat" because "rlatu" is
        ! in descending order)
@@ -251,6 +251,5 @@
     !use comconst_phy_mod, only : pi
     use nrtype, only: pi
-    use comgeomphy, only: rlatu
-
+    use regular_lonlat_mod, only : lat_reg
     integer, intent(in):: ncid_in, varid_in(:), n_plev
     character(len=*), intent(in):: name_out(:) ! of NetCDF variables
@@ -326,5 +325,5 @@
 
     ! Write one of the coordinate variables:
-    call nf95_put_var(ncid_out, varid_rlatu, rlatu(nbp_lat:1:-1) / pi * 180.)
+    call nf95_put_var(ncid_out, varid_rlatu, lat_reg(nbp_lat:1:-1) / pi * 180.)
     ! (convert from rad to degrees and sort in ascending order)
 
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/regular_lonlat_mod.f90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/regular_lonlat_mod.f90	(revision 3825)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/regular_lonlat_mod.f90	(revision 3825)
@@ -0,0 +1,53 @@
+MODULE regular_lonlat_mod
+
+  
+  INTEGER, PARAMETER :: north_east=1     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: north_west=2     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: south_west=3     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: south_east=4     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: east=1           ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: west=2           ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: north=1          ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: south=2          ! boundaries of regular lontlat
+
+  REAL,SAVE,ALLOCATABLE :: lon_reg(:)      ! value of longitude cell (rad)
+!$OMP THREADPRIVATE(lon_reg)
+
+  REAL,SAVE,ALLOCATABLE :: lat_reg(:)      ! value of longitude cell (rad)
+!$OMP THREADPRIVATE(lat_reg)
+
+  REAL,SAVE,ALLOCATABLE :: boundslon_reg(:,:)      ! value of boundaries cell (1=>east, 2=>west)(rad)
+!$OMP THREADPRIVATE(boundslon_reg)
+
+  REAL,SAVE,ALLOCATABLE :: boundslat_reg(:,:)      ! value of longitude cell (1=>north, 2=>south)(rad)
+!$OMP THREADPRIVATE(boundslat_reg)
+
+
+CONTAINS
+  
+  SUBROUTINE init_regular_lonlat(nbp_lon, nbp_lat, lon_reg_, lat_reg_, boundslon_reg_, boundslat_reg_)
+    IMPLICIT NONE
+    INTEGER,INTENT(IN) :: nbp_lon
+    INTEGER,INTENT(IN) :: nbp_lat
+    REAL,   INTENT(IN) :: lon_reg_(nbp_lon)  
+    REAL,   INTENT(IN) :: lat_reg_(nbp_lat) 
+    REAL,   INTENT(IN) :: boundslon_reg_(nbp_lon,2) 
+    REAL,   INTENT(IN) :: boundslat_reg_(nbp_lat,2)  
+
+    
+    ALLOCATE(lon_reg(nbp_lon))
+    lon_reg(:)=lon_reg_(:)
+
+    ALLOCATE(lat_reg(nbp_lat)) 
+    lat_reg(:)=lat_reg_(:)
+
+    ALLOCATE(boundslon_reg(nbp_lon,2)) 
+    boundslon_reg(:,:)=boundslon_reg_(:,:)
+
+    ALLOCATE(boundslat_reg(nbp_lat,2))  
+    boundslat_reg(:,:)=boundslat_reg_(:,:)
+
+  
+  END SUBROUTINE init_regular_lonlat    
+
+END MODULE regular_lonlat_mod
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/surf_land_bucket_mod.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/surf_land_bucket_mod.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/surf_land_bucket_mod.F90	(revision 3825)
@@ -24,5 +24,5 @@
     USE cpl_mod
     USE dimphy
-    USE comgeomphy
+    USE geometry_mod, ONLY : latitude
     USE mod_grid_phy_lmdz
     USE mod_phys_lmdz_para
@@ -99,5 +99,5 @@
     CALL calbeta(dtime, is_ter, knon, snow, qsol, beta, capsol, dif_grnd)
     if (type_veget=='betaclim') then
-       CALL calbeta_clim(knon,jour,rlatd(knindex(1:knon)),beta)
+       CALL calbeta_clim(knon,jour,latitude(knindex(1:knon)),beta)
     endif
        
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/surf_land_orchidee_mod.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/surf_land_orchidee_mod.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/surf_land_orchidee_mod.F90	(revision 3825)
@@ -16,5 +16,5 @@
   USE cpl_mod,      ONLY : cpl_send_land_fields
   USE surface_data, ONLY : type_ocean
-  USE comgeomphy,   ONLY : cuphy, cvphy
+  USE geometry_mod, ONLY : dx, dy
   USE mod_grid_phy_lmdz
   USE mod_phys_lmdz_para
@@ -328,6 +328,6 @@
        DO igrid = 1, knon
           ij = knindex(igrid)
-          resolution(igrid,1) = cuphy(ij)
-          resolution(igrid,2) = cvphy(ij)
+          resolution(igrid,1) = dx(ij)
+          resolution(igrid,2) = dy(ij)
        ENDDO
      
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/thermcellV0_main.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/thermcellV0_main.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/thermcellV0_main.F90	(revision 3825)
@@ -13,5 +13,4 @@
 
       USE dimphy
-      USE comgeomphy , ONLY:rlond,rlatd
       USE inifis_mod, ONLY: lunout, prt_level
       IMPLICIT NONE
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/thermcell_main.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/thermcell_main.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/thermcell_main.F90	(revision 3825)
@@ -23,5 +23,4 @@
       USE dimphy
       USE ioipsl
-      USE comgeomphy , ONLY:rlond,rlatd
       USE indice_sol_mod
       USE inifis_mod, ONLY: lunout, prt_level
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/tracinca_mod.F90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/tracinca_mod.F90	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/tracinca_mod.F90	(revision 3825)
@@ -44,5 +44,5 @@
     USE infotrac_phy
     USE vampir
-    USE comgeomphy
+    USE geometry_mod, ONLY : cell_area
     !USE control_phy_mod
     USE inifis_mod, ONLY: config_inca
@@ -144,5 +144,5 @@
          pplay,pdel,prfl,pmflxr,psfl,    &
          pmflxs,zpmfu,itop_con,ibas_con,  &
-         pphi,airephy,nstep,rneb,t_seri, &      
+         pphi,cell_area,nstep,rneb,t_seri, &      
          rh,tau_aero,piz_aero,cg_aero,   &
          rfname,ccm,lafin, config_inca)
@@ -161,5 +161,5 @@
          pplay,      & !pmid
          pdel,       & !pdel
-         airephy,    &
+         cell_area,    &
          pctsrf(1,1),& !oro
          ftsol,      & !tsurf
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/vertical_layers_mod.f90
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/vertical_layers_mod.f90	(revision 3825)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/vertical_layers_mod.f90	(revision 3825)
@@ -0,0 +1,36 @@
+MODULE vertical_layers_mod
+
+   REAL,SAVE             :: preff  ! reference surface pressure (Pa)
+   REAL,SAVE,ALLOCATABLE :: ap(:)
+   REAL,SAVE,ALLOCATABLE :: bp(:)
+   REAL,SAVE,ALLOCATABLE :: presnivs(:)
+   REAL,SAVE,ALLOCATABLE :: pseudoalt(:)
+!$OMP THREADPRIVATE(preff,ap,bp,presnivs,pseudoalt)
+
+
+CONTAINS
+
+  SUBROUTINE init_vertical_layers(nlayer,preff_,ap_,bp_,&
+                                 presnivs_, pseudoalt_)
+    IMPLICIT NONE
+    INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
+    REAL,INTENT(IN)    :: preff_ ! reference surface pressure (Pa)
+    REAL,INTENT(IN)    :: ap_(nlayer+1) ! hybrid coordinate at interfaces
+    REAL,INTENT(IN)    :: bp_(nlayer+1) ! hybrid coordinate at interfaces
+    REAL,INTENT(IN)    :: presnivs_(nlayer) ! Appproximative pressure of atm. layers (Pa)
+    REAL,INTENT(IN)    :: pseudoalt_(nlayer) ! pseudo-altitude of atm. layers (km)
+  
+    ALLOCATE(ap(nlayer+1))
+    ALLOCATE(bp(nlayer+1))
+    ALLOCATE(presnivs(nlayer))
+    ALLOCATE(pseudoalt(nlayer))
+  
+    preff = preff_
+    ap(:) = ap_(:)
+    bp(:) = bp_(:)
+    presnivs(:) = presnivs_(:)
+    pseudoalt(:) = pseudoalt_(:)
+
+  END SUBROUTINE init_vertical_layers
+
+END MODULE vertical_layers_mod
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/write_histday_seri.h
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/write_histday_seri.h	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/write_histday_seri.h	(revision 3825)
@@ -20,5 +20,5 @@
       ok_msk=.FALSE.
       msk(1:klon)=pctsrf(1:klon,is_ter)
-      CALL moyglo_pondaire(klon, zx_tmp_fi2d, airephy,  &
+      CALL moyglo_pondaire(klon, zx_tmp_fi2d, cell_area,  &
            ok_msk, msk, moyglo)
       zx_tmp_fi2d(1:klon)=moyglo
@@ -29,5 +29,5 @@
 !
       ok_msk=.FALSE.
-      CALL moyglo_pondaire(klon, bils, airephy,  &
+      CALL moyglo_pondaire(klon, bils, cell_area,  &
            ok_msk, msk, moyglo)
       zx_tmp_fi2d(1:klon)=moyglo
@@ -45,5 +45,5 @@
 !
       CALL moyglo_pondaima(klon, klev, zx_tmp_fi3d,  &
-           airephy, paprs, moyglo)
+           cell_area, paprs, moyglo)
       zx_tmp_fi2d(1:klon)=moyglo
 !
@@ -65,5 +65,5 @@
 !ym!
 !ym      CALL moyglo_pondaima(klon, klev, zx_tmp_fi3d,  &
-!ym           airephy, paprs, moyglo)
+!ym           cell_area, paprs, moyglo)
 !ym      zx_tmp_fi2d(1:klon)=moyglo
 !ym!
@@ -79,5 +79,5 @@
 !ym!
 !ym      ok_msk=.FALSE.
-!ym      CALL moyglo_pondaire(klon, zx_tmp_fi2d, airephy,  &
+!ym      CALL moyglo_pondaire(klon, zx_tmp_fi2d, cell_area,  &
 !ym           ok_msk, msk, moyglo)
 !ym      zx_tmp_fi2d(1:klon)=moyglo
@@ -90,5 +90,5 @@
 !ym!
 !ym!IM 190504 BEG
-!ym      CALL gr_fi_dyn(1,klon,nbp_lon+1,jjm+1,airephy,airedyn)
+!ym      CALL gr_fi_dyn(1,klon,nbp_lon+1,jjm+1,cell_area,airedyn)
 !ym      CALL gr_fi_dyn(klev+1,klon,nbp_lon+1,jjm+1,paprs,padyn)
 !ym      CALL gr_fi_dyn(1,klon,nbp_lon+1,jjm+1,rlat,rlatdyn)
@@ -126,8 +126,8 @@
 !ym!
 !ym
-!ym      CALL gr_fi_dyn(1,klon,nbp_lon+1,jjm+1,airephy,airedyn)
+!ym      CALL gr_fi_dyn(1,klon,nbp_lon+1,jjm+1,cell_area,airedyn)
 !#endif      
       
-      CALL grid1dTo2d_glo(airephy,zx_tmp_2d)
+      CALL grid1dTo2d_glo(cell_area,zx_tmp_2d)
       airetot=0.
 !     DO j = 1, jjmp1
@@ -145,5 +145,5 @@
       airetot=0.
       DO i=1, klon
-       airetot=airetot+airephy(i)
+       airetot=airetot+cell_area(i)
       ENDDO
 !     IF(itap.EQ.1) PRINT*,'airetotphy=',airetot
@@ -172,5 +172,5 @@
 !
       CALL moyglo_pondmass(klon, klev, t_seri, &
-           airephy, paprs, moyglo)
+           cell_area, paprs, moyglo)
       zx_tmp_fi2d(1:klon)=moyglo
 !
@@ -180,5 +180,5 @@
 !
       ok_msk=.FALSE.
-      CALL moyglo_pondaire(klon, paprs(:,1), airephy,  &
+      CALL moyglo_pondaire(klon, paprs(:,1), cell_area,  &
            ok_msk, msk, moyglo)
       zx_tmp_fi2d(1:klon)=moyglo
@@ -189,5 +189,5 @@
 !
       ok_msk=.FALSE.
-      CALL moyglo_pondaire(klon, evap, airephy,  &
+      CALL moyglo_pondaire(klon, evap, cell_area,  &
            ok_msk, msk, moyglo)
       zx_tmp_fi2d(1:klon)=moyglo
@@ -203,5 +203,5 @@
 !     ok_msk=.TRUE.
 !     msk(1:klon)=pctsrf(1:klon,is_ter)
-!     CALL moyglo_pondaire(klon, zx_tmp_fi2d, airephy, 
+!     CALL moyglo_pondaire(klon, zx_tmp_fi2d, cell_area, 
 !    .                     ok_msk, msk, moyglo)
 !     zx_tmp_fi2d(1:klon)=moyglo
@@ -219,5 +219,5 @@
 !     ok_msk=.TRUE.
 !     msk(1:klon)=pctsrf(1:klon,is_ter)
-!     CALL moyglo_pondaire(klon, zx_tmp_fi2d, airephy, 
+!     CALL moyglo_pondaire(klon, zx_tmp_fi2d, cell_area, 
 !    .     ok_msk, msk, moyglo)
 !     zx_tmp_fi2d(1:klon)=moyglo
@@ -233,5 +233,5 @@
       ok_msk=.TRUE.
       msk(1:klon)=pctsrf(1:klon,is_oce)
-      CALL moyglo_pondaire(klon, zx_tmp_fi2d, airephy,  &
+      CALL moyglo_pondaire(klon, zx_tmp_fi2d, cell_area,  &
            ok_msk, msk, moyglo)
       zx_tmp_fi2d(1:klon)=moyglo
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/write_histrac.h
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/write_histrac.h	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/write_histrac.h	(revision 3825)
@@ -8,5 +8,5 @@
      
      CALL histwrite_phy(nid_tra,.FALSE.,"phis",itau_w,pphis)
-     CALL histwrite_phy(nid_tra,.FALSE.,"aire",itau_w,airephy)
+     CALL histwrite_phy(nid_tra,.FALSE.,"aire",itau_w,cell_area)
      CALL histwrite_phy(nid_tra,.FALSE.,"zmasse",itau_w,zmasse)
 ! RomP >>>
Index: /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/write_paramLMDZ_phy.h
===================================================================
--- /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/write_paramLMDZ_phy.h	(revision 3824)
+++ /dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/write_paramLMDZ_phy.h	(revision 3825)
@@ -2,24 +2,24 @@
 ! calcul moyennes globales
 !
-       zx_tmp_fi2d=bils*airephy
-       CALL global_mean(zx_tmp_fi2d,airephy,.TRUE.,gbils)
-       zx_tmp_fi2d=evap*airephy
-       CALL global_mean(zx_tmp_fi2d,airephy,.TRUE.,gevap)
-       zx_tmp_fi2d(:)=fevap(:, is_ter)*airephy(:)
-       CALL global_mean(zx_tmp_fi2d,airephy,.TRUE.,gevapt)
-       zx_tmp_fi2d=zxfluxlat*airephy
-       CALL global_mean(zx_tmp_fi2d,airephy,.TRUE.,glat)
-       zx_tmp_fi2d=(topsw0-toplw0)*airephy
-       CALL global_mean(zx_tmp_fi2d,airephy,.TRUE.,gnet0)
-       zx_tmp_fi2d=(topsw-toplw)*airephy
-       CALL global_mean(zx_tmp_fi2d,airephy,.TRUE.,gnet)
-       zx_tmp_fi2d=(rain_fall+snow_fall)*airephy
-       CALL global_mean(zx_tmp_fi2d,airephy,.TRUE.,grain)
-       zx_tmp_fi2d=zxtsol*airephy
-       CALL global_mean(zx_tmp_fi2d,airephy,.TRUE.,gtsol)
-       zx_tmp_fi2d=zt2m*airephy
-       CALL global_mean(zx_tmp_fi2d,airephy,.TRUE.,gt2m)
-       zx_tmp_fi2d=prw*airephy
-       CALL global_mean(zx_tmp_fi2d,airephy,.TRUE.,gprw)
+       zx_tmp_fi2d=bils*cell_area
+       CALL global_mean(zx_tmp_fi2d,cell_area,.TRUE.,gbils)
+       zx_tmp_fi2d=evap*cell_area
+       CALL global_mean(zx_tmp_fi2d,cell_area,.TRUE.,gevap)
+       zx_tmp_fi2d(:)=fevap(:, is_ter)*cell_area(:)
+       CALL global_mean(zx_tmp_fi2d,cell_area,.TRUE.,gevapt)
+       zx_tmp_fi2d=zxfluxlat*cell_area
+       CALL global_mean(zx_tmp_fi2d,cell_area,.TRUE.,glat)
+       zx_tmp_fi2d=(topsw0-toplw0)*cell_area
+       CALL global_mean(zx_tmp_fi2d,cell_area,.TRUE.,gnet0)
+       zx_tmp_fi2d=(topsw-toplw)*cell_area
+       CALL global_mean(zx_tmp_fi2d,cell_area,.TRUE.,gnet)
+       zx_tmp_fi2d=(rain_fall+snow_fall)*cell_area
+       CALL global_mean(zx_tmp_fi2d,cell_area,.TRUE.,grain)
+       zx_tmp_fi2d=zxtsol*cell_area
+       CALL global_mean(zx_tmp_fi2d,cell_area,.TRUE.,gtsol)
+       zx_tmp_fi2d=zt2m*cell_area
+       CALL global_mean(zx_tmp_fi2d,cell_area,.TRUE.,gt2m)
+       zx_tmp_fi2d=prw*cell_area
+       CALL global_mean(zx_tmp_fi2d,cell_area,.TRUE.,gprw)
 !
 !$OMP MASTER
