Index: codes/icosagcm/trunk/src/physics.f90
===================================================================
--- codes/icosagcm/trunk/src/physics.f90	(revision 212)
+++ codes/icosagcm/trunk/src/physics.f90	(revision 213)
@@ -7,6 +7,7 @@
   INTEGER, PARAMETER :: phys_none=0, phys_HS94=1, phys_DCMIP=2
 
-  INTEGER :: nb_extra_physics_2D, nb_extra_physics_3D, phys_type
+  INTEGER :: phys_type
   TYPE(t_field),POINTER :: f_extra_physics_2D(:), f_extra_physics_3D(:)
+  TYPE(t_field),POINTER :: f_dulon(:), f_dulat(:)
 
   CHARACTER(LEN=255) :: physics_type="automatic"
@@ -14,4 +15,5 @@
 
   PUBLIC :: physics, init_physics
+
 CONTAINS
 
@@ -21,7 +23,12 @@
     USE icosa
     USE physics_interface_mod
-    USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics
-    IMPLICIT NONE
-
+    USE physics_dcmip_mod, ONLY : &
+         init_physics_dcmip=>init_physics, init_physics_dcmip_new=>init_physics_new
+    IMPLICIT NONE
+
+    CALL allocate_field(f_dulon,field_t,type_real,llm, name='dulon')
+    CALL allocate_field(f_dulat,field_t,type_real,llm, name='dulat')
+    CALL init_pack_before ! Compute physics_inout%ngrid and offsets used by pack/unpack
+    
     physics_type='automatic'
     CALL getin("physics",physics_type)
@@ -40,7 +47,6 @@
 
     CASE ('dcmip')
-       CALL init_physics_dcmip
-       nb_extra_physics_2D=1 ! precl
-       nb_extra_physics_3D=0
+       CALL init_physics_dcmip_new
+!       CALL init_physics_dcmip
        phys_type = phys_DCMIP
     CASE DEFAULT
@@ -55,7 +61,11 @@
        PRINT *, 'nb_extra_physics_3D = ', nb_extra_physics_3D
     END IF
-    IF(nb_extra_physics_2D>0) CALL allocate_field(f_extra_physics_2D,field_t,type_real,nb_extra_physics_2D)
-    IF(nb_extra_physics_3D>0) CALL allocate_field(f_extra_physics_3D,field_t,type_real,llm,nb_extra_physics_3D)
-
+
+    IF(.FALSE.) THEN ! draft interface
+       IF(nb_extra_physics_2D>0) CALL allocate_field(f_extra_physics_2D,field_t,type_real,nb_extra_physics_2D)
+       IF(nb_extra_physics_3D>0) CALL allocate_field(f_extra_physics_3D,field_t,type_real,llm,nb_extra_physics_3D)
+    ELSE
+       CALL init_pack_after ! Defines Ai, lon, lat in physics_inout
+    END IF
   END SUBROUTINE init_physics
 
@@ -63,5 +73,5 @@
     USE icosa
     USE physics_interface_mod
- !   USE etat0_mod
+    USE physics_dcmip_mod, ONLY : write_physics_dcmip => write_physics
     USE etat0_heldsz_mod
     IMPLICIT NONE
@@ -81,43 +91,123 @@
     LOGICAL:: firstcall,lastcall
     INTEGER :: ind
