source: LMDZ.3.3/trunk/libf/dyn3d/disvert0.F @ 3330

      SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)

c    Auteur :  P. Le Van .
c
      IMPLICIT NONE

#include "dimensions.h"
#include "paramet.h"
c
c=======================================================================
c
c
c    s = sigma ** kappa   :  coordonnee  verticale
c    dsig(l)            : epaisseur de la couche l ds la coord.  s
c    sig(l)             : sigma a l'interface des couches l et l-1
c    ds(l)              : distance entre les couches l et l-1 en coord.s
c
c=======================================================================
c
      REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
      REAL pa,preff,presnivs(llm)
c
c   declarations:
c   -------------
c
      REAL sig(llm+1),dsig(llm)
c
      INTEGER l
      REAL snorm
      REAL alpha,beta,gama,delta,deltaz,h
      INTEGER np,ierr
      REAL pi,x

      REAL SSUM
      EXTERNAL SSUM
c
c-----------------------------------------------------------------------
c
      pi=2.*ASIN(1.)

      OPEN(99,file='sigma.def',status='old',form='formatted',
     s   iostat=ierr)

c-----------------------------------------------------------------------
c   cas 1 on lit les options dans sigma.def:
c   ----------------------------------------

      IF (ierr.eq.0) THEN

      READ(99,*) deltaz
      READ(99,*) h
      READ(99,*) beta
      READ(99,*) gama
      READ(99,*) delta
      READ(99,*) np
      CLOSE(99)
      alpha=deltaz/(llm*h)
c

       DO 1  l = 1, llm
       dsig(l) = (alpha+(1.-alpha)*exp(-beta*(llm-l)))*
     $          ( (tanh(gama*l)/tanh(gama*llm))**np +
     $            (1.-l/FLOAT(llm))*delta )
   1   CONTINUE

       sig(1)=1.
       DO 101 l=1,llm-1
          sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
101    CONTINUE
       sig(llm+1)=0.

       DO 2  l = 1, llm
       dsig(l) = sig(l)-sig(l+1)
   2   CONTINUE
c

      ELSE
c-----------------------------------------------------------------------
c   cas 2 ancienne discretisation (LMD5...):
c   ----------------------------------------

      PRINT*,'WARNING!!! Ancienne discretisation verticale'

      h=7.
      snorm  = 0.
      DO l = 1, llm
         x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1)
         dsig(l) = 1.0 + 7.0 * SIN(x)**2
         snorm = snorm + dsig(l)
      ENDDO
      snorm = 1./snorm
      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

      ENDIF


      DO l=1,llm
        nivsigs(l) = FLOAT(l)
      ENDDO

      DO l=1,llmp1
        nivsig(l)= FLOAT(l)
      ENDDO

c ....  Calculs  de ap(l) et de bp(l)  ....
c .........................................
c
      pa        =  50 000.
      preff     = 101 325.
c                     .... Pascals ....

      bp(llmp1) =   0.

      DO l = 1, llm
c
ccc      ap(l) = 0.
ccc      bp(l) = sig(l)

       bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
       ap(l) = pa * ( sig(l) - bp(l) )

      ENDDO
       ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )

      PRINT *,' BP '
      PRINT *, bp
      PRINT *,' AP '
      PRINT *, ap

      DO l = 1, llm
       dpres(l) = bp(l) - bp(l+1)
       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
      ENDDO

      PRINT *,' PRESNIVS '
      PRINT *,presnivs
c-----------------------------------------------------------------------
c

c-----------------------------------------------------------------------
      RETURN
      END