MODULE interface_icosa_lmdz_mod

  USE field_mod, ONLY: t_field
  USE transfert_mod, ONLY: t_message 
  
 
  TYPE(t_message),SAVE :: req_u, req_z
  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_vort(:)   
  TYPE(t_field),POINTER,SAVE :: f_vortc(:)   
  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 :: start_clock
  INTEGER :: stop_clock
  INTEGER :: count_clock=0
  
  INTEGER,SAVE :: nbp_phys
  INTEGER,SAVE :: nbp_phys_glo


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
  USE nudging_mod, ONLY : lam_halo_scheme
  
! 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
  USE infotrac_phy, ONLY: init_infotrac_phy

   
  
  IMPLICIT NONE
  INTEGER  :: ind,i,j,ij,pos,h
  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(:,:)
  LOGICAL   ,ALLOCATABLE :: outside(:,:)
  LOGICAL   ,ALLOCATABLE :: outside_tmp(:,:)
  LOGICAL   ,POINTER     :: out(:,:)
!  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_vort,field_z,type_real,llm,name="vort_in")
    CALL allocate_field(f_vortc,field_t,type_real,llm,name="vortc_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)

      ALLOCATE(outside(ii_begin:ii_end,jj_begin:jj_end)) ! for limited area : don't take cells arround the border
      ALLOCATE(outside_tmp(ii_begin-1:ii_end+1,jj_begin-1:jj_end+1)) ! for limited area : don't take cells arround the border
      out=>domain(ind)%outside
      DO j=jj_begin,jj_end
        DO i=ii_begin,ii_end
          outside(i,j)=  out(i+1,j)     .OR. & ! right
                         out(i,j+1    ) .OR. & ! rup
                         out(i-1  ,j+1) .OR. & ! lup
                         out(i-1  ,j)   .OR. & !left
                         out(i    ,j-1) .OR. & !ldown
                         out(i+1,j-1)          !rdown    
        ENDDO
      ENDDO

      outside_tmp=.FALSE.
      outside_tmp(ii_begin:ii_end,jj_begin:jj_end)=outside
      
      DO h=1,lam_halo_scheme-1 ! do not compute physic on limited area halo
        DO j=jj_begin,jj_end
          DO i=ii_begin,ii_end
              outside(i,j) = outside_tmp(i,j)       .OR. &
                             outside_tmp(i+1,j)     .OR. & ! right
                             outside_tmp(i,j+1    ) .OR. & ! rup
                             outside_tmp(i-1  ,j+1) .OR. & ! lup
                             outside_tmp(i-1  ,j)   .OR. & !left
                             outside_tmp(i    ,j-1) .OR. & !ldown
                             outside_tmp(i+1,j-1)          !rdown
           ENDDO
        ENDDO
        outside_tmp(ii_begin:ii_end,jj_begin:jj_end)=outside
      ENDDO
      
      DO j=jj_begin,jj_end
        DO i=ii_begin,ii_end
          IF (domain(ind)%own(i,j) .AND. .NOT.outside(i,j)) nbp_phys=nbp_phys+1
        ENDDO
      ENDDO
      DEALLOCATE(outside)
      DEALLOCATE(outside_tmp)
    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

  CALL init_infotrac_phy

          
!$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 earth_const, ONLY: scale_height
  USE xios_mod
  USE time_mod , init_time_icosa=> init_time 
  USE omp_para

! from LMDZ
  USE mod_grid_phy_lmdz, ONLY : unstructured, klon_glo
  USE mod_phys_lmdz_para, ONLY: klon_omp, reduce_min_lmdz => reduce_min , gather_lmdz => gather , bcast_lmdz => bcast
  USE time_phylmdz_mod, ONLY: init_time_lmdz => init_time
  USE transfert_mod
!  USE physics_distribution_mod, ONLY : init_physics_distribution
  USE geometry_mod, ONLY : init_geometry
  USE vertical_layers_mod, ONLY : init_vertical_layers
  USE infotrac_phy, ONLY : init_infotrac_phy
  USE inifis_mod, ONLY : inifis 
  USE readTracFiles_mod, ONLY: trac_type, isot_type
  USE tracer_icosa_mod, ONLY : tracs
   USE readTracFiles_mod, ONLY: delPhase
!  USE phyaqua_mod, ONLY : iniaqua 
   
  
  IMPLICIT NONE



  INTEGER  :: ind,i,j,k,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_r(:)
  INTEGER,   ALLOCATABLE :: ind_cell_glo(:)
  INTEGER,   ALLOCATABLE :: ind_cell_glo_tot(:)
  INTEGER,   ALLOCATABLE :: cell_glo_tot(:)
  INTEGER :: ncell_glo_tot

  REAL(rstd) :: pseudoalt(llm)
  REAL(rstd) :: aps(llm)
  REAL(rstd) :: bps(llm)
  REAL(rstd) :: scaleheight

  INTEGER :: run_length  
  REAL :: day_length ! length of a day (s) ! SAVEd to be OpenMP shared <--- NO!!!!
  INTEGER :: annee_ref  
  INTEGER :: day_ref    
  INTEGER :: day_ini    
  REAL    :: start_time 
  REAL    :: physics_timestep    

  ! Tracer stuff (SAVEd when needed to be OpenMP shared)
  INTEGER :: nq
  INTEGER                       :: nqo, nbtr, nbtr_inca
  CHARACTER(len=256)              :: type_trac
  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, ALLOCATABLE, DIMENSION(:) :: hadv_inca  ! index of horizontal trasport schema
  INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv_inca  ! index of vertical trasport schema

   TYPE(trac_type) ::  tracers_ico2lmdz(nqtot)    !=== TRACERS DESCRIPTORS VECTOR
   TYPE(isot_type) :: isotopes_ico2lmdz(1)        !=== ISOTOPES PARAMETERS VECTOR
   INTEGER :: iq 

   CHARACTER(LEN=3)      :: descrq(30)            !--- Advection scheme description tags
   logical, save :: first = .TRUE. 


    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_r(klon_omp))
    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_r)
   CALL deallocate_field(f_ind_cell_glo)
   ind_cell_glo=INT(ind_cell_glo_r)
   DEALLOCATE(ind_cell_glo_r)   
   
    
   CALL reduce_min_lmdz(MINVAL(-ind_cell_glo),ncell_glo_tot) ! reduce_max does not exist in lmdz, use reduce_min
   CALL bcast_lmdz(ncell_glo_tot)
   ncell_glo_tot=-ncell_glo_tot
   ALLOCATE(cell_glo_tot(0:ncell_glo_tot))
   ALLOCATE(ind_cell_glo_tot(klon_glo))
   CALL gather_lmdz(ind_cell_glo,ind_cell_glo_tot)
   CALL bcast_lmdz(ind_cell_glo_tot)
   
   cell_glo_tot=-1
   DO i=1,klon_glo
     cell_glo_tot(ind_cell_glo_tot(i))= 0
   ENDDO
   
   pos=0
   DO i=0,ncell_glo_tot
     IF (cell_glo_tot(i)/=-1) THEN
       cell_glo_tot(i) = pos
       pos=pos + 1
     ENDIF
   ENDDO

   DO i=1,klon_omp
     ind_cell_glo(i)=cell_glo_tot(ind_cell_glo(i))
   ENDDO  

   ind_cell_glo = ind_cell_glo + 1 ! lmdz expect global indices begining to 1 not 0
 
   
