SUBROUTINE convadj(ngrid,nlay,ptimestep,
     S                   pplay,pplev,ppopsk,
     $                   pu,pv,ph,
     $                   pdufi,pdvfi,pdhfi,
     $                   pduadj,pdvadj,pdhadj)
      USE phys_const
      IMPLICIT NONE

c=======================================================================
c
c   ajustement convectif sec
c   on peut ajouter les tendances pdhfi au profil pdh avant l'ajustement
c
c=======================================================================

c-----------------------------------------------------------------------
c   declarations:
c   -------------

#include "dimensions.h"

c   arguments:
c   ----------

      INTEGER ngrid,nlay
      REAL ptimestep
      REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay)
      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay)
      REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay)
      REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay)

c   local:
c   ------

      INTEGER ig,i,l,l1,l2,jj
      INTEGER jcnt, jadrs(ngrid)

      REAL*8 sig(nlay+1),sdsig(nlay),dsig(nlay)
      REAL*8 zu(ngrid,nlay),zv(ngrid,nlay)
      REAL*8 zh(ngrid,nlay)
      REAL*8 zu2(ngrid,nlay),zv2(ngrid,nlay)
      REAL*8 zh2(ngrid,nlay)
      REAL*8 zhm,zsm,zum,zvm,zalpha

      LOGICAL vtest(ngrid),down

c
c-----------------------------------------------------------------------
c   initialisation:
c   ---------------
c
c
c-----------------------------------------------------------------------
c   detection des profils a modifier:
c   ---------------------------------
c   si le profil est a modifier
c   (i.e. ph(niv_sup) < ph(niv_inf) )
c   alors le tableau "vtest" est mis a .TRUE. ;
c   sinon, il reste a sa valeur initiale (.FALSE.)
c   cette operation est vectorisable
c   On en profite pour copier la valeur initiale de "ph"
c   dans le champ de travail "zh"


      DO 1010 l=1,nlay
         DO 1015 ig=1,ngrid
            zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep
            zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep
            zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep
1015     CONTINUE
1010  CONTINUE

      zu2(:,:)=zu(:,:)
      zv2(:,:)=zv(:,:)
      zh2(:,:)=zh(:,:)

      DO 1020 ig=1,ngrid
        vtest(ig)=.FALSE.
 1020 CONTINUE
c
      DO 1040 l=2,nlay
        DO 1060 ig=1,ngrid
CRAY      vtest(ig)=CVMGM(.TRUE. , vtest(ig),
CRAY .                      zh2(ig,l)-zh2(ig,l-1))
          IF(zh2(ig,l).LT.zh2(ig,l-1)) vtest(ig)=.TRUE.
 1060   CONTINUE
 1040 CONTINUE
c
CRAY  CALL WHENNE(ngrid, vtest, 1, 0, jadrs, jcnt)
      jcnt=0
      DO 1070 ig=1,ngrid
         IF(vtest(ig)) THEN
            jcnt=jcnt+1
            jadrs(jcnt)=ig
         ENDIF
1070  CONTINUE


c-----------------------------------------------------------------------
c  Ajustement des "jcnt" profils instables indices par "jadrs":
c  ------------------------------------------------------------
c
      DO 1080 jj = 1, jcnt
c
          i = jadrs(jj)
c
c   Calcul des niveaux sigma sur cette colonne
          DO l=1,nlay+1
            sig(l)=pplev(i,l)/pplev(i,1)
         ENDDO
         DO l=1,nlay
            dsig(l)=sig(l)-sig(l+1)
            sdsig(l)=ppopsk(i,l)*dsig(l)
         ENDDO
          l2 = 1
c
c      -- boucle de sondage vers le haut
c
cins$     Loop
 8000     CONTINUE
c
            l2 = l2 + 1
c
cins$       Exit
            IF (l2 .GT. nlay) Goto 8001
c
            IF (zh2(i, l2) .LT. zh2(i, l2-1)) THEN
c
c          -- l2 est le niveau le plus haut de la colonne instable
c
              l1 = l2 - 1
              l  = l1
              zsm = sdsig(l2)
              zhm = zh2(i, l2)
c
c          -- boucle de sondage vers le bas
c
cins$         Loop
 8020         CONTINUE
c
                zsm = zsm + sdsig(l)
                zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
c
c            -- doit on etendre la colonne vers le bas ?
c
c_EC (M1875) 20/6/87 : AND -> AND THEN
c
                down = .FALSE.
                IF (l1 .NE. 1) THEN    !--  and then
                  IF (zhm .LT. zh2(i, l1-1)) THEN
                    down = .TRUE.
                  END IF
                END IF
c
                IF (down) THEN
c
                  l1 = l1 - 1
                  l  = l1
c
                ELSE
c
c              -- peut on etendre la colonne vers le haut ?
c
cins$             Exit
                  IF (l2 .EQ. nlay) Goto 8021
c
cins$             Exit
                  IF (zh2(i, l2+1) .GE. zhm) Goto 8021
c
                  l2 = l2 + 1
                  l  = l2
c
                END IF
c
cins$         End Loop
              GO TO 8020
 8021         CONTINUE
c
c          -- nouveau profil : constant (valeur moyenne)
c
              zalpha=0.
              zum=0.
              zvm=0.
              DO 1100 l = l1, l2
                zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
                zh2(i, l) = zhm
                zum=zum+dsig(l)*zu(i,l)
                zvm=zvm+dsig(l)*zv(i,l)
 1100         CONTINUE
              zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
              zum=zum/(sig(l1)-sig(l2+1))
              zvm=zvm/(sig(l1)-sig(l2+1))
              IF(zalpha.GT.1.) THEN
                 PRINT*,'WARNING dans convadj zalpha=',zalpha
                 if(ig.eq.1) then
                    print*,'Au pole nord'
                 elseif (ig.eq.ngrid) then
                    print*,'Au pole sud'
                 else
                    print*,'Point i=',
     .              ig-((ig-1)/iim)*iim,'j=',(ig-1)/iim+1
                 endif
                 STOP
                 zalpha=1.
              ELSE
c                IF(zalpha.LT.0.) STOP'zalpha=0'
                 IF(zalpha.LT.1.e-5) zalpha=1.e-5
              ENDIF
              DO l=l1,l2
                 zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
                 zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
              ENDDO
 
              l2 = l2 + 1
c
            END IF
c
cins$     End Loop
          GO TO 8000
 8001     CONTINUE
c
 1080 CONTINUE
c
      DO 4000 l=1,nlay
        DO 4020 ig=1,ngrid
          pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep
          pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep
          pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep
c         pdhadj(ig,l)=0.
c         pduadj(ig,l)=0.
c         pdvadj(ig,l)=0.
 4020   CONTINUE
 4000 CONTINUE
c
      RETURN
      END