MODULE vertical_interp_mod
  USE icosa
  PRIVATE
  
  TYPE(t_field),SAVE, POINTER :: f_p(:)
  

  PUBLIC init_vertical_interp, vertical_interp, compute_vertical_interp

CONTAINS

  SUBROUTINE init_vertical_interp
  USE icosa
  IMPLICIT NONE
    
    CALL allocate_field(f_p,field_t,type_real,llm+1)
    
  END SUBROUTINE init_vertical_interp
    
  SUBROUTINE vertical_interp(f_ps,f_in,f_out,pval)
  USE icosa
  USE pression_mod
  USE omp_para
  IMPLICIT NONE
    TYPE(t_field),POINTER :: f_ps(:)
    TYPE(t_field),POINTER :: f_in(:)
    TYPE(t_field),POINTER :: f_out(:)
    REAL(rstd),INTENT(IN) :: pval

    REAL(rstd),POINTER :: in(:,:)
    REAL(rstd),POINTER :: out(:)
    REAL(rstd),POINTER :: p(:,:)
    
    INTEGER :: ind
        
    CALL pression(f_ps,f_p)
 
    DO ind=1,ndomain
      IF (.NOT. assigned_domain(ind)) CYCLE
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      p=f_p(ind)
      in=f_in(ind)
      out=f_out(ind)
      CALL compute_vertical_interp(p,in,out,pval)
    ENDDO
    
  END SUBROUTINE  vertical_interp

  SUBROUTINE compute_vertical_interp(p,in,out,pval)
  USE omp_para
  IMPLICIT NONE
    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
    REAL(rstd),INTENT(IN) :: in(iim*jjm,llm)
    REAL(rstd),INTENT(OUT) :: out(iim*jjm)
    REAL(rstd) :: pval
    REAL(rstd) :: coeff, pmid,pmidp1
    INTEGER :: i,j,ij,l
        
!$OMP BARRIER    
    IF (is_omp_level_master) THEN
    
      DO j=jj_begin-1,jj_end+1
        DO i=ii_begin-1,ii_end+1
          ij=(j-1)*iim+i
          l=llm-1
          DO WHILE(0.5*(p(ij,l)+p(ij,l+1))<pval .AND. l>1)
            l=l-1
          ENDDO
          pmid=0.5*(p(ij,l)+p(ij,l+1))
          pmidp1=0.5*(p(ij,l+1)+p(ij,l+2))

          coeff=(pval-pmid)/(pmid-pmidp1)
        
          out(ij)=in(ij,l)+coeff*(in(ij,l)-in(ij,l+1))
        ENDDO
      ENDDO
   
    ENDIF
!$OMP BARRIER    
        
  END SUBROUTINE compute_vertical_interp

END MODULE vertical_interp_mod