MODULE vertical_remap_mod
  USE icosa
  

CONTAINS

   
  SUBROUTINE vertical_remap(pressure_level,field_in,f_ps,field_out)
  USE icosa
  USE pression_mod
  USE omp_para
  IMPLICIT NONE
    REAL(rstd), INTENT(IN) :: pressure_level(:)
    TYPE(t_field),POINTER :: field_in(:)
    TYPE(t_field),POINTER :: f_ps(:)
    TYPE(t_field),POINTER :: field_out(:)

    TYPE(t_field),POINTER,SAVE :: f_p(:)
    REAL(rstd),POINTER :: in(:,:)
    REAL(rstd),POINTER :: out(:,:)
    REAL(rstd),POINTER :: p(:,:)
    
    INTEGER :: ind

    CALL allocate_field(f_p,field_t,type_real,llm+1)
    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=field_in(ind)
      out=field_out(ind)
      CALL compute_vertical_remap(pressure_level,in,p,out)
    ENDDO
    
  END SUBROUTINE  vertical_remap

  SUBROUTINE compute_vertical_remap(pressure_level,in,p,out)
  USE omp_para
  IMPLICIT NONE
    REAL(rstd),INTENT(IN)  :: pressure_level(:)
    REAL(rstd),INTENT(IN)  :: in(:,:)
    REAL(rstd),INTENT(IN)  :: p(iim*jjm,llm+1)
    REAL(rstd),INTENT(OUT) :: out(iim*jjm,llm)
    REAL(rstd) :: coeff, pmid
    INTEGER :: i,j,ij,l,n,nb_level
    INTEGER :: a
    INTEGER :: b
    LOGICAL :: positive
    
    nb_level=size(pressure_level)
    IF (pressure_level(1)>=pressure_level(nb_level)) THEN 
      positive=.FALSE.
    ELSE 
      positive=.TRUE.
    ENDIF
      
 !$OMP BARRIER    
    IF (is_omp_level_master) THEN
    
    DO l=1,llm
      DO j=jj_begin,jj_end
        DO i=ii_begin,ii_end
          ij=(j-1)*iim+i
          pmid=0.5*(p(ij,l)+p(ij,l+1))
          IF (positive) THEN
            a=0
            DO n=1,nb_level-1
              IF ( (pmid>=pressure_level(n) .AND. pmid<pressure_level(n+1))) THEN
               a=n ; b=n+1 ; EXIT
              ENDIF
            ENDDO
            IF (a==0) THEN
              IF (pmid<=pressure_level(1)) THEN
                a=1 ; b=2
              ELSE
                a=nb_level-1 ; b=nb_level
              ENDIF
            ENDIF
          ELSE
            a=0
            DO n=1,nb_level-1
              IF ( (pmid<=pressure_level(n) .AND. pmid>pressure_level(n+1))) THEN
               a=n ; b=n+1 ; EXIT
              ENDIF
            ENDDO
            
            IF (a==0) THEN
              IF (pmid >= pressure_level(1)) THEN
                a=1 ; b=2
              ELSE
                a=nb_level-1 ; b=nb_level
              ENDIF
            ENDIF
          ENDIF
                  
          coeff=(pmid-pressure_level(a))/(pressure_level(a)-pressure_level(b))
          out(ij,l)=in(ij,a)+coeff*(in(ij,a)-in(ij,b))
        ENDDO
      ENDDO
    ENDDO
 
    ENDIF
 !$OMP BARRIER    
        
  END SUBROUTINE compute_vertical_remap

END MODULE vertical_remap_mod