MODULE disvert_ncarl30_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)
  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presinter(:)
  !$OMP THREADPRIVATE(presinter)

CONTAINS

  SUBROUTINE init_disvert
  USE icosa
  USE mpipara
  IMPLICIT NONE
  
    IF (llm/=30) STOP 'ERROR, the initialisation of this vertical discretisation must be done with 30 levels'
    ALLOCATE(ap(llm+1))
    ALLOCATE(bp(llm+1))
    ALLOCATE(presnivs(llm))
    ALLOCATE(presinter(llm+1))
    CALL disvert(ap,bp,presnivs,presinter)

  END SUBROUTINE init_disvert  


  SUBROUTINE disvert(ap,bp,presnivs,presinter)
  USE icosa
  USE mpipara
  IMPLICIT NONE
  REAL(rstd),INTENT(OUT) :: ap(:)
  REAL(rstd),INTENT(OUT) :: bp(:)
  REAL(rstd),INTENT(OUT) :: presnivs(:)
  REAL(rstd),INTENT(OUT) :: presinter(:)
  INTEGER :: l
  REAL(rstd) :: a(31)=(/ 0.00225523952394724, &
                         0.00503169186413288, &
                         0.0101579474285245 , &
                         0.0185553170740604 , &
                         0.0306691229343414 , &
                         0.0458674766123295 , &
                         0.0633234828710556 , &
                         0.0807014182209969 , &
                         0.0949410423636436 , &
                         0.11169321089983   , &
                         0.131401270627975  , &
                         0.154586806893349  , &
                         0.181863352656364  , &
                         0.17459799349308   , &
                         0.166050657629967  , &
                         0.155995160341263  , &
                         0.14416541159153   , &
                         0.130248308181763  , &
                         0.113875567913055  , &
                         0.0946138575673103 , &
                         0.0753444507718086 , &
                         0.0576589405536652 , &
                         0.0427346378564835 , &
                         0.0316426791250706 , &
                         0.0252212174236774 , &
                         0.0191967375576496 , &
                         0.0136180268600583 , &
                         0.00853108894079924, &
                         0.00397881818935275, &
                         0.                 , &
                         0. /)  
                         
  REAL(rstd) :: b(31)=(/ 0.                 , &
                         0.                 , &
                         0.                 , &
                         0.                 , &
                         0.                 , &
                         0.                 , &
                         0.                 , &
                         0.                 , &
                         0.                 , &
                         0.                 , &
                         0.                 , &
                         0.                 , &
                         0.                 , &
                         0.0393548272550106 , &
                         0.0856537595391273 , &
                         0.140122056007385  , &
                         0.204201176762581  , &
                         0.279586911201477  , &
                         0.368274360895157  , &
                         0.47261056303978   , &
                         0.576988518238068  , &
                         0.672786951065063  , &
                         0.753628432750702  , &
                         0.813710987567902  , &
                         0.848494648933411  , &
                         0.881127893924713  , &
                         0.911346435546875  , &
                         0.938901245594025  , &
                         0.963559806346893  , &
                         0.985112190246582  , &
                         1. /)                         
                    
    ap(:)=a(31:1:-1)*preff
    bp(:)=b(31:1:-1) 
    DO l = 1, llm
      presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )

!$OMP MASTER  
     IF (is_mpi_root)  PRINT*, 'PRESNIVS(',l,')=',presnivs(l),'  Z ~ ',log(preff/presnivs(l))*8.,       &
                               ' DZ ~ ',8.*log((ap(l)+bp(l)*preff)/ max(ap(l+1)+bp(l+1)*preff,1.e-10))
!$OMP END MASTER
    ENDDO
    DO l=1, llm+1
      presinter(l)= ap(l)+bp(l)*preff
    END DO
    
  END SUBROUTINE disvert
  
END MODULE disvert_ncarl30_mod
