source: codes/icosagcm/trunk/src/disvert_strato.f90 @ 521

MODULE disvert_strato_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
  USE earth_const
  IMPLICIT NONE
  REAL(rstd),INTENT(OUT) :: ap(:)
  REAL(rstd),INTENT(OUT) :: bp(:)
  REAL(rstd),INTENT(OUT) :: presnivs(:)
  
  REAL(rstd) :: dsig(llm)
  REAL(rstd) :: sig(llm+1)
  REAL(rstd) :: snorm
  REAL(rstd) :: x
  INTEGER :: l
  
    snorm  = 0.
    DO l = 1, llm
      x = Pi*(l-0.5)/(llm+1)
      dsig(l) = (1.0 + 7.0 * SIN(x)**2) * (0.5* (1.-tanh( (x-Pi/2)/(Pi/2) ) ))**2
      snorm = snorm + dsig(l)
    ENDDO    
    
    DO l = 1, llm
      dsig(l) = dsig(l)/snorm
    ENDDO

    sig(llm+1) = 0.
    DO l = llm, 1, -1
      sig(l) = sig(l+1) + dsig(l)
    ENDDO

    bp(llm+1) =   0.
    DO l = 1, llm
      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
      ap(l) = pa * ( sig(l) - bp(l) )
    ENDDO
    bp(1)=1.
    ap(1)=0.
    ap(llm+1) = pa * ( sig(llm+1) - bp(llm+1) )
    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_strato_mod