!   CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_mpi,1,MPI_INTEGER,comm_icosa,ierr)
!
!   displ(0)=0
!   DO i=1,mpi_size-1
!     displ(i)=displ(i-1)+ncell_mpi(i-1)
!   ENDDO
!   
!   ALLOCATE(ind_glo_tot(ncell_tot))
!   ALLOCATE(cell_glo_tot(0:ncell_glo_tot-1))
!   
!   cell_glo_tot(:)= -1
!   CALL MPI_ALLGATHERV(ind_glo, ncell, MPI_INTEGER, ind_glo_tot, ncell_mpi, displ, MPI_INTEGER, comm_icosa,ierr)
!   
!   DO i=1,ncell_tot
!     cell_glo_tot(ind_glo_tot(i))= 0
!   ENDDO
! 
!   ncell_glo=0
!   DO i=0,ncell_glo_tot-1
!     IF (cell_glo_tot(i)/=-1) THEN
!       cell_glo_tot(i) = ncell_glo
!       ncell_glo=ncell_glo + 1
!     ENDIF
!   ENDDO
!
!   DO i=1,ncell
!     ind_glo(i)=cell_glo_tot(ind_glo(i))
!   ENDDO  

              
    CALL init_geometry(klon_omp,lonfi, latfi, bounds_lonfi, bounds_latfi, airefi, 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,presinter,pseudoalt)

    ! Initialize tracer names, numbers, etc. for physics
    !Config  Key  = type_trac
    !Config  Desc = Choix de couplage avec model de chimie INCA ou REPROBUS
    !Config  Def  = lmdz
    !Config  Help = 
    !Config  'lmdz' = pas de couplage, pur LMDZ
    !Config  'lmdz|inca' = model de chime INCA 
    !Config  'lmdz|repr' = model de chime REPROBUS
    type_trac = 'lmdz'
    CALL getin('type_trac',type_trac)

    descrq( 1: 2) = ['LMV','BAK']
    descrq(10:20) = ['VL1','VLP','FH1','FH2','VLH','   ','PPM','PPS','PPP','   ','SLP']
    descrq(30)    =  'PRA'

    nqo = 0 
    DO iq=1,nqtot

       tracers_ico2lmdz(iq)%name = tracs(iq)%name
       
       tracers_ico2lmdz(iq)%gen0Name = tracs(iq)%name
       tracers_ico2lmdz(iq)%phase = tracs(iq)%phase
       tracers_ico2lmdz(iq)%iadv = tracs(iq)%iadv
          
       IF (tracs(iq)%component .eq. "dynamico") then 
          tracers_ico2lmdz(iq)%component='lmdz'
       ELSE
          tracers_ico2lmdz(iq)%component=tracs(iq)%component
       ENDIF

       if (tracers_ico2lmdz(iq)%iadv .ne. 0 ) tracers_ico2lmdz(iq)%isAdvected=.true.
        
       tracers_ico2lmdz(iq)%longName   = tracers_ico2lmdz(iq)%name
       IF(tracers_ico2lmdz(iq)%iadv > 0) tracers_ico2lmdz(iq)%longName=TRIM(tracers_ico2lmdz(iq)%name)//descrq(tracers_ico2lmdz(iq)%iadv)

       tracers_ico2lmdz(iq)%isInPhysics= delPhase(tracers_ico2lmdz(iq)%gen0Name) /= 'H2O' .OR. tracers_ico2lmdz(iq)%component /= 'lmdz' 
       tracers_ico2lmdz(iq)%iGeneration = 0 

    ENDDO

    nqo = COUNT(delPhase(tracs(:)%name) == 'H2O' .AND. tracers_ico2lmdz(:)%component == 'lmdz') !--- Number of water phases

    isotopes_ico2lmdz(1)%parent='H2O'
    isotopes_ico2lmdz(1)%phase='gls'
    isotopes_ico2lmdz(1)%nphas=3
     

    nbtr=nqtot-nqo 

    ALLOCATE(conv_flg(nbtr))
    ALLOCATE(pbl_flg(nbtr))

    conv_flg(:) = 1 ! convection activated for all tracers
    pbl_flg(:)  = 1 ! boundary layer activated for all tracers
       

    CALL init_infotrac_phy

   ! Initialize physical constant
    day_length=86400
    CALL getin('day_length',day_length)
    CALL inifis(day_length,radius,g,kappa*cpp,cpp)
 

    
  ! 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, int(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

  

#ifdef INCA
    CONTAINS 

      SUBROUTINE init_chem_trac() 
        IMPLICIT NONE

        CALL  Init_chem_inca_trac(nbtr)

      END SUBROUTINE init_chem_trac

      SUBROUTINE init_chem_transport()

        IMPLICIT NONE

        CALL init_transport(solsym, conv_flg,pbl_flg, hadv_inca, vadv_inca)

      END SUBROUTINE init_chem_transport


#else 
    CONTAINS 
      SUBROUTINE init_chem_trac() 
        IMPLICIT NONE

      END SUBROUTINE init_chem_trac

      SUBROUTINE init_chem_transport()

        IMPLICIT NONE

      END SUBROUTINE init_chem_transport


#endif



  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
  USE vorticity_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 :: vort(:,:)
    REAL(rstd),POINTER :: vortc(:,:)
    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 :: rot_phy(:,:)
!$OMP THREADPRIVATE(rot_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)
      CALL init_message(f_vort,req_z1_scal,req_z)
      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(rot_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
      IF (assigned_domain(ind)) THEN
        CALL swap_dimensions(ind)
        CALL swap_geometry(ind)
        u=f_u(ind)
        vort=f_vort(ind)
        CALL compute_vorticity(u,vort)
      ENDIF
    ENDDO

    CALL transfert_message(f_vort,req_z)

    
    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)
        vort=f_vort(ind)
        vortc=f_vortc(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 writeField("vortc",f_vortc)

!    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_vortc   , rot_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
      wflux_phy(:,l) = - wflux_phy(:,l)*cell_area(:)
      phi_phy(:,l)=phi_phy(:,l)-phis_phy(:)
    ENDDO
    
    CALL wxios_set_context()
 
    ! Ehouarn: rot_phy() not implemented!! Set it to zero for now
!    rot_phy(:,:)=0
    CALL physiq(klon_omp, llm, debut, lafin, dtphy, &
                p_phy, p_layer_phy, phi_phy, phis_phy, presnivs, &
                ulon_phy, ulat_phy, rot_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)

!$OMP BARRIER
          
  ! 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


! compute centered vorticity
   
    DO l=ll_begin,ll_end
      DO j=jj_begin,jj_end
        DO i=ii_begin,ii_end
          ij=(j-1)*iim+i
          vortc(ij,l) =  Riv(ij,vup)    * vort(ij+z_up,l)    + & 
                         Riv(ij,vlup)  * vort(ij+z_lup,l)   + & 
                         Riv(ij,vldown)* vort(ij+z_ldown,l) + & 
                         Riv(ij,vdown) * vort(ij+z_down,l)  + & 
                         Riv(ij,vrdown)* vort(ij+z_rdown,l) + & 
                         Riv(ij,vrup)  * vort(ij+z_rup,l) 
      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