source: codes/icosagcm/devel/src/vertical/disvert_dcmip200.f90 @ 915

  MODULE disvert_dcmip200_mod
    USE prec, ONLY : rstd
    USE grid_param, ONLY : llm

  IMPLICIT NONE
  PRIVATE

  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)

  PUBLIC :: init_disvert, ap, bp, presnivs

CONTAINS
!=========================================================================

  SUBROUTINE init_disvert
    USE mpipara
    
    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 earth_const, ONLY : cpp, kappa, preff, g
    USE icosa, ONLY       : ncar_dz, ncar_T0, ncar_p0
    USE getin_mod, ONLY   : getin
    USE mpipara, ONLY     : is_mpi_root
    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) ::  eta_top, eta
    REAL(rstd),PARAMETER :: N=0.01         ! Brunt-Vaisala frequency (s-1)
    REAL(rstd),PARAMETER :: T0=300.       ! Surface temperature at the equator (K)
    REAL(rstd) :: Rd       
    REAL(rstd), PARAMETER :: Gamma=0.0065       
    
    Rd=cpp*kappa
    ncar_dz=400 ; CALL getin('ncar_dz',ncar_dz);
    cindx=1 ; CALL getin('ncar_disvert_c',cindx)  
    
    eta_top = (1-Gamma/T0*llm*ncar_dz)**(g/(Rd*Gamma))  
    do l = 1,llm+1
       eta = (1-Gamma/T0*(l-1)*ncar_dz)**(g/(Rd*Gamma))
       PRINT *,'eta ->', eta
       bp(l) = ((eta - eta_top)/(1 - eta_top))**cindx
       ap(l) = preff * ( eta - bp(l) )
    ENDDO
    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_dcmip200_mod