source: LMDZ6/trunk/libf/dyn3d_common/disvert_noterre.F90 @ 5248

! $Id: $
SUBROUTINE disvert_noterre

  !    Auteur :  F. Forget Y. Wanherdrick, P. Levan
  !    Nouvelle version 100% Mars !!
  !    On l'utilise aussi pour Venus et Titan, legerment modifiee.

#ifdef CPP_IOIPSL
  use IOIPSL
#else
  ! if not using IOIPSL, we still need to use (a local version of) getin
  use ioipsl_getincom
#endif
  USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt, &
        nivsig,nivsigs,pa,preff,scaleheight
  USE comconst_mod, ONLY: kappa
  USE logic_mod, ONLY: hybrid

  IMPLICIT NONE

  include "dimensions.h"
  include "paramet.h"
  include "iniprint.h"
  !
  !=======================================================================
  !    Discretisation verticale en coordonnée hybride (ou sigma)
  !
  !=======================================================================
  !
  !   declarations:
  !   -------------
  !
  !
  INTEGER :: l,ll
  REAL :: snorm
  REAL :: alpha,beta,gama,delta,deltaz
  real :: quoi,quand
  REAL :: zsig(llm),sig(llm+1)
  INTEGER :: np,ierr
  integer :: ierr1,ierr2,ierr3,ierr4
  REAL :: x

  REAL :: SSUM
  EXTERNAL SSUM
  real :: newsig
  REAL :: dz0,dz1,nhaut,sig1,esig,csig,zz
  real :: tt,rr,gg, prevz
  real :: s(llm),dsig(llm)

  integer :: iz
  real :: z, ps,p
  character(len=*),parameter :: modname="disvert_noterre"

  !
  !-----------------------------------------------------------------------
  !
  ! Initializations:
  !  pi=2.*ASIN(1.) ! already done in iniconst

  hybrid=.true. ! default value for hybrid (ie: use hybrid coordinates)
  CALL getin('hybrid',hybrid)
  write(lunout,*) trim(modname),': hybrid=',hybrid

  ! Ouverture possible de fichiers typiquement E.T.

     open(99,file="esasig.def",status='old',form='formatted', &
           iostat=ierr2)
     if(ierr2.ne.0) then
          close(99)
          open(99,file="z2sig.def",status='old',form='formatted', &
                iostat=ierr4)
     endif

  !-----------------------------------------------------------------------
  !   cas 1 on lit les options dans esasig.def:
  !   ----------------------------------------

  IF(ierr2.eq.0) then

     ! Lecture de esasig.def :
     ! Systeme peu souple, mais qui respecte en theorie
     ! La conservation de l'energie (conversion Energie potentielle
     ! <-> energie cinetique, d'apres la note de Frederic Hourdin...

     write(lunout,*)'*****************************'
     write(lunout,*)'WARNING reading esasig.def'
     write(lunout,*)'*****************************'
     READ(99,*) scaleheight
     READ(99,*) dz0
     READ(99,*) dz1
     READ(99,*) nhaut
     CLOSE(99)

     dz0=dz0/scaleheight
     dz1=dz1/scaleheight

     sig1=(1.-dz1)/tanh(.5*(llm-1)/nhaut)

     esig=1.

     do l=1,20
        esig=-log((1./sig1-1.)*exp(-dz0)/esig)/(llm-1.)
     enddo
     csig=(1./sig1-1.)/(exp(esig)-1.)

     DO L = 2, llm
        zz=csig*(exp(esig*(l-1.))-1.)
        sig(l) =1./(1.+zz) &
              * tanh(.5*(llm+1-l)/nhaut)
     ENDDO
     sig(1)=1.
     sig(llm+1)=0.
     quoi      = 1. + 2.* kappa
     s( llm )  = 1.
     s(llm-1) = quoi
     IF( llm.gt.2 )  THEN
        DO  ll = 2, llm-1
           l         = llm+1 - ll
           quand     = sig(l+1)/ sig(l)
           s(l-1)    = quoi * (1.-quand) * s(l)  + quand * s(l+1)
        ENDDO
     END IF
  !
     snorm=(1.-.5*sig(2)+kappa*(1.-sig(2)))*s(1)+.5*sig(2)*s(2)
     DO l = 1, llm
        s(l)    = s(l)/ snorm
     ENDDO

  !-----------------------------------------------------------------------
  !   cas 2 on lit les options dans z2sig.def:
  !   ----------------------------------------

  ELSE IF(ierr4.eq.0) then
     write(lunout,*)'****************************'
     write(lunout,*)'Reading z2sig.def'
     write(lunout,*)'****************************'

     READ(99,*) scaleheight
     do l=1,llm
        read(99,*) zsig(l)
     end do
     CLOSE(99)

     sig(1) =1
     do l=2,llm
       sig(l) = 0.5 * ( exp(-zsig(l)/scaleheight) + &
             exp(-zsig(l-1)/scaleheight) )
     end do
     sig(llm+1) =0

  !-----------------------------------------------------------------------
  ELSE
     write(lunout,*) 'didn t you forget something ??? '
     write(lunout,*) 'We need file  z2sig.def ! (OR esasig.def)'
     stop
  ENDIF
  !-----------------------------------------------------------------------

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

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


  !-----------------------------------------------------------------------
  !    ....  Calculs  de ap(l) et de bp(l)  ....
  !    .........................................
  !
  !   .....  pa et preff sont lus  sur les fichiers start par dynetat0 .....
  !-----------------------------------------------------------------------
  !

  if (hybrid) then  ! use hybrid coordinates
     write(lunout,*) "*********************************"
     write(lunout,*) "Using hybrid vertical coordinates"
     write(lunout,*)
     ! Coordonnees hybrides avec mod
     DO l = 1, llm

     call sig_hybrid(sig(l),pa,preff,newsig)
        bp(l) = EXP( 1. - 1./(newsig**2)  )
        ap(l) = pa * (newsig - bp(l) )
     enddo
     bp(llmp1) = 0.
     ap(llmp1) = 0.
  else ! use sigma coordinates
     write(lunout,*) "********************************"
     write(lunout,*) "Using sigma vertical coordinates"
     write(lunout,*)
     ! Pour ne pas passer en coordonnees hybrides
     DO l = 1, llm
        ap(l) = 0.
        bp(l) = sig(l)
     ENDDO
     ap(llmp1) = 0.
  endif

  bp(llmp1) =   0.

  write(lunout,*) trim(modname),': BP '
  write(lunout,*)  bp
  write(lunout,*) trim(modname),': AP '
  write(lunout,*)  ap

  ! Calcul au milieu des couches :
  ! WARNING : le choix de placer le milieu des couches au niveau de
  ! pression intermédiaire est arbitraire et pourrait etre modifié.
  ! Le calcul du niveau pour la derniere couche
  ! (on met la meme distance (en log pression)  entre P(llm)
  ! et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est
  ! Specifique.  Ce choix est spécifié ici ET dans exner_milieu.F

  DO l = 1, llm-1
   aps(l) =  0.5 *( ap(l) +ap(l+1))
   bps(l) =  0.5 *( bp(l) +bp(l+1))
  ENDDO

  if (hybrid) then
     aps(llm) = aps(llm-1)**2 / aps(llm-2)
     bps(llm) = 0.5*(bp(llm) + bp(llm+1))
  else
     bps(llm) = bps(llm-1)**2 / bps(llm-2)
     aps(llm) = 0.
  end if

  write(lunout,*) trim(modname),': BPs '
  write(lunout,*)  bps
  write(lunout,*) trim(modname),': APs'
  write(lunout,*)  aps

  DO l = 1, llm
   presnivs(l) = aps(l)+bps(l)*preff
   pseudoalt(l) = -scaleheight*log(presnivs(l)/preff)
  ENDDO

  write(lunout,*)trim(modname),' : PRESNIVS'
  write(lunout,*)presnivs
  write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ', &
        'height of ',scaleheight,' km)'
  write(lunout,*)pseudoalt

  ! --------------------------------------------------
  ! This can be used to plot the vertical discretization
  ! (> xmgrace -nxy testhybrid.tab )
  ! --------------------------------------------------
  ! open (53,file='testhybrid.tab')
  ! scaleheight=15.5
  ! do iz=0,34
  !   z = -5 + min(iz,34-iz)
  ! approximation of scale height for Venus
  !   scaleheight = 15.5 - z/55.*10.
  !   ps = preff*exp(-z/scaleheight)
  !   zsig(1)= -scaleheight*log((aps(1) + bps(1)*ps)/preff)
  !   do l=2,llm
  ! approximation of scale height for Venus
  !      if (zsig(l-1).le.55.) then
  !         scaleheight = 15.5 - zsig(l-1)/55.*10.
  !      else
  !         scaleheight = 5.5 - (zsig(l-1)-55.)/35.*2.
  !      endif
  !      zsig(l)= zsig(l-1)-scaleheight*
  !    .    log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps))
  !   end do
  !   write(53,'(I3,50F10.5)') iz, zsig
  !  end do
  !  close(53)
  ! --------------------------------------------------


  RETURN
