MODULE physics_lmdz5_mod
  USE field_mod
  USE transfert_mod
  
  INTEGER,SAVE :: nbp_phys
  TYPE(t_message) :: req_u

  TYPE(t_field),POINTER :: f_p(:) 
  TYPE(t_field),POINTER :: f_pks(:)  
  TYPE(t_field),POINTER :: f_pk(:)  
  TYPE(t_field),POINTER :: f_p_layer(:)    
  TYPE(t_field),POINTER :: f_theta(:)   
  TYPE(t_field),POINTER :: f_phi(:)   
  TYPE(t_field),POINTER :: f_Temp(:)    
  TYPE(t_field),POINTER :: f_ulon(:)    
  TYPE(t_field),POINTER :: f_ulat(:)   
  TYPE(t_field),POINTER :: f_dulon(:)
  TYPE(t_field),POINTER :: f_dulat(:)
  TYPE(t_field),POINTER :: f_dTemp(:)
  TYPE(t_field),POINTER :: f_dq(:)
  TYPE(t_field),POINTER :: f_dps(:)
  TYPE(t_field),POINTER :: f_duc(:)

  INTEGER :: start_clock
  INTEGER :: stop_clock
  INTEGER :: count_clock=0

CONTAINS

  SUBROUTINE init_physics
  USE icosa
  USE domain_mod
  USE dimensions
  USE mpi_mod
  USE mpipara
  IMPLICIT NONE
  INTEGER  :: distrib(0:mpi_size-1)
  INTEGER  :: ind,i,j,ij,pos

  REAL(rstd),ALLOCATABLE :: latfi(:)
  REAL(rstd),ALLOCATABLE :: lonfi(:)
  REAL(rstd),ALLOCATABLE :: airefi(:)
   

    CALL allocate_field(f_p,field_t,type_real,llm+1)
    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)
    CALL allocate_field(f_theta,field_t,type_real,llm)
    CALL allocate_field(f_phi,field_t,type_real,llm)
    CALL allocate_field(f_Temp,field_t,type_real,llm)
    CALL allocate_field(f_ulon,field_t,type_real,llm)
    CALL allocate_field(f_ulat,field_t,type_real,llm)
    CALL allocate_field(f_dulon,field_t,type_real,llm)
    CALL allocate_field(f_dulat,field_t,type_real,llm)
    CALL allocate_field(f_dTemp,field_t,type_real,llm)
    CALL allocate_field(f_dq,field_t,type_real,llm,nqtot)
    CALL allocate_field(f_dps,field_t,type_real)
    CALL allocate_field(f_duc,field_t,type_real,3,llm)    
    
    
    nbp_phys=0
    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      nbp_phys=nbp_phys+ii_nb*jj_nb
    ENDDO
    
    CALL MPI_ALLGATHER(nbp_phys,1,MPI_INTEGER,distrib,1,MPI_INTEGER,comm_icosa,ierr)
    
    ALLOCATE(latfi(nbp_phys))
    ALLOCATE(lonfi(nbp_phys))
    ALLOCATE(airefi(nbp_phys))
    
    pos=0
    DO ind=1,ndomain
      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
          pos=pos+1
          CALL xyz2lonlat(xyz_i(ij,:),lonfi(pos),latfi(pos))
          airefi(pos)=Ai(ij)  
        ENDDO
      ENDDO
    ENDDO
    
    CALL init_gcm_lmdz(nbp_phys,mpi_size,distrib,latfi,lonfi,airefi)
    
  END SUBROUTINE init_physics
  
  SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)
  USE ICOSA
  USE time_mod
  USE disvert_mod
  USE transfert_mod
  IMPLICIT NONE
    INTEGER,INTENT(IN)    :: it
    TYPE(t_field),POINTER :: f_phis(:)
    TYPE(t_field),POINTER :: f_ps(:)
    TYPE(t_field),POINTER :: f_theta_rhodz(:)
    TYPE(t_field),POINTER :: f_u(:)
    TYPE(t_field),POINTER :: f_wflux(:)
    TYPE(t_field),POINTER :: f_q(:)
  
    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
    
    REAL(rstd),ALLOCATABLE,SAVE :: ps_phy(:)
    REAL(rstd),ALLOCATABLE,SAVE :: p_phy(:,:)
    REAL(rstd),ALLOCATABLE,SAVE :: p_layer_phy(:,:)
    REAL(rstd),ALLOCATABLE,SAVE :: Temp_phy(:,:)
    REAL(rstd),ALLOCATABLE,SAVE :: phis_phy(:)
    REAL(rstd),ALLOCATABLE,SAVE :: phi_phy(:,:)
    REAL(rstd),ALLOCATABLE,SAVE :: ulon_phy(:,:)
    REAL(rstd),ALLOCATABLE,SAVE :: ulat_phy(:,:)
    REAL(rstd),ALLOCATABLE,SAVE :: q_phy(:,:,:)
    REAL(rstd),ALLOCATABLE,SAVE :: wflux_phy(:,:)
    REAL(rstd),ALLOCATABLE,SAVE :: dulon_phy(:,:)
    REAL(rstd),ALLOCATABLE,SAVE :: dulat_phy(:,:)
    REAL(rstd),ALLOCATABLE,SAVE :: dTemp_phy(:,:)
    REAL(rstd),ALLOCATABLE,SAVE :: dq_phy(:,:,:)
    REAL(rstd),ALLOCATABLE,SAVE :: dps_phy(:)
    REAL :: dtphy 
    INTEGER,SAVE :: offset
