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
