MODULE disvert_apbp_mod
  USE icosa
  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: ap(:)
!$OMP THREADPRIVATE(ap)
  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: bp(:)
!$OMP THREADPRIVATE(bp)
  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:)
!$OMP THREADPRIVATE(presnivs)

CONTAINS

  SUBROUTINE init_disvert
!  USE icosa
  IMPLICIT NONE
  
    ALLOCATE(ap(llm+1))
    ALLOCATE(bp(llm+1))
    ALLOCATE(presnivs(llm))

    CALL disvert(ap,bp,presnivs)    

  END SUBROUTINE init_disvert  


  SUBROUTINE disvert(ap,bp,presnivs)
!  USE icosa
  USE mpipara, ONLY: is_mpi_root
  USE omp_para, ONLY: omp_in_parallel
  USE transfert_omp_mod, ONLY: bcast_omp
  USE free_unit_mod, ONLY : free_unit
  IMPLICIT NONE
  REAL(rstd),INTENT(OUT) :: ap(:)
  REAL(rstd),INTENT(OUT) :: bp(:)
  REAL(rstd),INTENT(OUT) :: presnivs(:)
  INTEGER :: unit
  CHARACTER(len=255) :: filename
  INTEGER :: l,ok
  
    filename="apbp.txt" !file to read
    ! but users may want to use some other file name
    CALL getin('read_apbp_file',filename)
    
!$OMP MASTER
    unit=free_unit()
    OPEN(unit,file=filename,status="old",action="read",iostat=ok)
    IF (ok/=0) THEN
      WRITE(*,*) "disvert_ap_bp error: input file ",trim(filename)," not found!"
      CALL ABORT
    ENDIF
    ! read in ap() and b() line by line, starting from surface up
    ! to model top
    DO l=1,llm+1
      READ(unit,fmt=*,iostat=ok) ap(l),bp(l)
      IF (ok/=0) THEN
        WRITE(*,*) "disvert_ap_bp error: failed reading ap(l) and bp(l) for l=",l
        CALL ABORT
      ENDIF
    ENDDO
    
    CLOSE(unit)
!$OMP END MASTER
    IF (omp_in_parallel()) THEN
      CALL bcast_omp(ap)
      CALL bcast_omp(bp)
    ENDIF
    
    ! build presnivs(), approximative mid-layer pressures
    DO l=1,llm
      presnivs(l)=0.5*(ap(l)+bp(l)*preff+ap(l+1)+bp(l+1)*preff)
    ENDDO
    
    ! tell the world about it
    IF (is_mpi_root) THEN
!$OMP MASTER
      WRITE(*,*) "ap()=",ap
      WRITE(*,*) "bp()=",bp
      WRITE(*,*) "Approximative mid-layer pressure, assuming a surface pressure preff=",preff," Pa"
      WRITE(*,*) "and approximative mid-layer height, assuming an atmospheric scale height of ",scale_height/1000," (km)"
      DO l=1,llm
        WRITE(*,*) 'PRESNIVS(',l,')=',presnivs(l),'  Z ~ ',log(preff/presnivs(l))*scale_height/1000,       &
                   ' DZ ~ ',scale_height/1000*log((ap(l)+bp(l)*preff)/ max(ap(l+1)+bp(l+1)*preff,1.e-10))
      ENDDO
!$OMP END MASTER
    ENDIF
  
  END SUBROUTINE disvert
  
END MODULE disvert_apbp_mod