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_pmid,f_in,f_out,pval)
  USE icosa
  USE pression_mod
  USE omp_para
  IMPLICIT NONE
    TYPE(t_field),POINTER :: f_pmid(:)
    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 :: pmid(:,:)
    
    INTEGER :: ind
        
    DO ind=1,ndomain
      IF (.NOT. assigned_domain(ind)) CYCLE
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      pmid=f_pmid(ind)
      in=f_in(ind)
      out=f_out(ind)
      CALL compute_vertical_interp(pmid,in,out,pval)
    ENDDO
    
  END SUBROUTINE  vertical_interp

  SUBROUTINE compute_vertical_interp(pmid,in,out,pval)
  USE omp_para
  IMPLICIT NONE
    REAL(rstd),INTENT(IN) :: pmid(iim*jjm,llm)
    REAL(rstd),INTENT(IN) :: in(iim*jjm,llm)
    REAL(rstd),INTENT(OUT) :: out(iim*jjm)
    REAL(rstd) :: pval, coeff
    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(pmid(ij,l)<pval .AND. l>1)
            l=l-1
          ENDDO
          coeff=(pval-pmid(ij,l))/(pmid(ij,l)-pmid(ij,l+1))
          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