! $Id: disvert.F90 2040 2014-05-09 10:50:34Z evignon $

SUBROUTINE disvert()

#ifdef CPP_IOIPSL
  use ioipsl, only: getin
#else
  USE ioipsl_getincom, only: getin
#endif
  use new_unit_m, only: new_unit
  use assert_m, only: assert

  IMPLICIT NONE

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

!-------------------------------------------------------------------------------
! Purpose: Vertical distribution functions for LMDZ.
!          Triggered by the levels number llm.
!-------------------------------------------------------------------------------
! Read    in "comvert.h":
! pa                         !--- PURE PRESSURE COORDINATE FOR P<pa (in Pascals)
! preff                      !--- REFERENCE PRESSURE                 (101325 Pa)
! Written in "comvert.h":
! ap(llm+1), bp(llm+1)       !--- Ap, Bp HYBRID COEFFICIENTS AT INTERFACES
! aps(llm),  bps(llm)        !--- Ap, Bp HYBRID COEFFICIENTS AT MID-LAYERS
! dpres(llm)                 !--- PRESSURE DIFFERENCE FOR EACH LAYER
! presnivs(llm)              !--- PRESSURE AT EACH MID-LAYER
! scaleheight                !--- VERTICAL SCALE HEIGHT            (Earth: 8kms)
! nivsig(llm+1)              !--- SIGMA INDEX OF EACH LAYER INTERFACE
! nivsigs(llm)               !--- SIGMA INDEX OF EACH MID-LAYER
!-------------------------------------------------------------------------------
! Local variables:
  REAL sig(llm+1), dsig(llm)
  REAL sig0(llm+1), zz(llm+1)
  REAL zk, zkm1, dzk1, dzk2, z, k0, k1

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

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

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

  WRITE(lunout,*) TRIM(modname)//" starts"

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

  vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
  call getin('vert_sampling', vert_sampling)
  WRITE(lunout,*) TRIM(modname)//' 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("strato_custom") 
    ! custumize strato distribution with specific alpha & beta values and function
    ! depending on llm (experimental and temporary)!
    SELECT CASE (llm)
      CASE(55)
        alpha=0.45
        beta=4.0
      CASE(63)
        alpha=0.45
        beta=5.0
      CASE(71)
        alpha=3.05
        beta=65.
      CASE(79)
        alpha=3.20
        ! alpha=2.05 ! FLOTT 79 (PLANTE)
        beta=70.
    END SELECT
    ! Or used values provided by user in def file:
    CALL getin("strato_alpha",alpha)
    CALL getin("strato_beta",beta)
    
    ! Build geometrical distribution
    scaleheight=7.
    zz(1)=0.
    IF (llm==55.OR.llm==63) THEN
      DO l=1,llm
        z=zz(l)/scaleheight
        zz(l+1)=zz(l)+0.03+z*1.5*(1.-TANH(z-0.5))+alpha*(1.+TANH(z-1.5))     &
                            +5.0*EXP((l-llm)/beta)
      ENDDO
    ELSEIF (llm==71) THEN !.OR.llm==79) THEN      ! FLOTT 79 (PLANTE)
      DO l=1,llm
        z=zz(l)
        zz(l+1)=zz(l)+0.02+0.88*TANH(z/2.5)+alpha*(1.+TANH((z-beta)/15.))
      ENDDO
    ELSEIF (llm==79) THEN
      DO l=1,llm
        z=zz(l)
        zz(l+1)=zz(l)+0.02+0.80*TANH(z/3.8)+alpha*(1+TANH((z-beta)/17.))     &
                            +0.03*TANH(z/.25)
      ENDDO
    ENDIF ! of IF (llm==55.OR.llm==63) ...
    
    ! Build sigma distribution
    sig0=EXP(-zz(:)/scaleheight)
    sig0(llm+1)=0.
    sig=ridders(sig0)
    
    ! Compute ap() and bp()
    bp(1)=1.
    bp(2:llm)=EXP(1.-1./sig(2:llm)**2)
    bp(llm+1)=0.
    ap=pa*(sig-bp)

  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

CONTAINS

!-------------------------------------------------------------------------------
!
FUNCTION ridders(sig) RESULT(sg)
!
!-------------------------------------------------------------------------------
  IMPLICIT NONE
!-------------------------------------------------------------------------------
! Purpose: Search for s solving (Pa/Preff)*s+(1-Pa/Preff)*EXP(1-1./s**2)=sg
! Notes:   Uses Ridders' method, quite robust. Initial bracketing: 0<=sg<=1.
! Reference: Ridders, C. F. J. "A New Algorithm for Computing a Single Root of a
!       Real Continuous Function" IEEE Trans. Circuits Systems 26, 979-980, 1979
!-------------------------------------------------------------------------------
! Arguments:
  REAL, INTENT(IN)  :: sig(:)
  REAL              :: sg(SIZE(sig))
!-------------------------------------------------------------------------------
! Local variables:
  INTEGER :: it, ns, maxit
  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
!-------------------------------------------------------------------------------
  ns=SIZE(sig); maxit=9999
  c1=Pa/Preff; c2=1.-c1
  DO l=1,ns
    xx=HUGE(1.)
    x1=0.0; f1=distrib(x1,c1,c2,sig(l))
    x2=1.0; f2=distrib(x2,c1,c2,sig(l))
    DO it=1,maxit
      x3=0.5*(x1+x2); f3=distrib(x3,c1,c2,sig(l))
      s=SQRT(f3**2-f1*f2);                 IF(s==0.) EXIT
      x4=x3+(x3-x1)*(SIGN(1.,f1-f2)*f3/s); IF(ABS(10.*LOG(x4-xx))<=1E-5) EXIT
      xx=x4; f4=distrib(x4,c1,c2,sig(l));  IF(f4==0.) EXIT
      IF(SIGN(f3,f4)/=f3) THEN;      x1=x3; f1=f3; x2=xx; f2=f4
      ELSE IF(SIGN(f1,f4)/=f1) THEN; x2=xx; f2=f4
      ELSE IF(SIGN(f2,f4)/=f2) THEN; x1=xx; f1=f4
      ELSE; CALL abort_gcm("ridders",'Algorithm failed (which is odd...')
      END IF
      IF(ABS(10.*LOG(ABS(x2-x1)))<=1E-5) EXIT       !--- ERROR ON SIG <= 0.01m            
    END DO
    IF(it==maxit+1) WRITE(lunout,'(a,i3)')'WARNING in ridder: failed to converg&
     &e for level ',l
    sg(l)=xx
  END DO
  sg(1)=1.; sg(ns)=0.

END FUNCTION ridders

END SUBROUTINE disvert

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

FUNCTION distrib(x,c1,c2,x0) RESULT(res)
!
!-------------------------------------------------------------------------------
! Arguments:
  REAL, INTENT(IN) :: x, c1, c2, x0
  REAL             :: res
!-------------------------------------------------------------------------------
  res=c1*x+c2*EXP(1-1/(x**2))-x0

END FUNCTION distrib