source: LMDZ4/trunk/libf/dyn3d/disvert.F @ 5085

!
! $Id: disvert.F 1403 2010-07-01 09:02:53Z fairhead $
!
      SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)

c    Auteur :  P. Le Van .
c
      IMPLICIT NONE

#include "dimensions.h"
#include "paramet.h"
#include "iniprint.h"
#include "logic.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 pa,preff
      REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
      REAL presnivs(llm)
c
c   declarations:
c   -------------
c
      REAL sig(llm+1),dsig(llm)
       real zzz(1:llm+1)
       real dzz(1:llm)
      real zk,zkm1,dzk1,dzk2,k0,k1
c
      INTEGER l
      REAL snorm,dsigmin
      REAL alpha,beta,gama,delta,deltaz,h
      INTEGER np,ierr
      REAL pi,x

      REAL 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,*) h           ! hauteur d'echelle 8.
      READ(99,*) deltaz      ! epaiseur de la premiere couche 0.04
      READ(99,*) beta        ! facteur d'acroissement en haut 1.3
      READ(99,*) k0          ! nombre de couches dans la transition surf
      READ(99,*) k1          ! nombre de couches dans la transition haute
      CLOSE(99)
      alpha=deltaz/(llm*h)
      write(lunout,*)'h,alpha,k0,k1,beta'

c     read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2

      alpha=deltaz/tanh(1./k0)*2.
      zkm1=0.
      sig(1)=1.
      do l=1,llm
        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h)
     + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
        zk=-h*log(sig(l+1))

        dzk1=alpha*tanh(l/k0)
        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
        write(lunout,*)l,sig(l+1),zk,zk-zkm1,dzk1,dzk2
        zkm1=zk
      enddo

      sig(llm+1)=0.

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

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

      WRITE(LUNOUT,*)'WARNING!!! Ancienne discretisation verticale'

      if (ok_strato) then
         if (llm==39) then
            dsigmin=0.3
         else if (llm==50) then
            dsigmin=1.
         else
            WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
            dsigmin=1.
         endif
         WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
      endif

      h=7.
      snorm  = 0.
      DO l = 1, llm
         x = 2.*asin(1.) * (REAL(l)-0.5) / REAL(llm+1)

         IF (ok_strato) THEN
           dsig(l) =(dsigmin + 7.0 * SIN(x)**2)
     &            *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2        
         ELSE
           dsig(l) = 1.0 + 7.0 * SIN(x)**2
         ENDIF

         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) = REAL(l)
      ENDDO

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

c
c    ....  Calculs  de ap(l) et de bp(l)  ....
c    .........................................
c
c
c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
c

      bp(llmp1) =   0.

      DO l = 1, llm
cc
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) )
c
      ENDDO

      bp(1)=1.
      ap(1)=0.

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

      write(lunout,*)' BP '
      write(lunout,*)  bp
      write(lunout,*)' AP '
      write(lunout,*)  ap

      write(lunout,*)
     .'Niveaux de pressions approximatifs aux centres des'
      write(lunout,*)'couches calcules pour une pression de surface =',
     .                 preff
      write(lunout,*)
     .     'et altitudes equivalentes pour une hauteur d echelle de'
      write(lunout,*)'8km'
      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 )
       write(lunout,*)'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))
      ENDDO

      write(lunout,*)' PRESNIVS '
      write(lunout,*)presnivs

      RETURN
      END