!
! $Header: /home/cvsroot/LMDZ4/libf/phylmd/ajsec.F,v 1.1.1.1 2004/05/19 12:53:08 lmdzadmin Exp $
!
      SUBROUTINE ajsec(paprs, pplay, t, d_t)
      IMPLICIT none
c======================================================================
c inspiree de la base ecrite par:
c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
c Objet: ajustement sec (adaptation du GCM du LMD)
c
c Reecriture par S. Lebonnois en Novembre 2006
c pour prise en compte Cp variable dans atm de VENUS
c======================================================================
c Arguments:
c t-------input-R- Temperature
c
c d_t-----output-R-Incrementation de la temperature
c======================================================================
#include "dimensions.h"
#include "dimphy.h"
#include "YOMCST.h"
      REAL paprs(klon,klev+1), pplay(klon,klev)
      REAL t(klon,klev), d_t(klon,klev)
c
      INTEGER limbas, limhau ! les couches a ajuster
ccc      PARAMETER (limbas=klev-3, limhau=klev)
      PARAMETER (limbas=1, limhau=klev)
c
      REAL zt(klev),t_inter(klev),zh(klev)
      REAL zdp(klev),zdplay(klev)
      REAL alpha(klev),beta(klev)
      real somh,somhm1
      LOGICAL down
      INTEGER i, k, k1, k2

      real cpco2,cpstar(klev),cpdyn(klev)
      external cpco2
      save cpstar,cpdyn

      logical firstcall
      data    firstcall/.true./
      save    firstcall
c
c Initialisation:
c
      DO k = 1, klev
      DO i = 1, klon
         d_t(i,k) = 0.0
      ENDDO
      ENDDO

      DO i = 1, klon
c - - boucle sur les points de grille - - - - - -      

c------ initialisation des variables necessaires au calcul

        do k = 1, klev
	  zt(k)  = t(i,k)
	  zdp(k) = paprs(i,k)-paprs(i,k+1)

        if (firstcall) then
c CP considere: etre cohérent avec la dynamique ! (conservation...)
          cpdyn(k) = RCPD
c         cpdyn(k) = cpco2(t(i,k))
        endif

c enthalpie dans la couche k: c'est la grandeur qu'on mélange
	  zh(k)  = t(i,k)*zdp(k)*cpdyn(k)
	enddo
	
        do k = 2, klev
	  zdplay(k)  = pplay(i,k-1)-pplay(i,k)
	  t_inter(k) = (t(i,k-1)+t(i,k))/2.
        if (firstcall) then
          cpstar(k)  = cpco2(t_inter(k))
c----
c TEST si cp constant:
          cpstar(k)  = RCPD
          print*,"ATTENTION!! CP CONSTANT DANS AJSEC",cpstar
c----
c         print*,cpstar
        endif
	  beta(k)    = RD*zdplay(k)/(2.*cpstar(k)*paprs(i,k))
c       print*,beta(k),((zt(k-1)-zt(k))/(zt(k-1)+zt(k)))
	  alpha(k)   = (1.-beta(k))/(1.+beta(k))*(zdp(k)/zdp(k-1))
	enddo
c------------------------------------- correction des instabilites
c on part du bas:
          k2 = limbas

 8000     CONTINUE
c on cherche une zone de convection
            k2 = k2 + 1
	    somh   = 0.
	    somhm1 = 0.
            IF (k2 .GT. limhau) goto 8001
	    
c test instabilite initiale pour une zone
            IF (((zt(k2-1)-zt(k2))/(zt(k2-1)+zt(k2))).GT.beta(k2)) THEN
              k1 = k2 - 1
	      somh   = zh(k1)+zh(k2)
	      somhm1 = zh(k1)
              down = .FALSE.

 8020         CONTINUE
c calculs nouvelles temperatures entre k1 et k2 par recurrence
                if (down) then
                   zt(k1) = somh /( cpdyn(k1)*zdp(k1)*
     .      (1+somhm1/(cpdyn(k1+1)*zt(k1+1)*zdp(k1+1))*alpha(k1+1)) )
                   do k = k1+1, k2
		     zt(k) = zt(k-1)*(1-beta(k))/(1+beta(k))
		   enddo
		else
                   zt(k2) = somh /( cpdyn(k2)*zdp(k2)*
     .      (1+somhm1/(cpdyn(k2-1)*zt(k2-1)*zdp(k2-1))/alpha(k2)) )
                   do k = k2-1, k1, -1
		     zt(k) = zt(k+1)*(1+beta(k+1))/(1-beta(k+1))
		   enddo
		endif

                down = .FALSE.
                IF (k1 .ne. limbas) THEN
c test instabilite vers le bas
                 IF (((zt(k1-1)-zt(k1))/(zt(k1-1)+zt(k1))).GT.beta(k1))
     .                     down = .TRUE.
                ENDIF
                IF (down) THEN
                  k1 = k1 - 1
		  somhm1 = somh
		  somh   = somh + zh(k1)
                ELSE
                  IF ((k2 .EQ. limhau)) GOTO 8021
c test instabilite vers le haut
           IF (((zt(k2)-zt(k2+1))/(zt(k2)+zt(k2+1))).LE.beta(k2+1))
c                  (si vrai, stable)
     .                     GOTO 8021
                  k2 = k2 + 1
		  somhm1 = somh
		  somh   = somh + zh(k2)
                ENDIF
c on a etendu la zone, on recalcule:
              GOTO 8020

c sortie zone => on continue au-dessus pour chercher autre zone
 8021         CONTINUE
              k2 = k2 + 1
            ENDIF
          GOTO 8000

c sortie de la colonne:
 8001     CONTINUE

c ------------------------------------- tendances temperature

        DO k = 1,klev
         d_t(i,k) = zt(k) - t(i,k)
c        if(d_t(i,k).gt.0.) print*,zt(k),t(i,k)
        ENDDO

c - - fin boucle sur les points de grille - - - -      
      ENDDO

      firstcall  = .false.

      RETURN
      END