source: trunk/LMDZ.COMMON/libf/dyn3dpar/cpdet.F @ 847

! 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(ip1jmp1,llm, 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) :: ip1jmp1,llm
      real,intent(in) :: yt(ip1jmp1,llm)
      real,intent(out) :: yteta(ip1jmp1,llm)
      real,intent(in) :: ypk(ip1jmp1,llm)
! local variable:
      integer :: ij,l,ijb,ije
      
      !ijb=ij_begin
      !ije=ij_end
      ijb=1
      ije=ip1jmp1
     
      if (planet_type.eq."venus") then
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        do l=1,llm
          yteta(ijb:ije,l)=yt(ijb:ije,l)**nu_venus                      &
     &                     -nu_venus*t0_venus**nu_venus*                &
     &                          log(ypk(ijb:ije,l)/cpp)
          yteta(ijb:ije,l)=yteta(ijb:ije,l)**(1./nu_venus)
        enddo
!$OMP END DO 
      else
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        do l=1,llm
          yteta(ijb:ije,l)=yt(ijb:ije,l)*cpp/ypk(ijb:ije,l)
        enddo
!$OMP END DO 
      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(ip1jmp1,llm,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) :: ip1jmp1,llm
      real,intent(out) :: yt(ip1jmp1,llm)
      real,intent(in) :: yteta(ip1jmp1,llm)
      real,intent(in) :: ypk(ip1jmp1,llm)

! local variable:
      integer :: ij,l,ijb,ije

      !ijb=ij_begin
      !ije=ij_end
      ijb=1
      ije=ip1jmp1

      if (planet_type.eq."venus") then
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        do l=1,llm
          yt(ijb:ije,l)=yteta(ijb:ije,l)**nu_venus                      &
     &                  +nu_venus*t0_venus**nu_venus*                   &
     &                       log(ypk(ijb:ije,l)/cpp)
          yt(ijb:ije,l)=yt(ijb:ije,l)**(1./nu_venus)
        enddo
!$OMP END DO 
      else
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        do l=1,llm
          yt(ijb:ije,l)=yteta(ijb:ije,l)*ypk(ijb:ije,l)/cpp
        enddo
!$OMP END DO 
      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======================================================================