MODULE convection
  IMPLICIT NONE

CONTAINS

  SUBROUTINE convadj(ngrid,nlay,ptimestep, &
       &                   pplay,pplev,ppopsk,   &
       &                  pu,pv,ph,              &
       &                  pdufi,pdvfi,pdhfi,     &
       &                  pduadj,pdvadj,pdhadj)
    USE phys_const
    
    !=======================================================================
    !
    !   ajustement convectif sec
    !   on peut ajouter les tendances pdhfi au profil pdh avant l'ajustement
    !
    !=======================================================================
    
    !-----------------------------------------------------------------------
    !   declarations:
    !   -------------
    
#include "dimensions.h"
    
    !   arguments:
    !   ----------
    
    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)
    
    !   local:
    !   ------
    
    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
    
    !
    !-----------------------------------------------------------------------
    !   initialisation:
    !   ---------------
    !
    !
    !-----------------------------------------------------------------------
    !   detection des profils a modifier:
    !   ---------------------------------
    !   si le profil est a modifier
    !   (i.e. ph(niv_sup) < ph(niv_inf) )
    !   alors le tableau "vtest" est mis a .TRUE. ;
    !   sinon, il reste a sa valeur initiale (.FALSE.)
    !   cette operation est vectorisable
    !   On en profite pour copier la valeur initiale de "ph"
    !   dans le champ de travail "zh"
    
    DO l=1,nlay
       DO 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
       END DO
    END DO
    
    zu2(:,:)=zu(:,:)
    zv2(:,:)=zv(:,:)
    zh2(:,:)=zh(:,:)
    
    DO ig=1,ngrid
       vtest(ig)=.FALSE.
    END DO
    !
    DO l=2,nlay
       DO 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.
       END DO
    END DO
    !
    !CRAY  CALL WHENNE(ngrid, vtest, 1, 0, jadrs, jcnt)
    jcnt=0
    DO ig=1,ngrid
       IF(vtest(ig)) THEN
          jcnt=jcnt+1
          jadrs(jcnt)=ig
       ENDIF
    END DO
    
    !-----------------------------------------------------------------------
    !  Ajustement des "jcnt" profils instables indices par "jadrs":
    !  ------------------------------------------------------------
    !
    DO jj = 1, jcnt
       !
       i = jadrs(jj)
       !
       !   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
       !
       !      -- boucle de sondage vers le haut
       !
       DO WHILE(.TRUE.)
          ! 8000     CONTINUE
          !
          l2 = l2 + 1
          !
          IF (l2 .GT. nlay) EXIT ! Goto 8001
          !
          IF (zh2(i, l2) .LT. zh2(i, l2-1)) THEN
             !
             !          -- l2 est le niveau le plus haut de la colonne instable
             !
             l1 = l2 - 1
             l  = l1
             zsm = sdsig(l2)
             zhm = zh2(i, l2)
             !
             !          -- boucle de sondage vers le bas
             !
             DO WHILE(.TRUE.)                   
                ! 8020         CONTINUE
                zsm = zsm + sdsig(l)
                zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
                !
                !            -- doit on etendre la colonne vers le bas ?
                !
                !_EC (M1875) 20/6/87 : AND -> AND THEN
                !
                down = .FALSE.
                IF (l1 .NE. 1) THEN    !--  and then
                   IF (zhm .LT. zh2(i, l1-1)) THEN
                      down = .TRUE.
                   END IF
                END IF
                !
                IF (down) THEN
                   l1 = l1 - 1
                   l  = l1
                ELSE
                   !              -- peut on etendre la colonne vers le haut ?
                   IF (l2 .EQ. nlay) EXIT !Goto 8021
                   IF (zh2(i, l2+1) .GE. zhm) EXIT !Goto 8021
                   l2 = l2 + 1
                   l  = l2
                END IF
                
                ! GO TO 8020
             END DO
             ! 8021         CONTINUE
             !
             !          -- nouveau profil : constant (valeur moyenne)
             !
             zalpha=0.
             zum=0.
             zvm=0.
             DO 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)
             END DO
             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
                !                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
             !
          END IF
          !
          !          GO TO 8000
       END DO
       ! 8001      CONTINUE
       !
       !
       DO l=1,nlay
          DO 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
             !         pdhadj(ig,l)=0.
             !         pduadj(ig,l)=0.
             !         pdvadj(ig,l)=0.
          END DO
       END DO
       !
    END DO
    
  END SUBROUTINE convadj

END MODULE convection