-    TYPE(physics_inout) :: args
+    TYPE(t_physics_inout) :: args
+
+    IF(MOD(it+1,itau_physics)==0) THEN
+    
+       SELECT CASE(phys_type)
+       CASE (phys_none)
+          ! No physics, do nothing
+       CASE(phys_HS94)
+          CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 
+       CASE DEFAULT
+          IF(.FALSE.) THEN ! draft interface
+             args%dt_phys = dt*itau_physics
+             DO ind=1,ndomain
+                IF (.NOT. assigned_domain(ind)) CYCLE
+                CALL swap_dimensions(ind)
+                CALL swap_geometry(ind)
+                
+                phis=f_phis(ind)
+                ps=f_ps(ind)
+                theta_rhodz=f_theta_rhodz(ind)
+                ue=f_ue(ind)
+                q=f_q(ind)
+                
+                IF(nb_extra_physics_2D>0) args%extra_2D=f_extra_physics_2D(ind)
+                IF(nb_extra_physics_3D>0) args%extra_3D=f_extra_physics_3D(ind)
+                CALL physics_column(args, phis, ps, theta_rhodz, ue, q)
+             ENDDO
+             
+             IF (mod(it,itau_out)==0 ) THEN
+                IF(nb_extra_physics_2D>0) CALL writefield("extra_physics_2D",f_extra_physics_2D)
+                IF(nb_extra_physics_3D>0) CALL writefield("extra_physics_3D",f_extra_physics_3D)
+             ENDIF
+          ELSE ! new interface
+             physics_inout%dt_phys = dt*itau_physics
+             CALL physics_column_new(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
+          END IF
+       END SELECT
+
+       CALL transfert_request(f_theta_rhodz,req_i0)
+       CALL transfert_request(f_ue,req_e0_vect)
+       CALL transfert_request(f_q,req_i0)
+    END IF
+
+    IF (mod(it,itau_out)==0 ) THEN
+       SELECT CASE(phys_type)
+       CASE (phys_DCMIP)
+          CALL write_physics_dcmip
+       END SELECT
+    END IF
+    
+  END SUBROUTINE physics
+
+!--------------------------------- New interface --------------------------------------
+
+  SUBROUTINE physics_column_new(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
+    USE icosa
+    USE physics_interface_mod
+    USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics
+    IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_phis(:)
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_ue(:)
+    TYPE(t_field),POINTER :: f_q(:)
+    REAL(rstd),POINTER :: phis(:)
+    REAL(rstd),POINTER :: ps(:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:)
+    REAL(rstd),POINTER :: ue(:,:)
+    REAL(rstd),POINTER :: dulon(:,:)
+    REAL(rstd),POINTER :: dulat(:,:)
+    REAL(rstd),POINTER :: q(:,:,:)
+    INTEGER :: it, ind
+
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       phis=f_phis(ind)
+       ps=f_ps(ind)
+       theta_rhodz=f_theta_rhodz(ind)
+       ue=f_ue(ind)
+       q=f_q(ind)
+       CALL pack_physics(pack_info(ind), phis, ps, theta_rhodz, ue, q)
+    END DO
 
     SELECT CASE(phys_type)
-    CASE (phys_none)
-       ! No physics, do nothing
-    CASE(phys_HS94)
-       CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 
-    CASE DEFAULT
-       CALL transfert_request(f_ps,req_i1)
-       CALL transfert_request(f_theta_rhodz,req_i1)
-       CALL transfert_request(f_ue,req_e1_vect)
-       CALL transfert_request(f_q,req_i1)
-       
-       args%dt_phys = dt
-       DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind)) CYCLE
-          CALL swap_dimensions(ind)
-          CALL swap_geometry(ind)
-      
-          phis=f_phis(ind)
-          ps=f_ps(ind)
-          theta_rhodz=f_theta_rhodz(ind)
-          ue=f_ue(ind)
-          q=f_q(ind)
-
-          IF(nb_extra_physics_2D>0) args%extra_2D=f_extra_physics_2D(ind)
-          IF(nb_extra_physics_3D>0) args%extra_3D=f_extra_physics_3D(ind)
-          CALL physics_column(args, phis, ps, theta_rhodz, ue, q)
-       ENDDO
-
-       IF (mod(it,itau_out)==0 ) THEN
-          IF(nb_extra_physics_2D>0) CALL writefield("extra_physics_2D",f_extra_physics_2D)
-          IF(nb_extra_physics_3D>0) CALL writefield("extra_physics_3D",f_extra_physics_3D)
-       ENDIF
+    CASE (phys_DCMIP)
+       CALL full_physics_dcmip
     END SELECT
 
-  END SUBROUTINE physics
-
-  SUBROUTINE physics_column(args, phis, ps, theta_rhodz, ue, q)
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       ps=f_ps(ind)
+       theta_rhodz=f_theta_rhodz(ind)
+       q=f_q(ind)
+       dulon=f_dulon(ind)
+       dulat=f_dulat(ind)
+       CALL unpack_physics(pack_info(ind), ps, theta_rhodz, q, dulon, dulat)
+    END DO
+
+    ! Transfer dulon, dulat
+    CALL transfert_request(f_dulon,req_i0)
+    CALL transfert_request(f_dulat,req_i0)
+
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       ue=f_ue(ind)
+       dulon=f_dulon(ind)
+       dulat=f_dulat(ind)
+       CALL compute_update_velocity(dulon, dulat, ue)
+    END DO
+
+  END SUBROUTINE physics_column_new
+
+  SUBROUTINE pack_physics(info, phis, ps, theta_rhodz, ue, q)
     USE icosa
     USE wind_mod
@@ -125,7 +215,98 @@
     USE theta2theta_rhodz_mod
     USE physics_interface_mod
