source: LMDZ5/trunk/libf/dyn3d_common/disvert_noterre.F @ 2225

! $Id: $
      SUBROUTINE disvert_noterre

c    Auteur :  F. Forget Y. Wanherdrick, P. Levan
c    Nouvelle version 100% Mars !!
c    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

      IMPLICIT NONE

#include "dimensions.h"
#include "paramet.h"
#include "comvert.h"
#include "comconst.h"
#include "logic.h"
#include "iniprint.h"
c
c=======================================================================
c    Discretisation verticale en coordonnée hybride (ou sigma)
c
c=======================================================================
c
c   declarations:
c   -------------
c
c
      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"

c
c-----------------------------------------------------------------------
c
! 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',
     s   iostat=ierr2)
         if(ierr2.ne.0) then
              close(99)
              open(99,file="z2sig.def",status='old',form='formatted',
     s        iostat=ierr4)
         endif

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

      IF(ierr2.eq.0) then

c        Lecture de esasig.def :
c        Systeme peu souple, mais qui respecte en theorie
c        La conservation de l'energie (conversion Energie potentielle
c        <-> 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
c
         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

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

      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

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

      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   .....  pa et preff sont lus  sur les fichiers start par dynetat0 .....
c-----------------------------------------------------------------------
c

      if (hybrid) then  ! use hybrid coordinates
         write(lunout,*) "*********************************"
         write(lunout,*) "Using hybrid vertical coordinates"
         write(lunout,*) 
c        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,*) 
c        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

c     Calcul au milieu des couches :
c     WARNING : le choix de placer le milieu des couches au niveau de
c     pression intermédiaire est arbitraire et pourrait etre modifié.
c     Le calcul du niveau pour la derniere couche 
c     (on met la meme distance (en log pression)  entre P(llm)
c     et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est
c     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

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


      RETURN
      END

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