! ADAPTATION GCM POUR CP(T)
c======================================================================
c S. Lebonnois, 10/2010
c
c Cp doit être calculé par cpdet(t) pour être valable partout
c
c La fonction d'Exner reste pk = RCPD*(play/pref)**RKAPPA
c (RCPD=cpp, RKAPPA=kappa)
c
c On passe de T a teta (temperature potentielle) par t2tpot(t,teta,pk)
c On passe de teta a T par tpot2t(teta,t,pk)
c
c======================================================================

      SUBROUTINE ini_cpdet
      
      USE control_mod, ONLY: planet_type 
      IMPLICIT none
c======================================================================
c Initialisation de nu_venus et t0_venus
c======================================================================

! for cpp, nu_venus and t0_venus:
#include "comconst.h"

      if (planet_type.eq."venus") then
          nu_venus=0.35
          t0_venus=460.
      else
          nu_venus=0.
          t0_venus=0.
      endif

      return
      end

c======================================================================
c======================================================================

      FUNCTION cpdet(t)

      USE control_mod, ONLY: planet_type 
      IMPLICIT none

! for cpp, nu_venus and t0_venus:
#include "comconst.h"

      real cpdet,t

      if (planet_type.eq."venus") then
          cpdet = cpp*(t/t0_venus)**nu_venus
      else
          cpdet = cpp
      endif

      return
      end
      
c======================================================================
c======================================================================

      SUBROUTINE t2tpot(npoints, yt, yteta, ypk)
c======================================================================
c Arguments:
c
c yt   --------input-R- Temperature
c yteta-------output-R- Temperature potentielle
c ypk  --------input-R- Fonction d'Exner: RCPD*(pplay/pref)**RKAPPA
c
c======================================================================

      USE control_mod, ONLY: planet_type 
      IMPLICIT NONE
      
! for cpp, nu_venus and t0_venus:
#include "comconst.h"

      integer npoints
      REAL    yt(npoints), yteta(npoints), ypk(npoints)
      
      if (planet_type.eq."venus") then
          yteta = yt**nu_venus                                          &
     &            - nu_venus * t0_venus**nu_venus * log(ypk/cpp)
          yteta = yteta**(1./nu_venus)
      else
          yteta = yt * cpp/ypk
      endif

      return
      end

c======================================================================
c======================================================================

      SUBROUTINE t2tpot_p(nlon,nlev, yt, yteta, ypk)
! Parallel version of t2tpot
      USE parallel
      USE control_mod, only : planet_type
      IMPLICIT none

! for cpp, nu_venus and t0_venus:
#include "comconst.h"

      integer,intent(in) :: nlon,nlev
      real,intent(in) :: yt(nlon,nlev)
      real,intent(out) :: yteta(nlon,nlev)
      real,intent(in) :: ypk(nlon,nlev)
! local variable:
      integer :: l
      
      if (planet_type.eq."venus") then
        do l=1,nlev
          yteta(:,l)=yt(:,l)**nu_venus                                  &
     &                     -nu_venus*t0_venus**nu_venus*                &
     &                          log(ypk(:,l)/cpp)
          yteta(:,l)=yteta(:,l)**(1./nu_venus)
        enddo
      else
        do l=1,nlev
          yteta(:,l)=yt(:,l)*cpp/ypk(:,l)
        enddo
      endif ! of if (planet_type.eq."venus")

      END

c======================================================================

      SUBROUTINE tpot2t(npoints,yteta, yt, ypk)
c======================================================================
c Arguments:
c
c yteta--------input-R- Temperature potentielle
c yt   -------output-R- Temperature
c ypk  --------input-R- Fonction d'Exner: RCPD*(pplay/pref)**RKAPPA
c
c======================================================================

      USE control_mod, ONLY: planet_type 
      IMPLICIT NONE

! for cpp, nu_venus and t0_venus:
#include "comconst.h"

      integer npoints
      REAL    yt(npoints), yteta(npoints), ypk(npoints)
      
      if (planet_type.eq."venus") then
         yt = yteta**nu_venus                                           &
     &       + nu_venus * t0_venus**nu_venus * log(ypk/cpp)
         yt = yt**(1./nu_venus)
      else
          yt = yteta * ypk/cpp
      endif
  
      return
      end

c======================================================================
c======================================================================
      SUBROUTINE tpot2t_p(nlon,nlev,yteta,yt,ypk)
! Parallel version of tpot2t
      USE parallel
      USE control_mod, only : planet_type
      IMPLICIT none
! for cpp, nu_venus and t0_venus:
#include "comconst.h"

      integer,intent(in) :: nlon,nlev
      real,intent(out) :: yt(nlon,nlev)
      real,intent(in) :: yteta(nlon,nlev)
      real,intent(in) :: ypk(nlon,nlev)

! local variable:
      integer :: l

      if (planet_type.eq."venus") then
        do l=1,nlev
          yt(:,l)=yteta(:,l)**nu_venus                                  &
     &                  +nu_venus*t0_venus**nu_venus*                   &
     &                       log(ypk(:,l)/cpp)
          yt(:,l)=yt(:,l)**(1./nu_venus)
        enddo
      else
        do l=1,nlev
          yt(:,l)=yteta(:,l)*ypk(:,l)/cpp
        enddo
      endif ! of if (planet_type.eq."venus")
      END

c======================================================================
c======================================================================
c
c ATTENTION
c
c Si un jour on a besoin, il faudra coder les routines 
c    dt2dtpot / dtpto2dt 
c
c======================================================================
c======================================================================