+    IMPLICIT NONE
+    TYPE(t_pack_info) :: info
+    REAL(rstd) :: phis(iim*jjm)
+    REAL(rstd) :: ps(iim*jjm)
+    REAL(rstd) :: theta_rhodz(iim*jjm,llm)
+    REAL(rstd) :: ue(3*iim*jjm,llm)
+    REAL(rstd) :: q(iim*jjm,llm,nqtot)
+
+    REAL(rstd) :: p(iim*jjm,llm+1)
+    REAL(rstd) :: Temp(iim*jjm,llm)
+    REAL(rstd) :: uc(iim*jjm,3,llm)
+    REAL(rstd) :: ulon(iim*jjm,llm)
+    REAL(rstd) :: ulat(iim*jjm,llm)
+
+    CALL compute_pression(ps,p,0)
+    CALL compute_theta_rhodz2temperature(ps,theta_rhodz,Temp,0)
+    CALL compute_wind_centered(ue,uc)
+    CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
+
+    CALL pack_domain(info, phis, physics_inout%phis)
+    CALL pack_domain(info, p, physics_inout%p)
+    CALL pack_domain(info, Temp, physics_inout%Temp)
+    CALL pack_domain(info, ulon, physics_inout%ulon)
+    CALL pack_domain(info, ulat, physics_inout%ulat)
+    CALL pack_domain(info, q, physics_inout%q)
+  END SUBROUTINE pack_physics
+
+  SUBROUTINE unpack_physics(info, ps,theta_rhodz, q, dulon, dulat)
+    USE icosa
+    USE physics_interface_mod
+    USE theta2theta_rhodz_mod
+    IMPLICIT NONE
+    TYPE(t_pack_info) :: info
+    REAL(rstd) :: ps(iim*jjm)
+    REAL(rstd) :: theta_rhodz(iim*jjm,llm)
+    REAL(rstd) :: Temp(iim*jjm,llm)
+    REAL(rstd) :: q(iim*jjm,llm,nqtot)
+    REAL(rstd) :: dulon(iim*jjm,llm)
+    REAL(rstd) :: dulat(iim*jjm,llm)
+
+    REAL(rstd) :: dq(iim*jjm,llm,nqtot)
+    REAL(rstd) :: dTemp(iim*jjm,llm)
+    CALL unpack_domain(info, dulon, physics_inout%dulon)
+    CALL unpack_domain(info, dulat, physics_inout%dulat)
+    CALL unpack_domain(info, dq, physics_inout%dq)
+    CALL unpack_domain(info, Temp, physics_inout%Temp)
+    CALL unpack_domain(info, dTemp, physics_inout%dTemp)
+    q = q + physics_inout%dt_phys * dq
+    Temp = Temp + physics_inout%dt_phys * dTemp
+    CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
+  END SUBROUTINE unpack_physics
+
+  SUBROUTINE compute_update_velocity(dulon, dulat, ue)
+    USE icosa
+    USE physics_interface_mod
+    USE wind_mod
+    IMPLICIT NONE
+    REAL(rstd) :: dulon(iim*jjm,llm)
+    REAL(rstd) :: dulat(iim*jjm,llm)
+    REAL(rstd) :: ue(3*iim*jjm,llm)
+    REAL(rstd) :: duc(iim*jjm,3,llm)
+    REAL(rstd) :: dt2, due
+    INTEGER :: i,j,ij,l
+    ! Reconstruct wind tendencies at edges and add to normal wind
+    CALL compute_wind_centered_from_lonlat_compound(dulon,dulat,duc)
+    dt2=.5*physics_inout%dt_phys
+    DO l=1,llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          due = sum( (duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
+          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due
+
+          due = sum( (duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
+          ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due
+
+          due = sum( (duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
+          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due
+        ENDDO
+      ENDDO
+    ENDDO
+  END SUBROUTINE compute_update_velocity
+
+!--------------------------------- Draft interface --------------------------------------
+
+  SUBROUTINE physics_column(args, phis, ps, theta_rhodz, ue, q)
+    USE icosa
+    USE wind_mod
+    USE pression_mod
+    USE theta2theta_rhodz_mod
+    USE physics_interface_mod
     USE physics_dcmip_mod
     IMPLICIT NONE
-    TYPE(physics_inout) :: args    
+    TYPE(t_physics_inout) :: args    
     REAL(rstd) :: phis(iim*jjm)
     REAL(rstd) :: ps(iim*jjm)
@@ -190,5 +371,5 @@
         DO i=ii_begin,ii_end
           ij=(j-1)*iim+i
-          due = sum((uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
+          due = sum( (uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
           ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due
 
Index: codes/icosagcm/trunk/src/physics_dcmip.f90
===================================================================
--- codes/icosagcm/trunk/src/physics_dcmip.f90	(revision 212)
+++ codes/icosagcm/trunk/src/physics_dcmip.f90	(revision 213)
@@ -7,24 +7,81 @@
 
   TYPE(t_field),POINTER :: f_out_i(:)
-  TYPE(t_field),POINTER :: f_precl(:)
   REAL(rstd),POINTER :: out_i(:,:)
 
-  PUBLIC :: compute_phys_wrap, init_physics
+  TYPE(t_field),POINTER  :: f_precl(:)
+  REAL(rstd),ALLOCATABLE :: precl_packed(:)
+
+  PUBLIC :: compute_phys_wrap, init_physics, &
+       init_physics_new, full_physics, write_physics
 
 CONTAINS
 
+!-------------------------- New interface ----------------------
+
+  SUBROUTINE init_physics_new
+    USE physics_interface_mod
+    IMPLICIT NONE
+    INTEGER :: ngrid
+    testcase=1 ! OK for 4.2 (moist baroclinic instability)
+    CALL getin("dcmip_physics",testcase)
+    ngrid = physics_inout%ngrid
+    ! Input
+    ALLOCATE(physics_inout%Ai(ngrid))
+    ALLOCATE(physics_inout%lon(ngrid))
+    ALLOCATE(physics_inout%lat(ngrid))
+    ALLOCATE(physics_inout%phis(ngrid))
+    ALLOCATE(physics_inout%p(ngrid,llm+1))
+    ALLOCATE(physics_inout%Temp(ngrid,llm))
+    ALLOCATE(physics_inout%ulon(ngrid,llm))
+    ALLOCATE(physics_inout%ulat(ngrid,llm))
+    ALLOCATE(physics_inout%q(ngrid,llm,nqtot))
+    ! Output (tendencies)
+    ALLOCATE(physics_inout%dTemp(ngrid,llm))
+    ALLOCATE(physics_inout%dulon(ngrid,llm))
+    ALLOCATE(physics_inout%dulat(ngrid,llm))
+    ALLOCATE(physics_inout%dq(ngrid,llm,nqtot))
+    ! Physics-specific data
+    ALLOCATE(precl_packed(ngrid))
+    CALL allocate_field(f_precl, field_t,type_real)
+
+    PRINT *, 'init_physics_new', SIZE(physics_inout%Ai)
+  END SUBROUTINE init_physics_new
+
+  SUBROUTINE full_physics
+    USE physics_interface_mod
+    CALL compute_physics(physics_inout%ngrid, physics_inout%dt_phys, &
+         physics_inout%lat, physics_inout%p, physics_inout%Temp, & 
+         physics_inout%ulon, physics_inout%ulat, physics_inout%q(:,:,1), &
+         physics_inout%dTemp, physics_inout%dulon, physics_inout%dulat, &
+         physics_inout%dq(:,:,1), precl_packed)
+  END SUBROUTINE full_physics
+
+  SUBROUTINE write_physics
+    USE output_field_mod
+    USE physics_interface_mod
+    CALL unpack_field(f_precl, precl_packed)
+    CALL output_field("precl",f_precl)
+  END SUBROUTINE write_physics
+
+!------------------------ Draft interface -----------------------
+
   SUBROUTINE init_physics
+    USE physics_interface_mod
     IMPLICIT NONE
     testcase=1 ! OK for 4.2 (moist baroclinic instability)
     CALL getin("dcmip_physics",testcase)
+    nb_extra_physics_2D=1 ! precl
+    nb_extra_physics_3D=0
   END SUBROUTINE init_physics
 
   SUBROUTINE compute_phys_wrap(args)
     USE physics_interface_mod
-    TYPE(physics_inout) :: args
+    TYPE(t_physics_inout) :: args
     CALL compute_physics(args%ngrid, args%dt_phys, args%lat, &
          args%p, args%Temp, args%ulon, args%ulat, args%q(:,:,1), &
          args%dTemp, args%dulon, args%dulat, args%dq(:,:,1), args%extra_2D(:,1))
   END SUBROUTINE compute_phys_wrap
+
+!------------------ Interface-independent wrapper ---------------------------
 
   SUBROUTINE compute_physics(ngrid,dt_phys,lat, p,Temp,u,v,q, dTemp,du,dv,dq, precl)
@@ -54,5 +111,5 @@
     REAL(rstd) :: qfi(ngrid,llm)
 
-    INTEGER :: i,j,l,ll,ij
+    INTEGER :: l,ll,ij
     REAL(rstd) :: dt_phys, inv_dt
 
@@ -61,9 +118,6 @@
 
     DO l=1,llm+1
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
+      DO ij=1,ngrid
           pint(ij,l)=p(ij,llm+2-l)
-        ENDDO
       ENDDO
     ENDDO
@@ -71,14 +125,11 @@
     DO l=1,llm
        ll=llm+1-l
-       DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-             pmid(ij,l)=0.5*(pint(ij,l)+pint(ij,l+1)) ! Pressure inside layers
-             pdel(ij,l)=pint(ij,l+1)-pint(ij,l)       ! Pressure difference between two layers
-             ufi(ij,l)=u(ij,ll)
-             vfi(ij,l)=v(ij,ll)
-             qfi(ij,l)=q(ij,ll)
-             Tfi(ij,l)=Temp(ij,ll)/(1+0.608*qfi(ij,l))
-          ENDDO
+       DO ij=1,ngrid
+          pmid(ij,l)=0.5*(pint(ij,l)+pint(ij,l+1)) ! Pressure inside layers
+          pdel(ij,l)=pint(ij,l+1)-pint(ij,l)       ! Pressure difference between two layers
+          ufi(ij,l)=u(ij,ll)
+          vfi(ij,l)=v(ij,ll)
+          qfi(ij,l)=q(ij,ll)
+          Tfi(ij,l)=Temp(ij,ll)/(1+0.608*qfi(ij,l))  
        ENDDO
     ENDDO
@@ -90,12 +141,9 @@
     DO l=1,llm
        ll=llm+1-l
-       DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-             dTemp(ij,l) = inv_dt * ( Tfi(ij,ll)*(1+0.608*qfi(ij,ll)) - Temp(ij,l) )
-             du(ij,l) = inv_dt * (ufi(ij,ll) - u(ij,l))
-             dv(ij,l) = inv_dt * (vfi(ij,ll) - v(ij,l))
-             dq(ij,l)  = inv_dt * (qfi(ij,ll) - q(ij,l))
-          ENDDO
+       DO ij=1,ngrid
+          dTemp(ij,l) = inv_dt * ( Tfi(ij,ll)*(1+0.608*qfi(ij,ll)) - Temp(ij,l) )
+          du(ij,l) = inv_dt * (ufi(ij,ll) - u(ij,l))
+          dv(ij,l) = inv_dt * (vfi(ij,ll) - v(ij,l))
+          dq(ij,l)  = inv_dt * (qfi(ij,ll) - q(ij,l))
        ENDDO
     ENDDO
Index: codes/icosagcm/trunk/src/physics_interface.f90
===================================================================
--- codes/icosagcm/trunk/src/physics_interface.f90	(revision 212)
+++ codes/icosagcm/trunk/src/physics_interface.f90	(revision 213)
@@ -5,9 +5,12 @@
   PRIVATE
 
-  TYPE :: physics_inout
-     ! Input
+  INTEGER :: nb_extra_physics_2D, nb_extra_physics_3D
+
+  TYPE t_physics_inout
+     ! Input, time-independent
      INTEGER :: ngrid
      REAL(rstd) :: dt_phys
-     REAL(rstd), DIMENSION(:), POINTER :: lon, lat
+     REAL(rstd), DIMENSION(:), POINTER :: Ai, lon, lat, phis
+     ! Input, time-dependent
      REAL(rstd), DIMENSION(:,:), POINTER :: p, Temp, ulon, ulat
      REAL(rstd), DIMENSION(:,:,:), POINTER :: q
@@ -18,7 +21,385 @@
      REAL(rstd), DIMENSION(:,:), POINTER :: extra_2D
      REAL(rstd), DIMENSION(:,:,:), POINTER :: extra_3D
-  END TYPE physics_inout
-
-  PUBLIC :: physics_inout
+  END TYPE t_physics_inout
+
+!------------------------ (new interface) --------------------------
+! physics_inout is used to exchange information with physics
+! Field ngrid is initialized by physics.f90/init_physics. Its other fields
+! must be defined by XX/init_physics (where XX =  e.g. physics_dcmip.f90) 
+! by either pointing to internal data of the physics package
+! or by a specific allocation
+! size : (ngrid), (ngrid,llm) except p(ngrid,llm+1), (ngrid,llm,nqtot) 
+
+  TYPE(t_physics_inout), SAVE :: physics_inout
+ 
+!------------------------ (new interface) --------------------------
+! pack_info contains indices used by pack/unpack routines
+! to pack together the data of all the domains managed by the MPI process
+! It is initialized by physics.f90/init_physics 
+  TYPE t_pack_info
+     INTEGER :: ngrid, & ! number of non-halo points in that domain
+          nseg ! number of segments (contigous parts) in that domain
+     ! size and start of each segment : ij domain index, k packed index
+     INTEGER, ALLOCATABLE :: n(:), ij(:), k(:)
+  END TYPE t_pack_info
+
+  TYPE(t_pack_info), ALLOCATABLE, SAVE :: pack_info(:)
+  
+  INTERFACE pack_field
+     MODULE PROCEDURE pack_2D
+     MODULE PROCEDURE pack_3D
+     MODULE PROCEDURE pack_4D
+  END INTERFACE pack_field
+
+  INTERFACE unpack_field
+     MODULE PROCEDURE unpack_2D
+     MODULE PROCEDURE unpack_3D
+     MODULE PROCEDURE unpack_4D
+  END INTERFACE unpack_field
+
+  INTERFACE pack_domain
+     MODULE PROCEDURE pack_domain_2D
+     MODULE PROCEDURE pack_domain_3D
+     MODULE PROCEDURE pack_domain_4D
+  END INTERFACE pack_domain
+
+  INTERFACE unpack_domain
+     MODULE PROCEDURE unpack_domain_2D
+     MODULE PROCEDURE unpack_domain_3D
+     MODULE PROCEDURE unpack_domain_4D
+  END INTERFACE unpack_domain
+
+  PUBLIC :: nb_extra_physics_2D, nb_extra_physics_3D, &
+       t_physics_inout, physics_inout, &
+       t_pack_info, pack_info, init_pack_before, init_pack_after, &
+       pack_domain, pack_field, unpack_domain, unpack_field, &
+       garbage_3D
+
+CONTAINS
+  
+    SUBROUTINE init_pack_before
+    USE icosa
+    IMPLICIT NONE
+    INTEGER :: ind, offset
+    !$OMP MASTER
+    offset=0
+    ALLOCATE(pack_info(ndomain))
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       CALL count_segments(domain(ind)%own, pack_info(ind))
+       pack_info(ind)%k = pack_info(ind)%k + offset
+       offset = offset + pack_info(ind)%ngrid
+    END DO
+    physics_inout%ngrid = offset
+    !$OMP END MASTER
+    !$OMP BARRIER
+  END SUBROUTINE init_pack_before
+
+  SUBROUTINE count_segments(own, info)
+    USE icosa
+    IMPLICIT NONE
+    LOGICAL, DIMENSION(:,:) :: own
+    TYPE(t_pack_info) :: info
+
+    INTEGER, DIMENSION(jjm) :: n
+    INTEGER :: ngrid, nseg, i, j, jj, k
+    INTEGER, PARAMETER :: method=4
+    SELECT CASE(method)
+    CASE(1) ! Copy all points, including halo (works)
+       info%nseg=1
+       info%ngrid=iim*jjm
+       ALLOCATE(info%n(1))
+       ALLOCATE(info%ij(1))
+       ALLOCATE(info%k(1))
+       info%n(1)=iim*jjm
+       info%ij(1)=1
+       info%k(1)=1
+    CASE(2) ! Copy all points, including halo, one at a time (works, slow ?)
+       info%nseg=iim*jjm
+       info%ngrid=iim*jjm
+       ALLOCATE(info%n(iim*jjm))
+       ALLOCATE(info%ij(iim*jjm))
+       ALLOCATE(info%k(iim*jjm))
+       DO jj=1,iim*jjm
+          info%n(jj) =1
+          info%ij(jj)=jj
+          info%k(jj) =jj
+       END DO
+    CASE(3) ! Copy non-halo points only, one at a time (works, slow ?)
+       n=COUNT(own,1)
+       ngrid=SUM(n)
+       info%ngrid=ngrid
+       info%nseg=ngrid
+       ALLOCATE(info%n(ngrid))
+       ALLOCATE(info%ij(ngrid))
+       ALLOCATE(info%k(ngrid))
+       jj=1
+       DO j=1,jjm
+          DO i=1,iim
+             IF(own(i,j)) THEN
+                info%n(jj)=1
+                info%k(jj)=jj
+                info%ij(jj) = iim*(j-1)+i
+                jj=jj+1
+             END IF
+          END DO
+       END DO
+       
+    CASE DEFAULT ! Copy non-halo points only, as contiguous segments (works)
+       n=0
+       n=COUNT(own,1)
+       ngrid=SUM(n)
+       info%ngrid=ngrid
+       nseg=COUNT(n>0)
+       info%nseg=nseg
+       ALLOCATE(info%n(nseg))
+       ALLOCATE(info%ij(nseg))
+       ALLOCATE(info%k(nseg))
+
+       jj=1
+       k=1
+       DO j=1,jjm
+          IF(n(j)>0) THEN
+             ! find first .TRUE. value in own(:,j)
+             DO i=1,iim
+                IF(own(i,j)) THEN
+                   info%n(jj)=n(j)
+                   info%k(jj)=k
+                   info%ij(jj) = iim*(j-1)+i
+                   IF(COUNT(own(i:i+n(j)-1,j)) /= n(j)) STOP
+                   EXIT
+                END IF
+             END DO
+             k = k + n(j)
+             jj=jj+1
+          END IF
+       END DO
+
+       IF(k-1/=ngrid) THEN
+          PRINT *, 'Total number of grid points inconsistent', k-1, ngrid
+          STOP
+       END IF
+       IF(jj-1/=nseg) THEN
+          PRINT *, 'Number of segments inconsistent', jj-1, nseg
+          STOP
+       END IF
+
+    END SELECT
+    
+    PRINT *, 'count_segments', info%nseg, info%ngrid, SUM(info%n), COUNT(own), iim*jjm
+  END SUBROUTINE count_segments
+
+  SUBROUTINE init_pack_after
+    USE icosa
+    IMPLICIT NONE
+    INTEGER :: ind, offset
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       CALL pack_domain_2D(pack_info(ind), Ai, physics_inout%Ai)
+       CALL pack_lonlat(pack_info(ind))
+    END DO
+  END SUBROUTINE init_pack_after
+
+  SUBROUTINE pack_lonlat(info)
+    USE icosa
+    IMPLICIT NONE
+    TYPE(t_pack_info) :: info
+    REAL(rstd) :: lon(iim*jjm), lat(iim*jjm)
+    INTEGER :: ij
+    DO ij=1,iim*jjm
+       CALL xyz2lonlat(xyz_i(ij,:),lon(ij),lat(ij)) 
+    ENDDO
+    CALL pack_domain_2D(info, lon, physics_inout%lon)
+    CALL pack_domain_2D(info, lat, physics_inout%lat)
+  END SUBROUTINE pack_lonlat
+
+!-------------------------------- Pack / Unpack 2D ---------------------------
+
+  SUBROUTINE pack_2D(f_2D, packed)
+    USE icosa
+    IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_2D(:)
+    REAL(rstd)            :: packed(:)
+    REAL(rstd), POINTER   :: loc(:)
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       loc = f_2D(ind)
+       CALL pack_domain_2D(pack_info(ind), loc, packed)
+    END DO
+  END SUBROUTINE pack_2D
+
+  SUBROUTINE unpack_2D(f_2D, packed) 
+    USE icosa
+    IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_2D(:)
+    REAL(rstd)            :: packed(:)
+    REAL(rstd), POINTER   :: loc(:)
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       loc = f_2D(ind)
+       CALL unpack_domain_2D(pack_info(ind), loc, packed)
+    END DO
+  END SUBROUTINE unpack_2D
+
+  SUBROUTINE pack_domain_2D(info, loc, glob)
+    USE icosa
+    IMPLICIT NONE
+    TYPE(t_pack_info) :: info
+    REAL(rstd), DIMENSION(:) :: loc, glob
+    INTEGER :: jj,n,k,ij
+    DO jj=1, info%nseg
+       n = info%n(jj)-1
+       ij = info%ij(jj)
+       k = info%k(jj)
+       glob(k:k+n) = loc(ij:ij+n)
+    END DO
+  END SUBROUTINE pack_domain_2D
+
+  SUBROUTINE unpack_domain_2D(info, loc, glob)
+    IMPLICIT NONE
+    TYPE(t_pack_info) :: info
+    REAL(rstd), DIMENSION(:) :: loc, glob
+    INTEGER :: jj,n,k,ij
+    DO jj=1, info%nseg
+       n = info%n(jj)-1
+       ij = info%ij(jj)
+       k = info%k(jj)
+       loc(ij:ij+n) = glob(k:k+n)
+    END DO
+  END SUBROUTINE unpack_domain_2D
+
+!-------------------------------- Pack / Unpack 3D ---------------------------
+
+  SUBROUTINE pack_3D(f_3D, packed)
+    USE icosa
+    IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_3D(:)
+    REAL(rstd)            :: packed(:,:)
+    REAL(rstd), POINTER   :: loc(:,:)
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       loc = f_3D(ind)
+       CALL pack_domain_3D(pack_info(ind), loc, packed)
+    END DO
+  END SUBROUTINE pack_3D
+
+  SUBROUTINE unpack_3D(f_3D, packed) 
+    USE icosa
+    IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_3D(:)
+    REAL(rstd)            :: packed(:,:)
+    REAL(rstd), POINTER   :: loc(:,:)
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       loc = f_3D(ind)
+       CALL unpack_domain_3D(pack_info(ind), loc, packed)
+    END DO
+  END SUBROUTINE unpack_3D
+
+  SUBROUTINE pack_domain_3D(info, loc, glob)
+    IMPLICIT NONE
+    TYPE(t_pack_info) :: info
+    REAL(rstd), DIMENSION(:,:) :: loc, glob
+    INTEGER :: jj,n,k,ij
+    DO jj=1, info%nseg
+       n = info%n(jj)-1
+       ij = info%ij(jj)
+       k = info%k(jj)
+       glob(k:k+n,:) = loc(ij:ij+n,:)
+    END DO
+  END SUBROUTINE pack_domain_3D
+
+  SUBROUTINE unpack_domain_3D(info, loc, glob)
+    IMPLICIT NONE
+    TYPE(t_pack_info) :: info
+    REAL(rstd), DIMENSION(:,:) :: loc, glob
+    INTEGER :: jj,n,k,ij
+    DO jj=1, info%nseg
+       n = info%n(jj)-1
+       ij = info%ij(jj)
+       k = info%k(jj)
+       loc(ij:ij+n,:) = glob(k:k+n,:)
+    END DO    
+  END SUBROUTINE unpack_domain_3D
+
+  SUBROUTINE garbage_3D(loc,own)
+    USE icosa
+    IMPLICIT NONE
+    LOGICAL :: own(iim,jjm)
+    REAL(rstd) :: loc(iim*jjm,llm)
+    INTEGER :: i,j,ij
+    ! write garbage in non-owned points
+    DO j=1,jjm
+       DO i=1,iim
+          IF(.NOT.own(i,j)) THEN
+             ij=iim*(j-1)+i
+             loc(ij,:)=-1e30
+          END IF
+       END DO
+    END DO   
+  END SUBROUTINE garbage_3D
+
+!-------------------------------- Pack / Unpack 4D ---------------------------
+
+  SUBROUTINE pack_4D(f_4D, packed)
+    USE icosa
+    IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_4D(:)
+    REAL(rstd)            :: packed(:,:,:)
+    REAL(rstd), POINTER   :: loc(:,:,:)
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       loc = f_4D(ind)
+       CALL pack_domain_4D(pack_info(ind), loc, packed)
+    END DO
+  END SUBROUTINE pack_4D
+
+  SUBROUTINE unpack_4D(f_4D, packed) 
+    USE icosa
+    IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_4D(:)
+    REAL(rstd)            :: packed(:,:,:)
+    REAL(rstd), POINTER   :: loc(:,:,:)
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       loc = f_4D(ind)
+       CALL unpack_domain_4D(pack_info(ind), loc, packed)
+    END DO
+  END SUBROUTINE unpack_4D
+
+  SUBROUTINE pack_domain_4D(info, loc, glob)
+    IMPLICIT NONE
+    TYPE(t_pack_info) :: info
+    REAL(rstd), DIMENSION(:,:,:) :: loc, glob
+    INTEGER :: jj,n,k,ij
+    DO jj=1, info%nseg
+       n = info%n(jj)-1
+       ij = info%ij(jj)
+       k = info%k(jj)
+       glob(k:k+n,:,:) = loc(ij:ij+n,:,:)
+    END DO
+  END SUBROUTINE pack_domain_4D
+
+  SUBROUTINE unpack_domain_4D(info, loc, glob)
+    IMPLICIT NONE
+    TYPE(t_pack_info) :: info
+    REAL(rstd), DIMENSION(:,:,:) :: loc, glob
+    INTEGER :: jj,n,k,ij
+    DO jj=1, info%nseg
+       n = info%n(jj)-1
+       ij = info%ij(jj)
+       k = info%k(jj)
+       loc(ij:ij+n,:,:) = glob(k:k+n,:,:)
+    END DO
+  END SUBROUTINE unpack_domain_4D
 
 END MODULE physics_interface_mod