!$OMP THREADPRIVATE(offset)
    LOGICAL :: lafin=.FALSE.
    LOGICAL,SAVE :: first=.TRUE.
!$OMP THREADPRIVATE(first)

    
    IF (first) THEN
      first=.FALSE.
      CALL init_message(f_u,req_e1_vect,req_u)
!$OMP MASTER
      ALLOCATE(ps_phy(nbp_phys))
      ALLOCATE(p_phy(nbp_phys,llm+1))
      ALLOCATE(p_layer_phy(nbp_phys,llm))
      ALLOCATE(Temp_phy(nbp_phys,llm))
      ALLOCATE(phis_phy(nbp_phys))
      ALLOCATE(phi_phy(nbp_phys,llm))
      ALLOCATE(ulon_phy(nbp_phys,llm))
      ALLOCATE(ulat_phy(nbp_phys,llm))
      ALLOCATE(q_phy(nbp_phys,llm,nqtot))
      ALLOCATE(wflux_phy(nbp_phys,llm))
      ALLOCATE(dulon_phy(nbp_phys,llm))
      ALLOCATE(dulat_phy(nbp_phys,llm))
      ALLOCATE(dTemp_phy(nbp_phys,llm))
      ALLOCATE(dq_phy(nbp_phys,llm,nqtot))
      ALLOCATE(dps_phy(nbp_phys))
!$OMP END MASTER
!$OMP BARRIER
    ENDIF

!$OMP MASTER        
    CALL update_calendar(it)
!$OMP END MASTER
!$OMP BARRIER
    dtphy=itau_physics*dt
    IF(it==itaumax) THEN 
      lafin=.TRUE.
    ELSE 
      lafin=.FALSE.
    ENDIF
    CALL transfert_message(f_u,req_u)
    offset=0
    
    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      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
      
    ENDDO

!$OMP BARRIER
!$OMP MASTER
    CALL SYSTEM_CLOCK(start_clock)
!$OMP END MASTER
    CALL calfis_icosa(current_time, dtphy, lafin, presnivs,      &
                      p_phy, p_layer_phy, phi_phy, phis_phy, ulon_phy, ulat_phy, Temp_phy, q_phy, wflux_phy, &
                      dulon_phy, dulat_phy, dTemp_phy, dq_phy, dps_phy  )


!$OMP MASTER
    CALL SYSTEM_CLOCK(stop_clock)
    count_clock=count_clock+stop_clock-start_clock
!$OMP END MASTER

!$OMP BARRIER                       
    offset=0 

    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      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
    ENDDO

!$OMP BARRIER
     
  CONTAINS
   
    SUBROUTINE grid_physics_to_icosa
    USE theta2theta_rhodz_mod
    USE omp_para
    IMPLICIT NONE
      INTEGER :: i,j,ij,l,iq
          
      DO j=jj_begin,jj_end
        DO i=ii_begin,ii_end
          ij=(j-1)*iim+i
          offset=offset+1
          
          dulon(ij,ll_begin:ll_end)=dulon_phy(offset,ll_begin:ll_end)
          dulat(ij,ll_begin:ll_end)=dulat_phy(offset,ll_begin:ll_end)
          dTemp(ij,ll_begin:ll_end)=dTemp_phy(offset,ll_begin:ll_end)
          Temp(ij,ll_begin:ll_end) = Temp_phy(offset,ll_begin:ll_end)
          dq(ij,ll_begin:ll_end,:)=dq_phy(offset,ll_begin:ll_end,:)
          if (omp_first) dps(ij)=dps_phy(offset)
        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
            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 (omp_first) 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 (omp_first) 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) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp )
        ENDDO
      ENDDO
    ENDDO
    
    END SUBROUTINE grid_physics_to_icosa
    
    
    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 (omp_first) 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) / ((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 (omp_first) 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

   DO l = 2, llm
!$OMP DO
     DO j=jj_begin,jj_end
       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



! 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


      
      DO j=jj_begin,jj_end
        DO i=ii_begin,ii_end
          ij=(j-1)*iim+i
          offset=offset+1

          IF (omp_first) ps_phy(offset) = ps(ij)
          p_phy(offset,ll_begin:ll_endp1) = p(ij,ll_begin:ll_endp1)
          p_layer_phy(offset,ll_begin:ll_end) = p_layer(ij,ll_begin:ll_end)
          Temp_phy(offset,ll_begin:ll_end) = Temp(ij,ll_begin:ll_end)
          IF (omp_first) phis_phy(offset) = phis(ij)
          phi_phy(offset,ll_begin:ll_end) = phi(ij,ll_begin:ll_end)-phis(ij)
          ulon_phy(offset,ll_begin:ll_end) = ulon(ij,ll_begin:ll_end)
          ulat_phy(offset,ll_begin:ll_end) = ulat(ij,ll_begin:ll_end)
          q_phy(offset,ll_begin:ll_end,:) = q(ij,ll_begin:ll_end,:)
          wflux_phy(offset,ll_begin:ll_end) = wflux(ij,ll_begin:ll_end)
        ENDDO
      ENDDO
     
     END SUBROUTINE grid_icosa_to_physics

  END SUBROUTINE physics
  
  
END MODULE physics_lmdz5_mod