Index: src/convmix_gfs.f90
===================================================================
--- src/convmix_gfs.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/convmix_gfs.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
@@ -0,0 +1,303 @@
+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
+!                                                                     *
+! This file is part of FLEXPART.                                      *
+!                                                                     *
+! FLEXPART is free software: you can redistribute it and/or modify    *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or   *
+! (at your option) any later version.                                 *
+!                                                                     *
+! FLEXPART is distributed in the hope that it will be useful,         *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
+! GNU General Public License for more details.                        *
+!                                                                     *
+! You should have received a copy of the GNU General Public License   *
+! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
+!**********************************************************************
+
+subroutine convmix(itime)
+  !                     i
+  !**************************************************************
+  !handles all the calculations related to convective mixing
+  !Petra Seibert, Bernd C. Krueger, Feb 2001
+  !nested grids included, Bernd C. Krueger, May 2001
+  !
+  !Changes by Caroline Forster, April 2004 - February 2005:
+  !  convmix called every lsynctime seconds
+  !CHANGES by A. Stohl:
+  !  various run-time optimizations - February 2005
+  !CHANGES by C. Forster, November 2005, NCEP GFS version
+  !      in the ECMWF version convection is calculated on the
+  !      original eta-levels
+  !      in the GFS version convection is calculated on the
+  !      FLEXPART levels
+  !**************************************************************
+
+  use par_mod
+  use com_mod
+  use conv_mod
+
+  implicit none
+
+  integer :: igr,igrold, ipart, itime, ix, j, inest
+  integer :: ipconv
+  integer :: jy, kpart, ktop, ngrid,kz
+  integer :: igrid(maxpart), ipoint(maxpart), igridn(maxpart,maxnests)
+  ! itime [s]                 current time
+  ! igrid(maxpart)            horizontal grid position of each particle
+  ! igridn(maxpart,maxnests)  dto. for nested grids
+  ! ipoint(maxpart)           pointer to access particles according to grid position
+
+  logical :: lconv
+  real :: x, y, xtn,ytn, ztold, delt
+  real :: dt1,dt2,dtt
+  integer :: mind1,mind2
+  ! dt1,dt2,dtt,mind1,mind2       variables used for time interpolation
+  integer :: itage,nage
+
+  !monitoring variables
+  !real sumconv,sumall
+
+
+  ! Calculate auxiliary variables for time interpolation
+  !*****************************************************
+
+  dt1=real(itime-memtime(1))
+  dt2=real(memtime(2)-itime)
+  dtt=1./(dt1+dt2)
+  mind1=memind(1)
+  mind2=memind(2)
+  delt=real(abs(lsynctime))
+
+
+  lconv = .false.
+
+  ! if no particles are present return after initialization
+  !********************************************************
+
+  if (numpart.le.0) return
+
+  ! Assign igrid and igridn, which are pseudo grid numbers indicating particles
+  ! that are outside the part of the grid under consideration
+  ! (e.g. particles near the poles or particles in other nests).
+  ! Do this for all nests but use only the innermost nest; for all others
+  ! igrid shall be -1
+  ! Also, initialize index vector ipoint
+  !************************************************************************
+
+  do ipart=1,numpart
+    igrid(ipart)=-1
+    do j=numbnests,1,-1
+      igridn(ipart,j)=-1
+    end do
+    ipoint(ipart)=ipart
+  ! do not consider particles that are (yet) not part of simulation
+    if (itra1(ipart).ne.itime) goto 20
+    x = xtra1(ipart)
+    y = ytra1(ipart)
+
+  ! Determine which nesting level to be used
+  !**********************************************************
+
+    ngrid=0
+    do j=numbnests,1,-1
+      if ( x.gt.xln(j) .and. x.lt.xrn(j) .and. &
+           y.gt.yln(j) .and. y.lt.yrn(j) ) then
+        ngrid=j
+        goto 23
+      endif
+    end do
+ 23   continue
+
+  ! Determine nested grid coordinates
+  !**********************************
+
+    if (ngrid.gt.0) then
+  ! nested grids
+      xtn=(x-xln(ngrid))*xresoln(ngrid)
+      ytn=(y-yln(ngrid))*yresoln(ngrid)
+      ix=nint(xtn)
+      jy=nint(ytn)
+      igridn(ipart,ngrid) = 1 + jy*nxn(ngrid) + ix
+    else if(ngrid.eq.0) then
+  ! mother grid
+      ix=nint(x)
+      jy=nint(y)
+      igrid(ipart) = 1 + jy*nx + ix
+    endif
+
+ 20 continue
+  end do
+
+  !sumall = 0.
+  !sumconv = 0.
+
+  !*****************************************************************************
+  ! 1. Now, do everything for the mother domain and, later, for all of the nested domains
+  ! While all particles have to be considered for redistribution, the Emanuel convection
+  ! scheme only needs to be called once for every grid column where particles are present.
+  ! Therefore, particles are sorted according to their grid position. Whenever a new grid
+  ! cell is encountered by looping through the sorted particles, the convection scheme is called.
+  !*****************************************************************************
+
+  ! sort particles according to horizontal position and calculate index vector IPOINT
+
+  call sort2(numpart,igrid,ipoint)
+
+  ! Now visit all grid columns where particles are present
+  ! by going through the sorted particles
+
+  igrold = -1
+  do kpart=1,numpart
+    igr = igrid(kpart)
+    if (igr .eq. -1) goto 50
+    ipart = ipoint(kpart)
+  !  sumall = sumall + 1
+    if (igr .ne. igrold) then
+  ! we are in a new grid column
+      jy = (igr-1)/nx
+      ix = igr - jy*nx - 1
+
+  ! Interpolate all meteorological data needed for the convection scheme
+      psconv=(ps(ix,jy,1,mind1)*dt2+ps(ix,jy,1,mind2)*dt1)*dtt
+      tt2conv=(tt2(ix,jy,1,mind1)*dt2+tt2(ix,jy,1,mind2)*dt1)*dtt
+      td2conv=(td2(ix,jy,1,mind1)*dt2+td2(ix,jy,1,mind2)*dt1)*dtt
+!!$      do kz=1,nconvlev+1    !old
+      do kz=1,nuvz-1           !bugfix
+        pconv(kz)=(pplev(ix,jy,kz,mind1)*dt2+ &
+             pplev(ix,jy,kz,mind2)*dt1)*dtt
+        tconv(kz)=(tt(ix,jy,kz,mind1)*dt2+ &
+             tt(ix,jy,kz,mind2)*dt1)*dtt
+        qconv(kz)=(qv(ix,jy,kz,mind1)*dt2+ &
+             qv(ix,jy,kz,mind2)*dt1)*dtt
+      end do
+
+  ! Calculate translocation matrix
+      call calcmatrix(lconv,delt,cbaseflux(ix,jy))
+      igrold = igr
+      ktop = 0
+    endif
+
+  ! treat particle only if column has convection
+    if (lconv .eqv. .true.) then
+  ! assign new vertical position to particle
+
+      ztold=ztra1(ipart)
+      call redist(ipart,ktop,ipconv)
+  !    if (ipconv.le.0) sumconv = sumconv+1
+
+  ! Calculate the gross fluxes across layer interfaces
+  !***************************************************
+
+      if (iflux.eq.1) then
+        itage=abs(itra1(ipart)-itramem(ipart))
+        do nage=1,nageclass
+          if (itage.lt.lage(nage)) goto 37
+        end do
+ 37     continue
+
+        if (nage.le.nageclass) &
+             call calcfluxes(nage,ipart,real(xtra1(ipart)), &
+             real(ytra1(ipart)),ztold)
+      endif
+
+    endif   !(lconv .eqv. .true)
+50  continue
+  end do
+
+
+  !*****************************************************************************
+  ! 2. Nested domains
+  !*****************************************************************************
+
+  ! sort particles according to horizontal position and calculate index vector IPOINT
+
+  do inest=1,numbnests
+    do ipart=1,numpart
+      ipoint(ipart)=ipart
+      igrid(ipart) = igridn(ipart,inest)
+    enddo
+    call sort2(numpart,igrid,ipoint)
+
+  ! Now visit all grid columns where particles are present
+  ! by going through the sorted particles
+
+    igrold = -1
+    do kpart=1,numpart
+      igr = igrid(kpart)
+      if (igr .eq. -1) goto 60
+      ipart = ipoint(kpart)
+  !    sumall = sumall + 1
+      if (igr .ne. igrold) then
+  ! we are in a new grid column
+        jy = (igr-1)/nxn(inest)
+        ix = igr - jy*nxn(inest) - 1
+
+  ! Interpolate all meteorological data needed for the convection scheme
+        psconv=(psn(ix,jy,1,mind1,inest)*dt2+ &
+             psn(ix,jy,1,mind2,inest)*dt1)*dtt
+        tt2conv=(tt2n(ix,jy,1,mind1,inest)*dt2+ &
+             tt2n(ix,jy,1,mind2,inest)*dt1)*dtt
+        td2conv=(td2n(ix,jy,1,mind1,inest)*dt2+ &
+             td2n(ix,jy,1,mind2,inest)*dt1)*dtt
+!!$        do kz=1,nconvlev+1    !old
+        do kz=1,nuvz-1           !bugfix
+          tconv(kz)=(tthn(ix,jy,kz+1,mind1,inest)*dt2+ &
+               tthn(ix,jy,kz+1,mind2,inest)*dt1)*dtt
+          qconv(kz)=(qvhn(ix,jy,kz+1,mind1,inest)*dt2+ &
+               qvhn(ix,jy,kz+1,mind2,inest)*dt1)*dtt
+        end do
+
+  ! calculate translocation matrix
+  !*******************************
+        call calcmatrix(lconv,delt,cbasefluxn(ix,jy,inest))
+        igrold = igr
+        ktop = 0
+      endif
+
+  ! treat particle only if column has convection
+      if (lconv .eqv. .true.) then
+  ! assign new vertical position to particle
+        ztold=ztra1(ipart)
+        call redist(ipart,ktop,ipconv)
+  !      if (ipconv.le.0) sumconv = sumconv+1
+
+  ! Calculate the gross fluxes across layer interfaces
+  !***************************************************
+
+        if (iflux.eq.1) then
+          itage=abs(itra1(ipart)-itramem(ipart))
+          do nage=1,nageclass
+            if (itage.lt.lage(nage)) goto 47
+          end do
+ 47       continue
+
+          if (nage.le.nageclass) &
+               call calcfluxes(nage,ipart,real(xtra1(ipart)), &
+               real(ytra1(ipart)),ztold)
+        endif
+
+      endif !(lconv .eqv. .true.)
+
+
+60    continue
+    end do
+  end do
+  !--------------------------------------------------------------------------
+  !write(*,*)'############################################'
+  !write(*,*)'TIME=',
+  !    &  itime
+  !write(*,*)'fraction of particles under convection',
+  !    &  sumconv/(sumall+0.001)
+  !write(*,*)'total number of particles',
+  !    &  sumall
+  !write(*,*)'number of particles under convection',
+  !    &  sumconv
+  !write(*,*)'############################################'
+
+  return
+end subroutine convmix
