MODULE disvert_dcmip31_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
  IMPLICIT NONE
  REAL(rstd),INTENT(OUT) :: ap(:)
  REAL(rstd),INTENT(OUT) :: bp(:)
  REAL(rstd),INTENT(OUT) :: presnivs(:)
  ! reads from run.def : ncar_dz, ncar_T0, ncar_p0, ncar_disvert_c
  INTEGER :: l,cindx
  REAL(rstd) :: GG, eta_top, eta
  REAL(rstd),PARAMETER :: N=0.01         ! Brunt-Vaisala frequency (s-1)
  REAL(rstd),PARAMETER :: Teq=300.       ! Surface temperature at the equator (K)
  REAL(rstd),PARAMETER :: Peq=1e5        ! Reference surface pressure at the equator (hPa)

  ncar_dz=400 ; CALL getin('ncar_dz',ncar_dz);
  cindx=1 ; CALL getin('ncar_disvert_c',cindx)

  GG=(g/N)**2/cpp
  
  eta_top = (GG/Teq*exp(-N**2*llm*ncar_dz/g)+1-GG/Teq)**(1./kappa)
  IF (is_mpi_root) PRINT *,'eta_top ->', eta_top
  do l = 1,llm+1
     eta = (GG/Teq*exp(-N**2*(l-1)*ncar_dz/g)+1-GG/Teq)**(1./kappa)
     IF (is_mpi_root) PRINT *,'eta ->', eta
     bp(l) = ((eta - eta_top)/(1 - eta_top))**cindx
     ap(l) = preff * ( eta - bp(l) )
  ENDDO
  IF (is_mpi_root) PRINT *,'eta ->', eta
  bp(1)=1.
  ap(1)=0.
  bp(llm+1) = 0
  
  DO l = 1, llm
     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
  ENDDO

!$OMP MASTER  
  IF (is_mpi_root) PRINT *, 'Vertical placement of model levels according to DCMIP Appendix E.3'
  IF (is_mpi_root) PRINT *, 'Parameters : ncar_dz=', ncar_dz, '  ncar_p0=',ncar_p0, '  ncar_disvert_c=',cindx
  IF (is_mpi_root) PRINT *, 'Isothermal amtosphere with ncar_T0=',ncar_T0 
!$OMP END MASTER  

END SUBROUTINE disvert

END  MODULE disvert_dcmip31_mod