module cpdet_mod

implicit none

! ADAPTATION OF GCM TO CP(T)
!======================================================================
! S. Lebonnois, 10/2010
!
! Cp must be computed using cpdet(t) to be valid
!
! The Exner function is still pk = RCPD*(play/pref)**RKAPPA
! (RCPD=cpp, RKAPPA=kappa)
!
! One goes from T to teta (potential temperature) using t2tpot(t,teta,pk)
! One goes from teta to T using tpot2t(teta,t,pk)
!
!======================================================================

contains

      SUBROUTINE ini_cpdet
      
      USE control_mod, ONLY: planet_type 
      IMPLICIT none
!======================================================================
! Initialization of nu_venus and t0_venus
!======================================================================

! 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 subroutine ini_cpdet

!======================================================================
!======================================================================

      FUNCTION cpdet(t)

      USE control_mod, ONLY: planet_type 
      IMPLICIT none

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

      real,intent(in) :: t
      real cpdet

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

      return
      end function cpdet
      
!======================================================================
!======================================================================

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

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

      integer,intent(in) :: npoints
      REAL,intent(in) :: yt(npoints), ypk(npoints)
      REAL,intent(out) :: yteta(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 subroutine t2tpot

!======================================================================
!======================================================================

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

      USE control_mod, ONLY: planet_type 
      IMPLICIT NONE

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

      integer,intent(in) :: npoints
      REAL,intent(in) :: yteta(npoints), ypk(npoints)
      REAL,intent(out) :: yt(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 subroutine tpot2t

!======================================================================
!======================================================================
! Routines pour les calculs paralleles
!======================================================================
#ifdef CPP_PARA
!======================================================================

      SUBROUTINE t2tpot_p(nlon,nlev, yt, yteta, ypk)
! Parallel version of t2tpot, for an arbitrary number of columns
      USE control_mod, only : planet_type
      USE parallel_lmdz, only : OMP_CHUNK
      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
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        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
!$OMP END DO
      else
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        do l=1,nlev
          yteta(:,l)=yt(:,l)*cpp/ypk(:,l)
        enddo
!$OMP END DO
      endif ! of if (planet_type.eq."venus")

      end subroutine t2tpot_p

!======================================================================
!======================================================================

      SUBROUTINE t2tpot_glo_p(yt, yteta, ypk)
! Parallel version of t2tpot, over the full dynamics (scalar) grid
! (more efficient than multiple calls to t2tpot_p() with slices of data)
      USE parallel_lmdz, only : jj_begin,jj_end,OMP_CHUNK
      USE control_mod, only : planet_type
      IMPLICIT none
! for iip1, jjp1 and llm
#include "dimensions.h"
#include "paramet.h"
! for cpp, nu_venus and t0_venus:
#include "comconst.h"

      real,intent(in) :: yt(iip1,jjp1,llm)
      real,intent(out) :: yteta(iip1,jjp1,llm)
      real,intent(in) :: ypk(iip1,jjp1,llm)
! local variable:
      integer :: j,l
      integer :: jjb,jje
      
      jjb=jj_begin
      jje=jj_end

      if (planet_type.eq."venus") then
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        do l=1,llm
          yteta(:,jjb:jje,l)=yt(:,jjb:jje,l)**nu_venus                  &
     &                     -nu_venus*t0_venus**nu_venus*                &
     &                          log(ypk(:,jjb:jje,l)/cpp)
          yteta(:,jjb:jje,l)=yteta(:,jjb:jje,l)**(1./nu_venus)
        enddo
!$OMP END DO
      else
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        do l=1,llm
          yteta(:,jjb:jje,l)=yt(:,jjb:jje,l)*cpp/ypk(:,jjb:jje,l)
        enddo
!$OMP END DO
      endif ! of if (planet_type.eq."venus")

      end subroutine t2tpot_glo_p

!======================================================================
!======================================================================

      SUBROUTINE tpot2t_p(nlon,nlev,yteta,yt,ypk)
! Parallel version of tpot2t, for an arbitrary number of columns
      USE control_mod, only : planet_type
      USE parallel_lmdz, only : OMP_CHUNK
      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
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        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
!$OMP END DO
      else
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        do l=1,nlev
          yt(:,l)=yteta(:,l)*ypk(:,l)/cpp
        enddo
!$OMP END DO
      endif ! of if (planet_type.eq."venus")
      end subroutine tpot2t_p

!======================================================================
!======================================================================

      SUBROUTINE tpot2t_glo_p(yteta,yt,ypk)
! Parallel version of tpot2t, over the full dynamics (scalar) grid
! (more efficient than multiple calls to tpot2t_p() with slices of data)
      USE parallel_lmdz, only : jj_begin,jj_end,OMP_CHUNK
      USE control_mod, only : planet_type
      IMPLICIT none
! for iip1, jjp1 and llm
#include "dimensions.h"
#include "paramet.h"
! for cpp, nu_venus and t0_venus:
#include "comconst.h"

      real,intent(out) :: yt(iip1,jjp1,llm)
      real,intent(in) :: yteta(iip1,jjp1,llm)
      real,intent(in) :: ypk(iip1,jjp1,llm)
! local variable:
      integer :: j,l
      integer :: jjb,jje
      
      jjb=jj_begin
      jje=jj_end

      if (planet_type.eq."venus") then
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        do l=1,llm
          yt(:,jjb:jje,l)=yteta(:,jjb:jje,l)**nu_venus                  &
     &                  +nu_venus*t0_venus**nu_venus*                   &
     &                       log(ypk(:,jjb:jje,l)/cpp)
          yt(:,jjb:jje,l)=yt(:,jjb:jje,l)**(1./nu_venus)
        enddo
!$OMP END DO
      else
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        do l=1,llm
          yt(:,jjb:jje,l)=yteta(:,jjb:jje,l)*ypk(:,jjb:jje,l)/cpp
        enddo
!$OMP END DO
      endif ! of if (planet_type.eq."venus")
      end subroutine tpot2t_glo_p

!======================================================================
#endif
!======================================================================
! Fin routines specifiques parallele
!======================================================================
!======================================================================
!
! ATTENTION
!
! Si un jour on a besoin, il faudra coder les routines 
!    dt2dtpot / dtpto2dt 
!
!======================================================================
!======================================================================
end module cpdet_mod