subroutine convadj(ngrid,nlay,nq,ptimestep,
     &                   pplay,pplev,ppopsk,lmax_th,
     &                   pu,pv,ph,pq,
     &                   pdufi,pdvfi,pdhfi,pdqfi,
     &                   pduadj,pdvadj,pdhadj,
     &                   pdqadj)

      implicit none

!==================================================================
!     
!     Purpose
!     -------
!     Calculates dry convective adjustment. If one tracer is CO2,
!     we take into account the molecular mass variation (e.g. when
!     CO2 condenses) to trigger convection (F. Forget 01/2005)
!
!     Authors
!     -------
!     Original author unknown.
!     Modif. 2005 by F. Forget.
!     Modif. 2010 by R. Wordsworth
!     Modif. 2011 by A. Colaitis
!
!==================================================================

!     ------------
!     Declarations
!     ------------

#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"
#include "tracer.h"


!     Arguments
!     ---------

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

!     Tracers
      integer,intent(in) :: nq
      real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
      real,intent(out) :: pdqadj(ngrid,nlay,nq)


!     Local
!     -----

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

      REAL sig(nlayermx+1),sdsig(nlayermx),dsig(nlayermx)
      REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
      REAL zh(ngridmx,nlayermx)
      REAL zu2(ngridmx,nlayermx),zv2(ngridmx,nlayermx)
      REAL zh2(ngridmx,nlayermx), zhc(ngridmx,nlayermx)
      REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc

!     Tracers
      INTEGER iq,ico2
      save ico2
      REAL zq(ngridmx,nlayermx,nqmx), zq2(ngridmx,nlayermx,nqmx)
      REAL zqm(nqmx),zqco2m
      real m_co2, m_noco2, A , B
      save A, B

      real mtot1, mtot2 , mm1, mm2
       integer l1ref, l2ref
      LOGICAL vtest(ngridmx),down,firstcall
      save firstcall
      data firstcall/.true./

!     for conservation test
      real masse,cadjncons

!     --------------
!     Initialisation
!     --------------

      IF (firstcall) THEN 
        IF(ngrid.NE.ngridmx) THEN
           PRINT*
           PRINT*,'STOP in convadj'
           PRINT*,'ngrid    =',ngrid
           PRINT*,'ngridmx  =',ngridmx
        ENDIF
        ico2=0
        if (tracer) then
!     Prepare Special treatment if one of the tracers is CO2 gas
           do iq=1,nqmx
             if (noms(iq).eq."co2") then
                ico2=iq
                m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
                m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
!               Compute A and B coefficient use to compute
!               mean molecular mass Mair defined by
!               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
!               1/Mair = A*q(ico2) + B
                A =(1/m_co2 - 1/m_noco2)
                B=1/m_noco2
             end if
           enddo
        endif
        firstcall=.false.
      ENDIF ! of IF (firstcall)

      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
         ENDDO
      ENDDO

      if(tracer) then      
        DO iq =1, nq
         DO l=1,nlay
           DO ig=1,ngrid
              zq(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
           ENDDO
         ENDDO
        ENDDO
      end if

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

!     -----------------------------
!     Detection of unstable columns
!     -----------------------------
!     If ph(above) < ph(below) we set vtest=.true.

      DO ig=1,ngrid
        vtest(ig)=.false.
      ENDDO

      if (ico2.ne.0) then
!     Special case if one of the tracers is CO2 gas
         DO l=1,nlay
           DO ig=1,ngrid
             zhc(ig,l) = zh2(ig,l)*(A*zq2(ig,l,ico2)+B)
           ENDDO
         ENDDO
       else
          zhc(:,:)=zh2(:,:)
       end if

!     Find out which grid points are convectively unstable
      DO l=2,nlay
        DO ig=1,ngrid
          IF((zhc(ig,l).LT.zhc(ig,l-1))                                 $
     $  .AND. (l .GT. lmax_th(ig))) vtest(ig)=.true.
       ENDDO
      ENDDO

!     Make a list of them
      jcnt=0
      DO ig=1,ngrid
         IF(vtest(ig)) THEN
            jcnt=jcnt+1
            jadrs(jcnt)=ig
         ENDIF
      ENDDO


!     ---------------------------------------------------------------
!     Adjustment of the "jcnt" unstable profiles indicated by "jadrs"
!     ---------------------------------------------------------------

      DO jj = 1, jcnt   ! loop on every convective grid point

          i = jadrs(jj)
 
!     Calculate sigma in this column
          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

!     Test loop upwards on l2

      DO
       l2 = l2 + 1
       IF (l2 .GT. nlay) EXIT
       IF ((zhc(i, l2) .LT. zhc(i, l2-1)).AND.(l2 .GT. lmax_th(i))) THEN

!     l2 is the highest level of the unstable column
 
              l1 = l2 - 1
              l  = l1
              zsm = sdsig(l2)
              zdsm = dsig(l2)
              zhm = zh2(i, l2)
              if(ico2.ne.0) zqco2m = zq2(i,l2,ico2)

!     Test loop downwards

              DO
                zsm = zsm + sdsig(l)
                zdsm = zdsm + dsig(l)
                zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
                if(ico2.ne.0) then
                  zqco2m = 
     &            zqco2m + dsig(l) * (zq2(i,l,ico2) - zqco2m) / zdsm
                  zhmc = zhm*(A*zqco2m+B)
                else 
                  zhmc = zhm
                end if
 
!     do we have to extend the column downwards?
 
            down = .false.
            IF (l1 .ne. 1) then    !--  and then
              IF ((zhmc .lt. zhc(i, l1-1)).and.(l1.gt.lmax_th(i))) then
                down = .true.
              END IF
            END IF
 
                ! this could be a problem...

                if (down) then
 
                  l1 = l1 - 1
                  l  = l1
 
                else
 
!     can we extend the column upwards?
 
                  if (l2 .eq. nlay) exit
 
                  if (zhc(i, l2+1) .ge. zhmc) exit

                  l2 = l2 + 1
                  l  = l2

                end if ! of if (down)

              enddo ! of do ! Test loop forward

!     New constant profile (average value)


              zalpha=0.
              zum=0.
              zvm=0.
              do iq=1,nq
                zqm(iq) = 0.
              end do
              DO l = l1, l2
                if(ico2.ne.0) then
                  zalpha=zalpha+
     &            ABS(zhc(i,l)/(A+B*zqco2m) -zhm)*dsig(l)
                else
                  zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
                endif
                zh2(i, l) = zhm
!     modifs by RDW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                zum=zum+dsig(l)*zu2(i,l)
                zvm=zvm+dsig(l)*zv2(i,l)
!                zum=zum+dsig(l)*zu(i,l)
!                zvm=zvm+dsig(l)*zv(i,l)
                do iq=1,nq
                   zqm(iq) = zqm(iq)+dsig(l)*zq2(i,l,iq)
!                   zqm(iq) = zqm(iq)+dsig(l)*zq(i,l,iq)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!     to conserve tracers/ KE, we must calculate zum, zvm and zqm using 
!     the up-to-date column values. If we do not do this, there are cases 
!     where convection stops at one level and starts at the next where we
!     can break conservation of stuff (particularly tracers) significantly.

                end do
              ENDDO ! of DO l = l1, l2
              zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
              zum=zum/(sig(l1)-sig(l2+1))
              zvm=zvm/(sig(l1)-sig(l2+1))
              do iq=1,nq
                 zqm(iq) = zqm(iq)/(sig(l1)-sig(l2+1))
              end do

              IF(zalpha.GT.1.) THEN
                 zalpha=1.
              ELSE
!                IF(zalpha.LT.0.) STOP
                 IF(zalpha.LT.1.e-4) zalpha=1.e-4
              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))
                 do iq=1,nq