END SUBROUTINE disvert_noterre

! ************************************************************
subroutine sig_hybrid(sig,pa,preff,newsig)
  ! ----------------------------------------------
  ! Subroutine utilisee pour calculer des valeurs de sigma modifie
  ! pour conserver les coordonnees verticales decrites dans
  ! esasig.def/z2sig.def lors du passage en coordonnees hybrides
  ! F. Forget 2002
  ! Connaissant sig (niveaux "sigma" ou on veut mettre les couches)
  ! L'objectif est de calculer newsig telle que
  !   (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig
  ! Cela ne se résoud pas analytiquement:
  ! => on résoud par iterration bourrine
  ! ----------------------------------------------
  ! Information  : where exp(1-1./x**2) become << x
  !       x      exp(1-1./x**2) /x
  !       1           1
  !       0.68       0.5
  !       0.5        1.E-1
  !       0.391      1.E-2
  !       0.333      1.E-3
  !       0.295      1.E-4
  !       0.269      1.E-5
  !       0.248      1.E-6
  !    => on peut utiliser newsig = sig*preff/pa si sig*preff/pa < 0.25


  implicit none
  real :: x1, x2, sig,pa,preff, newsig, F
  integer :: j

  newsig = sig
  x1=0
  x2=1
  if (sig.ge.1) then
        newsig= sig
  else if (sig*preff/pa.ge.0.25) then
    DO J=1,9999  ! nombre d''iteration max
      F=((1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig
      ! write(0,*) J, ' newsig =', newsig, ' F= ', F
      if (F.gt.1) then
          X2 = newsig
          newsig=(X1+newsig)*0.5
      else
          X1 = newsig
          newsig=(X2+newsig)*0.5
      end if
      ! Test : on arete lorsque on approxime sig à moins de 0.01 m près
      ! (en pseudo altitude) :
      IF(abs(10.*log(F)).LT.1.E-5) goto 999
    END DO
   else   !    if (sig*preff/pa.le.0.25) then
         newsig= sig*preff/pa
   end if
 999   continue
   Return
END SUBROUTINE sig_hybrid