Index: trunk/ICOSA_LMDZ/src/phymars/interface_icosa_lmdz.f90
===================================================================
--- trunk/ICOSA_LMDZ/src/phymars/interface_icosa_lmdz.f90	(revision 2271)
+++ trunk/ICOSA_LMDZ/src/phymars/interface_icosa_lmdz.f90	(revision 2271)
@@ -0,0 +1,926 @@
+MODULE interface_icosa_lmdz_mod
+
+  USE field_mod, ONLY: t_field
+  USE transfert_mod, ONLY: t_message 
+  
+ 
+  TYPE(t_message),SAVE :: req_u
+  TYPE(t_message),SAVE :: req_dps0, req_dulon0, req_dulat0, req_dTemp0, req_dq0
+
+  TYPE(t_field),POINTER,SAVE :: f_p(:) 
+  TYPE(t_field),POINTER,SAVE :: f_pks(:)  
+  TYPE(t_field),POINTER,SAVE :: f_pk(:)  
+  TYPE(t_field),POINTER,SAVE :: f_p_layer(:)    
+  TYPE(t_field),POINTER,SAVE :: f_theta(:)   
+  TYPE(t_field),POINTER,SAVE :: f_phi(:)   
+  TYPE(t_field),POINTER,SAVE :: f_Temp(:)    
+  TYPE(t_field),POINTER,SAVE :: f_ulon(:)    
+  TYPE(t_field),POINTER,SAVE :: f_ulat(:)   
+  TYPE(t_field),POINTER,SAVE :: f_dulon(:)
+  TYPE(t_field),POINTER,SAVE :: f_dulat(:)
+  TYPE(t_field),POINTER,SAVE :: f_dTemp(:)
+  TYPE(t_field),POINTER,SAVE :: f_dq(:)
+  TYPE(t_field),POINTER,SAVE :: f_dps(:)
+  TYPE(t_field),POINTER,SAVE :: f_duc(:)
+  TYPE(t_field),POINTER,SAVE :: f_bounds_lon(:)
+  TYPE(t_field),POINTER,SAVE :: f_bounds_lat(:)
+
+  INTEGER,SAVE :: start_clock
+  INTEGER,SAVE :: stop_clock
+  INTEGER,SAVE :: count_clock=0
+  
+  REAL,SAVE :: day_length ! length of a day (s)
+  REAL,SAVE :: year_length ! length of a year (s)
+!!!  INTEGER,SAVE :: start_day ! reference sol value at beginning of the run wrt Ls=0
+  INTEGER,SAVE :: day_ini ! reference sol value at beginning of the run wrt Ls=0
+  REAL,SAVE :: hour_ini ! reference sol value at beginning of the run wrt Ls=0
+  
+  INTEGER,SAVE :: nbp_phys
+  INTEGER,SAVE :: nbp_phys_glo
+  
+  CHARACTER(len=30),SAVE,ALLOCATABLE :: tname(:) ! tracer names
+  REAL,SAVE :: pday ! number of ellapsed sols since Ls=0
+  REAL,SAVE :: ptime ! "universal time" as fraction of sol (e.g. 0.5 for noon)
+
+
+CONTAINS
+
+  SUBROUTINE initialize_physics
+  USE distrib_icosa_lmdz_mod, ONLY : init_distrib_icosa_lmdz, transfer_icosa_to_lmdz
+! from dynamico
+  USE domain_mod
+  USE dimensions
+  USE mpi_mod
+  USE mpipara
+  USE disvert_mod
+  USE xios_mod
+  USE time_mod , init_time_icosa=> init_time 
+  USE transfert_mod
+  
+! from LMDZ
+  USE mod_grid_phy_lmdz, ONLY : unstructured
+  USE mod_phys_lmdz_omp_data, ONLY: klon_omp
+  USE transfert_mod
+  USE physics_distribution_mod, ONLY : init_physics_distribution
+   
+  
+  IMPLICIT NONE
+  INTEGER  :: ind,i,j,ij,pos
+  REAL(rstd),POINTER :: bounds_lon(:,:)
+  REAL(rstd),POINTER :: bounds_lat(:,:)
+  
+  REAL(rstd),ALLOCATABLE :: latfi(:)
+  REAL(rstd),ALLOCATABLE :: lonfi(:)
+  REAL(rstd),ALLOCATABLE :: airefi(:)
+  REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:)
+  REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:)
+  REAL(rstd) :: pseudoalt(llm)
+
+  INTEGER :: nbp_phys, nbp_phys_glo
+  
+!$OMP PARALLEL
+    CALL allocate_field(f_bounds_lon,field_t,type_real,6)
+    CALL allocate_field(f_bounds_lat,field_t,type_real,6)
+    CALL allocate_field(f_p,field_t,type_real,llm+1,name="p_in")
+    CALL allocate_field(f_pks,field_t,type_real)
+    CALL allocate_field(f_pk,field_t,type_real,llm)
+    CALL allocate_field(f_p_layer,field_t,type_real,llm,name="p_layer_in")
+    CALL allocate_field(f_theta,field_t,type_real,llm)
+    CALL allocate_field(f_phi,field_t,type_real,llm,name="phi_in")
+    CALL allocate_field(f_Temp,field_t,type_real,llm,name="Temp_in")
+    CALL allocate_field(f_ulon,field_t,type_real,llm,name="ulon_in")
+    CALL allocate_field(f_ulat,field_t,type_real,llm,name="ulat_in")
+    CALL allocate_field(f_dulon,field_t,type_real,llm,name="dulon_out")
+    CALL allocate_field(f_dulat,field_t,type_real,llm,name="dulat_out")
+    CALL allocate_field(f_dTemp,field_t,type_real,llm,name="dTemp_out")
+    CALL allocate_field(f_dq,field_t,type_real,llm,nqtot,name="dq_out")
+    CALL allocate_field(f_dps,field_t,type_real,name="dps_out")
+    CALL allocate_field(f_duc,field_t,type_real,3,llm)    
+
+    CALL init_message(f_dps,req_i0,req_dps0)
+    CALL init_message(f_dulon,req_i0,req_dulon0)
+    CALL init_message(f_dulat,req_i0,req_dulat0)
+    CALL init_message(f_dTemp,req_i0,req_dTemp0)
+    CALL init_message(f_dq,req_i0,req_dq0)
+!$OMP END PARALLEL    
+
+    nbp_phys=0
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          IF (domain(ind)%own(i,j)) nbp_phys=nbp_phys+1
+        ENDDO
+      ENDDO
+    ENDDO
+    
+
+!initialize LMDZ5 physic mpi decomposition
+    CALL MPI_ALLREDUCE(nbp_phys,nbp_phys_glo,1,MPI_INTEGER,MPI_SUM,comm_icosa,ierr)
+    CALL init_physics_distribution(unstructured, 6, nbp_phys, 1, nbp_phys_glo, llm, comm_icosa)
+    
+    DO ind=1,ndomain
+        CALL swap_dimensions(ind)
+        CALL swap_geometry(ind)
+        bounds_lon=f_bounds_lon(ind)
+        bounds_lat=f_bounds_lat(ind)
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            CALL xyz2lonlat(xyz_v(ij+z_rup,:), bounds_lon(ij,1), bounds_lat(ij,1))
+            CALL xyz2lonlat(xyz_v(ij+z_up,:), bounds_lon(ij,2), bounds_lat(ij,2))
+            CALL xyz2lonlat(xyz_v(ij+z_lup,:), bounds_lon(ij,3), bounds_lat(ij,3))
+            CALL xyz2lonlat(xyz_v(ij+z_ldown,:), bounds_lon(ij,4), bounds_lat(ij,4))
+            CALL xyz2lonlat(xyz_v(ij+z_down,:), bounds_lon(ij,5), bounds_lat(ij,5))
+            CALL xyz2lonlat(xyz_v(ij+z_rdown,:), bounds_lon(ij,6), bounds_lat(ij,6))
+         ENDDO
+       ENDDO            
+    ENDDO
+          
+!$OMP PARALLEL
+    CALL initialize_physics_omp
+!$OMP END PARALLEL            
+
+    CALL xios_set_context    
+
+
+     
+
+  END SUBROUTINE initialize_physics
+
+
+  SUBROUTINE initialize_physics_omp
+  USE distrib_icosa_lmdz_mod, ONLY : init_distrib_icosa_lmdz, transfer_icosa_to_lmdz
+! from dynamico
+  USE domain_mod
+  USE dimensions
+  USE mpi_mod
+  USE mpipara
+  USE disvert_mod
+  USE xios_mod
+  USE time_mod , ONLY: init_time_icosa=> init_time, dt, itaumax, itau_physics
+  USE omp_para
+
+! from LMDZ
+  USE mod_grid_phy_lmdz, ONLY : unstructured
+  USE mod_phys_lmdz_omp_data, ONLY: klon_omp
+  USE time_phylmdz_mod, ONLY: init_time_lmdz => init_time
+  USE transfert_mod
+  USE physics_distribution_mod, ONLY : init_physics_distribution
+  USE dimphy, ONLY: init_dimphy
+  USE geometry_mod, ONLY : init_geometry
+  USE vertical_layers_mod, ONLY : init_vertical_layers
+  USE phys_state_var_init_mod, ONLY : phys_state_var_init
+  use comgeomfi_h, only: ini_fillgeom
+ 
+  USE netcdf  
+  
+  IMPLICIT NONE
+
+
+
+  INTEGER  :: ind,i,j,ij,pos
+  REAL(rstd),POINTER :: bounds_lon(:,:)
+  REAL(rstd),POINTER :: bounds_lat(:,:)
+  
+  REAL(rstd),ALLOCATABLE :: latfi(:)
+  REAL(rstd),ALLOCATABLE :: lonfi(:)
+  REAL(rstd),ALLOCATABLE :: airefi(:)
+  REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:)
+  REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:)
+  REAL(rstd),ALLOCATABLE :: ind_cell_glo(:)
+
+  REAL(rstd) :: pseudoalt(llm)
+  REAL(rstd) :: aps(llm)
+  REAL(rstd) :: bps(llm)
+  real(rstd) :: scaleheight
+
+  INTEGER :: run_length  
+  INTEGER :: annee_ref  
+  INTEGER :: day_ref    
+  INTEGER :: day_ini    
+  REAL    :: start_time 
+  REAL    :: physics_timestep    
+
+
+!  INTEGER                       :: nqo, nbtr
+!  CHARACTER(len=20),ALLOCATABLE :: tname(:)    ! tracer short name for restart and diagnostics
+!  CHARACTER(len=23),ALLOCATABLE :: ttext(:)     ! tracer long name for diagnostics
+!  INTEGER,ALLOCATABLE           :: niadv(:)    ! equivalent dyn / physique
+!  INTEGER,ALLOCATABLE           :: conv_flg(:) ! conv_flg(it)=0 : convection desactivated for tracer number it 
+!  INTEGER,ALLOCATABLE           :: pbl_flg(:)  ! pbl_flg(it)=0  : boundary layer diffusion desactivaded for tracer number it 
+!  CHARACTER(len=8),ALLOCATABLE  :: solsym(:)  ! tracer name from inca
+  TYPE(t_field),POINTER,SAVE    :: f_ind_cell_glo(:)
+  
+  INTEGER :: iflag_phys    
+  INTEGER :: nq
+
+  !! to get starting date
+  !! --------------------
+  logical :: startphy_file
+  ! NetCDF stuff
+  integer :: status ! NetCDF return code
+  integer :: ncid ! NetCDF file ID
+  integer :: varid ! NetCDF variable ID
+  real :: tab_cntrl(100)
+
+    CALL init_distrib_icosa_lmdz
+    
+    ALLOCATE(latfi(klon_omp))
+    ALLOCATE(lonfi(klon_omp))
+    ALLOCATE(airefi(klon_omp))
+    ALLOCATE(bounds_latfi(klon_omp,6))
+    ALLOCATE(bounds_lonfi(klon_omp,6))
+    ALLOCATE(ind_cell_glo(klon_omp))
+
+    CALL transfer_icosa_to_lmdz(geom%lat_i,latfi)
+    CALL transfer_icosa_to_lmdz(geom%lon_i,lonfi)
+    CALL transfer_icosa_to_lmdz(f_bounds_lat,bounds_latfi)
+    CALL transfer_icosa_to_lmdz(f_bounds_lon,bounds_lonfi)
+    CALL transfer_icosa_to_lmdz(geom%Ai,airefi)
+
+    CALL allocate_field(f_ind_cell_glo,field_t,type_real)
+    
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          f_ind_cell_glo(ind)%rval2d(ij)=domain(ind)%assign_cell_glo(i,j)
+        ENDDO
+      ENDDO
+    ENDDO
+
+      
+    CALL transfer_icosa_to_lmdz(f_ind_cell_glo,ind_cell_glo)
+    CALL deallocate_field(f_ind_cell_glo)
+      
+              
+    ! Initialize dimphy module
+    CALL init_dimphy(klon_omp,llm)
+
+    CALL init_geometry(klon_omp, lonfi, latfi, bounds_lonfi, bounds_latfi, airefi, INT(ind_cell_glo))
+
+    scaleheight=scale_height/1000. ! Atmospheric scale height (km)
+    aps(1:llm)=0.5*(ap(1:llm)+ap(2:llm+1))
+    bps(1:llm)=0.5*(bp(1:llm)+bp(2:llm+1))
+    pseudoalt(:)=-scaleheight*log(presnivs(:)/preff)
+    CALL init_vertical_layers(llm,preff,scaleheight,ap,bp,aps,bps,presnivs,pseudoalt)
+
+!!!    ! Initialize planet_mod (quite redundant wrt vertical_levels...)
+!!!    CALL ini_planete_mod(llm,preff,ap,bp)
+
+!$OMP MASTER
+    ALLOCATE(tname(nqtot))
+!$OMP END MASTER
+!$OMP BARRIER    
+
+! read tname() from traceur.def file
+    IF (is_mpi_root) THEN
+!$OMP MASTER
+    OPEN(unit=42,file="traceur.def",form="formatted",status="old",iostat=ierr)
+    IF (ierr==0) THEN
+      READ(42,*) nq ! should be the same as nqtot
+      IF (nq /= nqtot) THEN
+        WRITE(*,*) "Error: number of tracers in tracer.def should match nqtot!"
+        WRITE(*,*) "       will just use nqtot=",nqtot," tracers"
+      ENDIF
+      DO i=1,nqtot
+        READ(42,*) tname(i)
+      ENDDO
+      CLOSE(42)
+    ENDIF
+!$OMP END MASTER
+!$OMP BARRIER
+    ENDIF ! of (is_mpi_root)
+
+    DO i=1,nqtot
+      CALL bcast(tname(i))
+    ENDDO
+
+    startphy_file=.true.
+    CALL getin('startphy_file',startphy_file)
+
+    IF (startphy_file) THEN
+
+      IF (is_mpi_root) THEN
+!$OMP MASTER      
+      status=nf90_open('startfi.nc',NF90_NOWRITE,ncid)
+      if (status.ne.nf90_noerr) then
+        write(*,*)"Failed to open startfi.nc"
+        write(*,*)trim(nf90_strerror(status))
+        stop
+      endif
+
+      status=nf90_inq_varid(ncid,"controle",varid)
+      if (status.ne.nf90_noerr) then
+        write(*,*)"Failed to find controle variable"
+        write(*,*)trim(nf90_strerror(status))
+        stop
+      endif
+
+      status=nf90_get_var(ncid,varid,tab_cntrl)
+      day_ini=tab_cntrl(4)
+      print*,"initialize_physics_omp: day_ini",day_ini
+      hour_ini=tab_cntrl(29)
+      print*,"initialize_physics_omp: hour_ini",hour_ini
+      status=nf90_close(ncid)
+!$OMP END MASTER      
+!$OMP BARRIER
+      ENDIF ! of !IF (is_mpi_root)
+      
+      CALL bcast(day_ini)
+      CALL bcast(hour_ini)
+
+    ELSE
+
+      day_ini=0
+      hour_ini=0.
+      CALL getin('day_ini',day_ini)
+      CALL getin('hour_ini',hour_ini)
+
+    ENDIF
+
+    day_length=88775
+    CALL getin('day_length',day_length)
+!!!    year_length = 31557600 ! 365.25 * 86400
+!!!    CALL getin('year_length',year_length)
+    ndays=nint(itaumax*(dt/day_length))! number of days to run
+    physics_timestep=dt*itau_physics
+
+!!!    CALL inifis(klon_omp,llm,nqtot,start_day,day_length,ndays,physics_timestep, &
+!!!            latfi,lonfi,airefi,radius,g,kappa*cpp,cpp)
+
+   CALL phys_state_var_init(klon_omp,llm,nqtot,tname, &
+                       day_ini,hour_ini,day_length,physics_timestep, &
+                       radius,g,kappa*cpp,cpp)
+
+   CALL ini_fillgeom(klon_omp,latfi,lonfi,airefi)
+
+   CALL conf_phys(klon_omp,llm,nqtot)
+    
+  ! init time
+!    annee_ref=2015
+!    CALL getin("anneeref",annee_ref)
+    
+!    day_ref=1
+!    CALL getin("dayref",day_ref)
+    
+!    physics_timestep=dt*itau_physics
+!    run_length=itaumax*dt
+!    ndays=NINT(run_length/day_length)
+    
+!    day_ini=INT(itau0*dt/day_length)+day_ref
+!    start_time= itau0*dt/day_length-INT(itau0*dt/day_length)
+
+!    CALL init_time_lmdz(annee_ref, day_ref, day_ini, start_time, ndays, physics_timestep)
+
+!  Additional initializations for aquaplanets
+!    CALL getin("iflag_phys",iflag_phys)
+!    IF (iflag_phys>=100) THEN
+!      CALL iniaqua(klon_omp, iflag_phys)
+!    END IF
+
+    ! Initialize "calendar"
+!$OMP MASTER
+    ! initialize pday and ptime
+    pday = day_ini
+    ptime = hour_ini
+!$OMP END MASTER
+!$OMP BARRIER
+  
+  END SUBROUTINE  initialize_physics_omp 
+  
+  
+
+
+  SUBROUTINE physics
+  USE icosa
+  USE time_mod
+  USE disvert_mod
+  USE transfert_mod
+  USE mpipara
+  USE xios_mod
+  USE wxios
+  USE trace
+  USE distrib_icosa_lmdz_mod, ONLY : transfer_icosa_to_lmdz, transfer_lmdz_to_icosa
+  USE physics_external_mod, ONLY : it, f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q 
+  USE write_field_mod
+  USE checksum_mod
+! from LMDZ
+  USE mod_phys_lmdz_omp_data, ONLY: klon_omp
+  USE geometry_mod, ONLY : cell_area
+  USE physiq_mod, ONLY: physiq
+  IMPLICIT NONE
+  
+    REAL(rstd),POINTER :: phis(:)
+    REAL(rstd),POINTER :: ps(:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:,:)
+    REAL(rstd),POINTER :: u(:,:)
+    REAL(rstd),POINTER :: wflux(:,:)
+    REAL(rstd),POINTER :: q(:,:,:)
+    REAL(rstd),POINTER :: p(:,:)
+    REAL(rstd),POINTER :: pks(:)
+    REAL(rstd),POINTER :: pk(:,:)
+    REAL(rstd),POINTER :: p_layer(:,:)
+    REAL(rstd),POINTER :: theta(:,:)
+    REAL(rstd),POINTER :: phi(:,:)
+    REAL(rstd),POINTER :: Temp(:,:)
+    REAL(rstd),POINTER :: ulon(:,:)
+    REAL(rstd),POINTER :: ulat(:,:)
+    REAL(rstd),POINTER :: dulon(:,:)
+    REAL(rstd),POINTER :: dulat(:,:)
+    REAL(rstd),POINTER :: dTemp(:,:)
+    REAL(rstd),POINTER :: dq(:,:,:)
+    REAL(rstd),POINTER :: dps(:)
+    REAL(rstd),POINTER :: duc(:,:,:)
+
+
+    INTEGER :: ind,l
+    
+    REAL(rstd),ALLOCATABLE,SAVE :: ps_phy(:)
+!$OMP THREADPRIVATE(ps_phy)
+    REAL(rstd),ALLOCATABLE,SAVE :: p_phy(:,:)
+!$OMP THREADPRIVATE(p_phy)
+    REAL(rstd),ALLOCATABLE,SAVE :: p_layer_phy(:,:)
+!$OMP THREADPRIVATE(p_layer_phy)
+    REAL(rstd),ALLOCATABLE,SAVE :: Temp_phy(:,:)
+!$OMP THREADPRIVATE(Temp_phy)
+    REAL(rstd),ALLOCATABLE,SAVE :: phis_phy(:)
+!$OMP THREADPRIVATE(phis_phy)
+    REAL(rstd),ALLOCATABLE,SAVE :: phi_phy(:,:)
+!$OMP THREADPRIVATE(phi_phy)
+    REAL(rstd),ALLOCATABLE,SAVE :: ulon_phy(:,:)
+!$OMP THREADPRIVATE(ulon_phy)
+    REAL(rstd),ALLOCATABLE,SAVE :: ulat_phy(:,:)
+!$OMP THREADPRIVATE(ulat_phy)
+    REAL(rstd),ALLOCATABLE,SAVE :: q_phy(:,:,:)
+!$OMP THREADPRIVATE(q_phy)
+    REAL(rstd),ALLOCATABLE,SAVE :: wflux_phy(:,:)
+!$OMP THREADPRIVATE(wflux_phy)
+    REAL(rstd),ALLOCATABLE,SAVE :: dulon_phy(:,:)
+!$OMP THREADPRIVATE(dulon_phy)
+    REAL(rstd),ALLOCATABLE,SAVE :: dulat_phy(:,:)
+!$OMP THREADPRIVATE(dulat_phy)
+    REAL(rstd),ALLOCATABLE,SAVE :: dTemp_phy(:,:)
+!$OMP THREADPRIVATE(dTemp_phy)
+    REAL(rstd),ALLOCATABLE,SAVE :: dq_phy(:,:,:)
+!$OMP THREADPRIVATE(dq_phy)
+    REAL(rstd),ALLOCATABLE,SAVE :: dps_phy(:)
+!$OMP THREADPRIVATE(dps_phy)
+    REAL(rstd)   :: dtphy 
+    LOGICAL      :: debut
+    LOGICAL      :: lafin
+    LOGICAL,SAVE :: first=.TRUE.
+!$OMP THREADPRIVATE(first)
+
+    
+    IF(first) THEN 
+      debut=.TRUE.
+    ELSE 
+      debut=.FALSE.
+    ENDIF
+
+
+    IF(it-itau0>=itaumax) THEN 
+      lafin=.TRUE.
+    ELSE 
+      lafin=.FALSE.
+    ENDIF
+
+    IF (first) THEN
+      first=.FALSE.
+      CALL init_message(f_u,req_e1_vect,req_u)
+      ALLOCATE(ps_phy(klon_omp))
+      ALLOCATE(p_phy(klon_omp,llm+1))
+      ALLOCATE(p_layer_phy(klon_omp,llm))
+      ALLOCATE(Temp_phy(klon_omp,llm))
+      ALLOCATE(phis_phy(klon_omp))
+      ALLOCATE(phi_phy(klon_omp,llm))
+      ALLOCATE(ulon_phy(klon_omp,llm))
+      ALLOCATE(ulat_phy(klon_omp,llm))
+      ALLOCATE(q_phy(klon_omp,llm,nqtot))
+      ALLOCATE(wflux_phy(klon_omp,llm))
+      ALLOCATE(dulon_phy(klon_omp,llm))
+      ALLOCATE(dulat_phy(klon_omp,llm))
+      ALLOCATE(dTemp_phy(klon_omp,llm))
+      ALLOCATE(dq_phy(klon_omp,llm,nqtot))
+      ALLOCATE(dps_phy(klon_omp))
+!$OMP BARRIER
+    ENDIF
+
+
+!$OMP MASTER        
+!    CALL update_calendar(it)
+!$OMP END MASTER
+!$OMP BARRIER
+    dtphy=itau_physics*dt
+    
+    
+    
+    CALL transfert_message(f_u,req_u)
+    
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      IF (assigned_domain(ind)) THEN
+        CALL swap_geometry(ind)
+      
+        phis=f_phis(ind)
+        ps=f_ps(ind)
+        theta_rhodz=f_theta_rhodz(ind)
+        u=f_u(ind)
+        q=f_q(ind)
+        wflux=f_wflux(ind)
+        p=f_p(ind)
+        pks=f_pks(ind)
+        pk=f_pk(ind)
+        p_layer=f_p_layer(ind)
+        theta=f_theta(ind)
+        phi=f_phi(ind)
+        Temp=f_Temp(ind)
+        ulon=f_ulon(ind)
+        ulat=f_ulat(ind)
+            
+        CALL grid_icosa_to_physics
+
+      ENDIF
+    ENDDO
+   
+!$OMP BARRIER
+!$OMP MASTER
+    CALL SYSTEM_CLOCK(start_clock)
+!$OMP END MASTER
+    CALL trace_start("physic")
+!    CALL trace_off()
+
+
+!    CALL writeField("p_in",f_p)
+!    CALL writeField("p_layer_in",f_p_layer)
+!    CALL writeField("phi_in",f_phi)
+!    CALL writeField("phis_in",f_phis)
+!    CALL writeField("ulon_in",f_ulon)
+!    CALL writeField("ulat_in",f_ulat)
+!    CALL writeField("Temp_in",f_Temp)
+!    CALL writeField("q_in",f_q)
+!    CALL writeField("wflux_in",f_wflux)
+
+!    CALL checksum(f_p)
+!    CALL checksum(f_p_layer)
+!    CALL checksum(f_phi)
+!    CALL checksum(f_phis)
+!    CALL checksum(f_ulon)
+!    CALL checksum(f_ulat)
+!    CALL checksum(f_Temp)
+!    CALL checksum(f_q)
+!    CALL checksum(f_wflux)
+
+    CALL transfer_icosa_to_lmdz(f_p      , p_phy)
+    CALL transfer_icosa_to_lmdz(f_p_layer, p_layer_phy)
+    CALL transfer_icosa_to_lmdz(f_phi    , phi_phy)
+    CALL transfer_icosa_to_lmdz(f_phis   , phis_phy )
+    CALL transfer_icosa_to_lmdz(f_ulon   , ulon_phy )
+    CALL transfer_icosa_to_lmdz(f_ulat   , ulat_phy)
+    CALL transfer_icosa_to_lmdz(f_Temp   , Temp_phy)
+    CALL transfer_icosa_to_lmdz(f_q      , q_phy)
+    CALL transfer_icosa_to_lmdz(f_wflux  , wflux_phy)
+
+    DO l=1,llm
+      ! Warning: In the physics, vertical flux convention is positive if downwards!
+      wflux_phy(:,l)= - wflux_phy(:,l)*cell_area(:)
+      ! Compute relative geopotential
+      phi_phy(:,l)=phi_phy(:,l)-phis_phy(:)
+    ENDDO
+    
+    CALL wxios_set_context()
+
+    ! Update pday and ptime to send to physics
+!$OMP MASTER
+    ptime=ptime+dtphy/day_length
+    IF (ptime >= 1.) THEN
+      ptime=ptime-1
+      pday=pday+1
+    ENDIF
+!$OMP END MASTER
+!$OMP BARRIER    
+    CALL physiq(klon_omp, llm, nqtot, &
+                debut, lafin, &
+                pday, ptime, dtphy, &
+                p_phy, p_layer_phy, phi_phy, &
+                ulon_phy, ulat_phy, Temp_phy, q_phy, &
+                wflux_phy, &
+                dulon_phy, dulat_phy, dTemp_phy, dq_phy, dps_phy)
+    
+    CALL transfer_lmdz_to_icosa(dulon_phy, f_dulon )
+    CALL transfer_lmdz_to_icosa(dulat_phy, f_dulat )
+    CALL transfer_lmdz_to_icosa(dTemp_phy, f_dTemp )
+    CALL transfer_lmdz_to_icosa(dq_phy   , f_dq )
+    CALL transfer_lmdz_to_icosa(dps_phy  , f_dps )
+ 
+!    CALL writeField("dulon_out",f_dulon)
+!    CALL writeField("dulat_out",f_dulat)
+!    CALL writeField("dTemp_out",f_dTemp)
+!    CALL writeField("dq_out",f_dq)
+!    CALL writeField("dps_out",f_dps)
+
+!    CALL checksum(f_dulon)
+!    CALL checksum(f_dulat)
+!    CALL checksum(f_dTemp)
+!    CALL checksum(f_dq)
+!    CALL checksum(f_dps)
+    
+    CALL send_message(f_dps,req_dps0)
+    CALL send_message(f_dulon,req_dulon0)
+    CALL send_message(f_dulat,req_dulat0)
+    CALL send_message(f_dTemp,req_dTemp0)
+    CALL send_message(f_dq,req_dq0)
+
+    CALL wait_message(req_dps0)
+    CALL wait_message(req_dulon0)
+    CALL wait_message(req_dulat0)
+    CALL wait_message(req_dTemp0)
+    CALL wait_message(req_dq0)
+
+
+!    CALL trace_on()
+    CALL trace_end("physic")
+!$OMP MASTER
+    CALL SYSTEM_CLOCK(stop_clock)
+    count_clock=count_clock+stop_clock-start_clock
+!$OMP END MASTER
+
+!$OMP BARRIER                       
+
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      IF (assigned_domain(ind)) THEN
+        CALL swap_geometry(ind)
+
+        theta_rhodz=f_theta_rhodz(ind)
+        u=f_u(ind)
+        q=f_q(ind)
+        ps=f_ps(ind)
+        dulon=f_dulon(ind)
+        dulat=f_dulat(ind)
+        Temp=f_temp(ind)
+        dTemp=f_dTemp(ind)
+        dq=f_dq(ind)
+        dps=f_dps(ind)
+        duc=f_duc(ind)
+        p=f_p(ind)
+        pks=f_pks(ind)
+        pk=f_pk(ind)
+      
+        CALL grid_physics_to_icosa
+      ENDIF
+    ENDDO
+
+!$OMP BARRIER
+    CALL xios_set_context    
+   
+ 
+  CONTAINS
+
+    SUBROUTINE grid_icosa_to_physics
+    USE pression_mod
+    USE exner_mod
+    USE theta2theta_rhodz_mod
+    USE geopotential_mod
+    USE wind_mod
+    USE omp_para
+    IMPLICIT NONE
+    
+    REAL(rstd) :: uc(3)
+    INTEGER :: i,j,ij,l
+    
+
+! compute pression
+
+      DO    l    = ll_begin,ll_endp1
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            p(ij,l) = ap(l) + bp(l) * ps(ij)
+          ENDDO
+        ENDDO
+      ENDDO
+
+!$OMP BARRIER
+
+! compute exner
+       
+       IF (is_omp_first_level) THEN
+         DO j=jj_begin,jj_end
+            DO i=ii_begin,ii_end
+               ij=(j-1)*iim+i
+               pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+            ENDDO
+         ENDDO
+       ENDIF
+
+       ! 3D : pk
+       DO l = ll_begin,ll_end
+          DO j=jj_begin,jj_end
+             DO i=ii_begin,ii_end
+                ij=(j-1)*iim+i
+                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
+             ENDDO
+          ENDDO
+       ENDDO
+
+!$OMP BARRIER
+
+!   compute theta, temperature and pression at layer
+    DO    l    = ll_begin, ll_end
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          theta(ij,l) = theta_rhodz(ij,l,1) / ((p(ij,l)-p(ij,l+1))/g)
+          Temp(ij,l) = theta(ij,l) * pk(ij,l) / cpp
+          p_layer(ij,l)=preff*(pk(ij,l)/cpp)**(1./kappa) 
+        ENDDO
+      ENDDO
+    ENDDO
+
+
+!!! Compute geopotential
+       
+  ! for first layer
+  IF (is_omp_first_level) THEN
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
+      ENDDO
+    ENDDO
+  ENDIF
+!!-> implicit flush on phi(:,1)
+          
+  ! for other layers
+  DO l = ll_beginp1, ll_end
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        phi(ij,l) =  0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
+                         * (  pk(ij,l-1) -  pk(ij,l)    )
+      ENDDO
+    ENDDO
+  ENDDO       
+
+!$OMP BARRIER
+
+
+  IF (is_omp_first_level) THEN
+    DO l = 2, llm
+      DO j=jj_begin,jj_end
+! ---> Bug compilo intel ici en openmp
+! ---> Couper la boucle
+       IF (j==jj_end+1) PRINT*,"this message must not be printed"
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          phi(ij,l) = phi(ij,l)+ phi(ij,l-1)
+        ENDDO
+      ENDDO
+    ENDDO
+! --> IMPLICIT FLUSH on phi --> non 
+  ENDIF 
+
+! compute wind centered lon lat compound
+    DO l=ll_begin,ll_end
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          uc(:)=1/Ai(ij)*                                                                                                &
+                        ( ne(ij,right)*u(ij+u_right,l)*le(ij+u_right)*((xyz_v(ij+z_rdown,:)+xyz_v(ij+z_rup,:))/2-centroid(ij,:))  &
+                         + ne(ij,rup)*u(ij+u_rup,l)*le(ij+u_rup)*((xyz_v(ij+z_rup,:)+xyz_v(ij+z_up,:))/2-centroid(ij,:))          &
+                         + ne(ij,lup)*u(ij+u_lup,l)*le(ij+u_lup)*((xyz_v(ij+z_up,:)+xyz_v(ij+z_lup,:))/2-centroid(ij,:))          &
+                         + ne(ij,left)*u(ij+u_left,l)*le(ij+u_left)*((xyz_v(ij+z_lup,:)+xyz_v(ij+z_ldown,:))/2-centroid(ij,:))    &
+                         + ne(ij,ldown)*u(ij+u_ldown,l)*le(ij+u_ldown)*((xyz_v(ij+z_ldown,:)+xyz_v(ij+z_down,:))/2-centroid(ij,:))&
+                         + ne(ij,rdown)*u(ij+u_rdown,l)*le(ij+u_rdown)*((xyz_v(ij+z_down,:)+xyz_v(ij+z_rdown,:))/2-centroid(ij,:)))
+          ulon(ij,l)=sum(uc(:)*elon_i(ij,:))
+          ulat(ij,l)=sum(uc(:)*elat_i(ij,:)) 
+        ENDDO
+      ENDDO
+    ENDDO
+
+!$OMP BARRIER
+    END SUBROUTINE grid_icosa_to_physics
+
+
+    SUBROUTINE grid_physics_to_icosa
+    USE theta2theta_rhodz_mod
+    USE omp_para
+    IMPLICIT NONE
+      INTEGER :: i,j,ij,l,iq
+          
+      DO l=ll_begin,ll_end
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            duc(ij,:,l)=dulon(ij,l)*elon_i(ij,:)+dulat(ij,l)*elat_i(ij,:)
+          ENDDO
+        ENDDO
+      ENDDO
+
+      DO l=ll_begin,ll_end
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            u(ij+u_right,l) = u(ij+u_right,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
+            u(ij+u_lup,l) = u(ij+u_lup,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
+            u(ij+u_ldown,l) = u(ij+u_ldown,l) + dtphy*sum( 0.5*(duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
+          ENDDO
+        ENDDO
+      ENDDO          
+
+      DO l=ll_begin,ll_end
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            Temp(ij,l)=Temp(ij,l)+ dtphy * dTemp(ij,l)
+          ENDDO
+        ENDDO
+      ENDDO          
+      
+      DO iq=1,nqtot
+        DO l=ll_begin,ll_end
+          DO j=jj_begin,jj_end
+            DO i=ii_begin,ii_end
+              ij=(j-1)*iim+i
+              q(ij,l,iq)=q(ij,l,iq)+ dtphy * dq(ij,l,iq)
+            ENDDO
+          ENDDO
+        ENDDO 
+      ENDDO
+
+!$OMP BARRIER
+      
+     IF (is_omp_first_level) THEN
+       DO j=jj_begin,jj_end
+         DO i=ii_begin,ii_end
+           ij=(j-1)*iim+i
+           ps(ij)=ps(ij)+ dtphy * dps(ij)
+          ENDDO
+       ENDDO
+     ENDIF
+
+!     CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
+
+! compute pression
+!$OMP BARRIER
+      DO    l    = ll_begin,ll_endp1
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            p(ij,l) = ap(l) + bp(l) * ps(ij)
+          ENDDO
+        ENDDO
+      ENDDO
+
+!$OMP BARRIER
+
+! compute exner
+       
+       IF (is_omp_first_level) THEN
+         DO j=jj_begin,jj_end
+            DO i=ii_begin,ii_end
+               ij=(j-1)*iim+i
+               pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+            ENDDO
+         ENDDO
+       ENDIF
+
+       ! 3D : pk
+       DO l = ll_begin,ll_end
+          DO j=jj_begin,jj_end
+             DO i=ii_begin,ii_end
+                ij=(j-1)*iim+i
+                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
+             ENDDO
+          ENDDO
+       ENDDO
+
+!$OMP BARRIER
+
+!   compute theta, temperature and pression at layer
+    DO    l    = ll_begin, ll_end
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          theta_rhodz(ij,l,1) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp )
+        ENDDO
+      ENDDO
+    ENDDO
+    
+    END SUBROUTINE grid_physics_to_icosa
+
+
+
+  END SUBROUTINE physics
+
+
+
+
+
+END MODULE interface_icosa_lmdz_mod
