! $Id: disvert.F90 1959 2014-02-07 09:37:38Z idelkadi $

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

  ! Auteur : P. Le Van

  use new_unit_m, only: new_unit
  use ioipsl, only: getin
  use assert_m, only: assert

  IMPLICIT NONE

  include "dimensions.h"
  include "paramet.h"
  include "iniprint.h"
  include "logic.h"

  ! s = sigma ** kappa : coordonnee verticale
  ! dsig(l) : epaisseur de la couche l ds la coord. s
  ! sig(l) : sigma a l'interface des couches l et l-1
  ! ds(l) : distance entre les couches l et l-1 en coord.s

  real,intent(in) :: pa, preff
  real,intent(out) :: ap(llmp1) ! in Pa
  real, intent(out):: bp(llmp1)
  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
  real,intent(out) :: presnivs(llm)
  real,intent(out) :: scaleheight

  REAL sig(llm+1), dsig(llm)
  real zk, zkm1, dzk1, dzk2, k0, k1

  INTEGER l, unit
  REAL dsigmin
  REAL alpha, beta, deltaz
  REAL x
  character(len=*),parameter :: modname="disvert"

  character(len=6):: vert_sampling
  ! (allowed values are "param", "tropo", "strato" and "read")

  !-----------------------------------------------------------------------

  print *, "Call sequence information: disvert"

  ! default scaleheight is 8km for earth
  scaleheight=8.

  vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
  call getin('vert_sampling', vert_sampling)
  print *, 'vert_sampling = ' // vert_sampling
  if (llm==39 .and. vert_sampling=="strato") then
     dsigmin=0.3 ! Vieille option par défaut pour CMIP5
  else
     dsigmin=1.
  endif
  call getin('dsigmin', dsigmin)
  WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin


  select case (vert_sampling)
  case ("param")
     ! On lit les options dans sigma.def:
     OPEN(99, file='sigma.def', status='old', form='formatted')
     READ(99, *) scaleheight ! 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*scaleheight)
     write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
                               scaleheight, alpha, k0, k1, beta

     alpha=deltaz/tanh(1./k0)*2.
     zkm1=0.
     sig(1)=1.
     do l=1, llm
        sig(l+1)=(cosh(l/k0))**(-alpha*k0/scaleheight) &
             *exp(-alpha/scaleheight*tanh((llm-k1)/k0) &
                  *beta**(l-(llm-k1))/log(beta))
        zk=-scaleheight*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.

     bp(: llm) = EXP(1. - 1. / sig(: llm)**2)
     bp(llmp1) = 0.

     ap = pa * (sig - bp)
  case("sigma")
     DO l = 1, llm
        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
        dsig(l) = dsigmin + 7.0 * SIN(x)**2
     ENDDO
     dsig = dsig / sum(dsig)
     sig(llm+1) = 0.
     DO l = llm, 1, -1
        sig(l) = sig(l+1) + dsig(l)
     ENDDO

     bp(1)=1.
     bp(2: llm) = sig(2:llm)
     bp(llmp1) = 0.
     ap(:)=0.
  case("tropo")
     DO l = 1, llm
        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
        dsig(l) = dsigmin + 7.0 * SIN(x)**2
     ENDDO
     dsig = dsig / sum(dsig)
     sig(llm+1) = 0.
     DO l = llm, 1, -1
        sig(l) = sig(l+1) + dsig(l)
     ENDDO

     bp(1)=1.
     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
     bp(llmp1) = 0.

     ap(1)=0.
     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
  case("strato")
     DO l = 1, llm
        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
        dsig(l) =(dsigmin + 7. * SIN(x)**2) &
             *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
     ENDDO
     dsig = dsig / sum(dsig)
     sig(llm+1) = 0.
     DO l = llm, 1, -1
        sig(l) = sig(l+1) + dsig(l)
     ENDDO

     bp(1)=1.
     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
     bp(llmp1) = 0.

     ap(1)=0.
     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
  case("read")
     ! Read "ap" and "bp". First line is skipped (title line). "ap"
     ! should be in Pa. First couple of values should correspond to
     ! the surface, that is : "bp" should be in descending order.
     call new_unit(unit)
     open(unit, file="hybrid.txt", status="old", action="read", &
          position="rewind")
     read(unit, fmt=*) ! skip title line
     do l = 1, llm + 1
        read(unit, fmt=*) ap(l), bp(l)
     end do
     close(unit)
     call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
          bp(llm + 1) == 0., "disvert: bad ap or bp values")
  case default
     call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
  END select

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

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

  write(lunout, *)  trim(modname),': BP '
  write(lunout, *) bp
  write(lunout, *)  trim(modname),': 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, *) scaleheight,' km'
  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))*scaleheight &
          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
  ENDDO

  write(lunout, *) trim(modname),': PRESNIVS '
  write(lunout, *) presnivs

END SUBROUTINE disvert