!                  zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq)) 
                   zq2(i,l,iq)=zqm(iq)
                 end do
              ENDDO
              if (ico2.ne.0) then
                DO l=l1, l2
                  zhc(i,l) = zh2(i,l)*(A*zq2(i,l,ico2)+B)
                ENDDO
              end if


              l2 = l2 + 1

            END IF   ! End of l1 to l2 instability treatment
                     ! We now continue to test from l2 upwards

          ENDDO   ! End of upwards loop on l2


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     check conservation
         cadjncons=0.0
         if(water)then
          iq    = igcm_h2o_vap
          do l = 1, nlay
            masse = (pplev(i,l) - pplev(i,l+1))/g
            cadjncons = cadjncons + 
     &           masse*(zq2(i,l,iq)-zq(i,l,iq))/ptimestep 
          end do
         endif

         if(cadjncons.lt.-1.e-6)then
            print*,'convadj has just crashed...'
            print*,'i  = ',i
            print*,'l1 = ',l1
            print*,'l2 = ',l2
            print*,'cadjncons        = ',cadjncons
         do l = 1, nlay
            print*,'dsig         = ',dsig(l)
         end do         
         do l = 1, nlay
            print*,'dsig         = ',dsig(l)
         end do
         do l = 1, nlay
            print*,'sig         = ',sig(l)
         end do
         do l = 1, nlay
            print*,'pplay(ig,:)         = ',pplay(i,l)
         end do
         do l = 1, nlay+1
            print*,'pplev(ig,:)         = ',pplev(i,l)
         end do
         do l = 1, nlay
            print*,'ph(ig,:)         = ',ph(i,l)
         end do
         do l = 1, nlay
            print*,'zh(ig,:)         = ',zh(i,l)
         end do
         do l = 1, nlay
            print*,'zh2(ig,:)        = ',zh2(i,l)
         end do
         do l = 1, nlay
            print*,'zq(ig,:,vap)     = ',zq(i,l,igcm_h2o_vap)
         end do
         do l = 1, nlay
            print*,'zq2(ig,:,vap)    = ',zq2(i,l,igcm_h2o_vap)
         end do
            print*,'zqm(vap)         = ',zqm(igcm_h2o_vap)
            print*,'jadrs=',jadrs

            call abort
         endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      ENDDO ! of DO jj = 1, jcnt   ! loop on every convective grid point


      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
        ENDDO
      ENDDO

      if(tracer) then 
        do iq=1, nq
          do  l=1,nlay
            do ig=1, ngrid
              pdqadj(ig,l,iq)=(zq2(ig,l,iq)-zq(ig,l,iq))/ptimestep 
            end do 
          end do 
        end do 
      end if


!     output
!      if (ngrid.eq.1) then
!         ig=1
!         iq =1 
!         write(*,*)'**convadj: l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)'  
!         do l=nlay,1,-1
!           write(*,*) l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)
!         end do
!      end if


      return
      end