Index: c/calcmatrix_gfs.f90
===================================================================
--- src/calcmatrix_gfs.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ 	(revision )
@@ -1,139 +1,0 @@
-!**********************************************************************
-! 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 calcmatrix(lconv,delt,cbmf)
-  !                        o    i    o
-  !*****************************************************************************
-  !                                                                            *
-  !  This subroutine calculates the matrix describing convective               *
-  !  redistribution of mass in a grid column, using the subroutine             *
-  !  convect43c.f provided by Kerry Emanuel.                                   *
-  !                                                                            *
-  !  Petra Seibert, Bernd C. Krueger, 2000-2001                                *
-  !                                                                            *
-  !  changed by C. Forster, November 2003 - February 2004                      *
-  !  array fmassfrac(nconvlevmax,nconvlevmax) represents                       *
-  !  the convective redistribution matrix for the particles                    *
-  !                                                                            *
-  !  Changes by C. Forster, November 2005, NCEP GFS version                    *
-  !                                                                            *
-  !  lconv        indicates whether there is convection in this cell, or not   *
-  !  delt         time step for convection [s]                                 *
-  !  cbmf         cloud base mass flux                                         *
-  !                                                                            *
-  !*****************************************************************************
-
-  use par_mod
-  use com_mod
-  use conv_mod
-
-  implicit none
-
-  real :: rlevmass,summe
-
-  integer :: iflag, k, kk, kuvz
-
-  !1-d variables for convection
-  !variables for redistribution matrix
-  real :: cbmfold, precip, qprime
-  real :: tprime, wd, f_qvsat
-  real :: delt,cbmf
-  logical :: lconv
-
-  lconv = .false.
-
-
-  ! calculate pressure at eta levels for use in convect
-  ! and assign temp & spec. hum. to 1D workspace
-  ! -------------------------------------------------------
-
-  ! pconv(1) is the pressure at the first level above ground
-  ! phconv(k) is the pressure between levels k-1 and k
-  ! dpr(k) is the pressure difference "around" tconv(k)
-  ! phconv(kmax) must also be defined 1/2 level above pconv(kmax)
-  ! Therefore, we define k = kuvz-1 and let kuvz start from 2
-  ! top layer cannot be used for convection because p at top of this layer is
-  ! not given
-
-  phconv(1) = psconv
-  do kuvz = 2,nuvz
-    k = kuvz-1
-    phconv(kuvz) =  0.5*(pconv(kuvz)+pconv(k))
-    dpr(k) = phconv(k) - phconv(kuvz)
-    qsconv(k) = f_qvsat( pconv(k), tconv(k) )
-  ! initialize mass fractions
-    do kk=1,nconvlev
-      fmassfrac(k,kk)=0.
-    enddo
-  end do
-
-  !note that Emanuel says it is important
-  !a. to set this =0. every grid point
-  !b. to keep this value in the calling programme in the iteration
-
-  ! CALL CONVECTION
-  !******************
-
-    cbmfold = cbmf
-  ! Convert pressures to hPa, as required by Emanuel scheme
-  !********************************************************
-!!$    do k=1,nconvlev     !old
-    do k=1,nconvlev+1      !bugfix
-      pconv_hpa(k)=pconv(k)/100.
-      phconv_hpa(k)=phconv(k)/100.
-    end do
-    phconv_hpa(nconvlev+1)=phconv(nconvlev+1)/100.
-    call convect(nconvlevmax, nconvlev, delt, iflag, &
-         precip, wd, tprime, qprime, cbmf)
-
-  ! do not update fmassfrac and cloudbase massflux
-  ! if no convection takes place or
-  ! if a CFL criterion is violated in convect43c.f
-   if (iflag .ne. 1 .and. iflag .ne. 4) then
-     cbmf=cbmfold
-     goto 200
-   endif
-
-  ! do not update fmassfrac and cloudbase massflux
-  ! if the old and the new cloud base mass
-  ! fluxes are zero
-   if (cbmf.le.0..and.cbmfold.le.0.) then
-     cbmf=cbmfold
-     goto 200
-   endif
-
-  ! Update fmassfrac
-  ! account for mass displaced from level k to level k
-
-   lconv = .true.
-    do k=1,nconvtop
-      rlevmass = dpr(k)/ga
-      summe = 0.
-      do kk=1,nconvtop
-        fmassfrac(k,kk) = delt*fmass(k,kk)
-        summe = summe + fmassfrac(k,kk)
-      end do
-      fmassfrac(k,k)=fmassfrac(k,k) + rlevmass - summe
-    end do
-
-200   continue
-
-end subroutine calcmatrix
Index: c/calcpar_gfs.f90
===================================================================
--- src/calcpar_gfs.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ 	(revision )
@@ -1,243 +1,0 @@
-!**********************************************************************
-! 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 calcpar(n,uuh,vvh,pvh)
-  !                   i  i   i   o
-  !*****************************************************************************
-  !                                                                            *
-  !     Computation of several boundary layer parameters needed for the        *
-  !     dispersion calculation and calculation of dry deposition velocities.   *
-  !     All parameters are calculated over the entire grid.                    *
-  !                                                                            *
-  !     Author: A. Stohl                                                       *
-  !                                                                            *
-  !     21 May 1995                                                            *
-  !                                                                            *
-  ! ------------------------------------------------------------------         *
-  !     Petra Seibert, Feb 2000:                                               *
-  !     convection scheme:                                                     *
-  !     new variables in call to richardson                                    *
-  !                                                                            *
-  !*****************************************************************************
-  !  Changes, Bernd C. Krueger, Feb. 2001:
-  !   Variables tth and qvh (on eta coordinates) in common block
-  !*****************************************************************************
-  !                                                                            *
-  !  CHANGE 17/11/2005 Caroline Forster NCEP GFS version                       *
-  !                                                                            *
-  ! Variables:                                                                 *
-  ! n                  temporal index for meteorological fields (1 to 3)       *
-  !                                                                            *
-  ! Constants:                                                                 *
-  !                                                                            *
-  !                                                                            *
-  ! Functions:                                                                 *
-  ! scalev             computation of ustar                                    *
-  ! obukhov            computatio of Obukhov length                            *
-  !                                                                            *
-  !*****************************************************************************
-
-  use par_mod
-  use com_mod
-
-  implicit none
-
-  integer :: n,ix,jy,i,kz,lz,kzmin,llev
-  real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
-  real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
-  real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax),hmixdummy
-  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
-  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
-  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
-  real,parameter :: const=r_air/ga
-
-
-  ! Loop over entire grid
-  !**********************
-
-  do jy=0,nymin1
-
-  ! Set minimum height for tropopause
-  !**********************************
-
-    ylat=ylat0+real(jy)*dy
-    if ((ylat.ge.-20.).and.(ylat.le.20.)) then
-      altmin = 5000.
-    else
-      if ((ylat.gt.20.).and.(ylat.lt.40.)) then
-        altmin=2500.+(40.-ylat)*125.
-      else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
-        altmin=2500.+(40.+ylat)*125.
-      else
-        altmin=2500.
-      endif
-    endif
-
-    do ix=0,nxmin1
-
-  ! 1) Calculation of friction velocity
-  !************************************
-
-      ustar(ix,jy,1,n)=scalev(ps(ix,jy,1,n),tt2(ix,jy,1,n), &
-           td2(ix,jy,1,n),surfstr(ix,jy,1,n))
-      if (ustar(ix,jy,1,n).le.1.e-8) ustar(ix,jy,1,n)=1.e-8
-
-  ! 2) Calculation of inverse Obukhov length scale
-  !***********************************************
-
-  ! NCEP version: find first level above ground
-      llev = 0
-      do i=1,nuvz
-       if (ps(ix,jy,1,n).lt.akz(i)) llev=i
-      end do
-       llev = llev+1
-       if (llev.gt.nuvz) llev = nuvz-1
-  ! NCEP version
-
-  ! calculate inverse Obukhov length scale with tth(llev)
-      ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
-           tth(ix,jy,llev,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akz(llev))
-      if (ol.ne.0.) then
-        oli(ix,jy,1,n)=1./ol
-      else
-        oli(ix,jy,1,n)=99999.
-      endif
-
-
-  ! 3) Calculation of convective velocity scale and mixing height
-  !**************************************************************
-
-      do i=1,nuvz
-        ulev(i)=uuh(ix,jy,i)
-        vlev(i)=vvh(ix,jy,i)
-        ttlev(i)=tth(ix,jy,i,n)
-        qvlev(i)=qvh(ix,jy,i,n)
-      end do
-
-  ! NCEP version hmix has been read in in readwind.f, is therefore not calculated here
-      call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
-           ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), &
-           td2(ix,jy,1,n),hmixdummy,wstar(ix,jy,1,n),hmixplus)
-
-      if(lsubgrid.eq.1) then
-        subsceff=min(excessoro(ix,jy),hmixplus)
-      else
-        subsceff=0
-      endif
-  !
-  ! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
-  !
-      hmix(ix,jy,1,n)=hmix(ix,jy,1,n)+subsceff
-      hmix(ix,jy,1,n)=max(hmixmin,hmix(ix,jy,1,n)) ! set minimum PBL height
-      hmix(ix,jy,1,n)=min(hmixmax,hmix(ix,jy,1,n)) ! set maximum PBL height
-
-  ! 4) Calculation of dry deposition velocities
-  !********************************************
-
-      if (DRYDEP) then
-  ! Sabine Eckhardt, Dec 06: use new index for z0 for water depending on
-  ! windspeed
-        z0(7)=0.016*ustar(ix,jy,1,n)*ustar(ix,jy,1,n)/ga
-
-  ! Calculate relative humidity at surface
-  !***************************************
-        rh=ew(td2(ix,jy,1,n))/ew(tt2(ix,jy,1,n))
-
-        call getvdep(n,ix,jy,ustar(ix,jy,1,n), &
-             tt2(ix,jy,1,n),ps(ix,jy,1,n),1./oli(ix,jy,1,n), &
-             ssr(ix,jy,1,n),rh,lsprec(ix,jy,1,n)+convprec(ix,jy,1,n), &
-             sd(ix,jy,1,n),vd)
-
-        do i=1,nspec
-          vdep(ix,jy,i,n)=vd(i)
-        end do
-
-      endif
-
-  !******************************************************
-  ! Calculate height of thermal tropopause (Hoinka, 1997)
-  !******************************************************
-
-  ! 1) Calculate altitudes of NCEP model levels
-  !*********************************************
-
-      tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
-           ps(ix,jy,1,n))
-      pold=ps(ix,jy,1,n)
-      zold=0.
-      do kz=llev,nuvz
-        pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)  ! pressure on model layers
-        tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
-
-        if (abs(tv-tvold).gt.0.2) then
-         zlev(kz)=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
-        else
-          zlev(kz)=zold+const*log(pold/pint)*tv
-        endif
-        tvold=tv
-        pold=pint
-        zold=zlev(kz)
-      end do
-
-  ! 2) Define a minimum level kzmin, from which upward the tropopause is
-  !    searched for. This is to avoid inversions in the lower troposphere
-  !    to be identified as the tropopause
-  !************************************************************************
-
-      do kz=llev,nuvz
-        if (zlev(kz).ge.altmin) then
-          kzmin=kz
-          goto 45
-        endif
-      end do
-45    continue
-
-  ! 3) Search for first stable layer above minimum height that fulfills the
-  !    thermal tropopause criterion
-  !************************************************************************
-
-      do kz=kzmin,nuvz
-        do lz=kz+1,nuvz
-          if ((zlev(lz)-zlev(kz)).gt.2000.) then
-            if (((tth(ix,jy,kz,n)-tth(ix,jy,lz,n))/ &
-                 (zlev(lz)-zlev(kz))).lt.0.002) then
-              tropopause(ix,jy,1,n)=zlev(kz)
-              goto 51
-            endif
-            goto 50
-          endif
-        end do
-50      continue
-      end do
-51    continue
-
-
-    end do
-  end do
-
-
-  ! Calculation of potential vorticity on 3-d grid
-  !***********************************************
-
-  call calcpv(n,uuh,vvh,pvh)
-
-
-end subroutine calcpar
Index: c/convmix_gfs.f90
===================================================================
--- src/convmix_gfs.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ 	(revision )
@@ -1,303 +1,0 @@
-!**********************************************************************
-! 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
Index: c/ecmwf_mod.f90
===================================================================
--- src/ecmwf_mod.f90	(revision 842074e7416f59098921221d89faf6edc3639279)
+++ 	(revision )
@@ -1,68 +1,0 @@
-!**********************************************************************
-! 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/>.   *
-!**********************************************************************
-
-!*******************************************************************************
-!   Include file for calculation of particle trajectories (Program FLEXPART)   *
-!                                                                              *
-!       This file contains ECMWF specific parameters used in FLEXPART          *
-!       Note that module name differs from file name.                          *
-!       The makefile selects either this file, or gfs_mod.f90, depending       *
-!       on target.                                                             *
-!                                                                              *
-!       Author: ESO                                                            *
-!                                                                              *
-!       2015                                                                   *
-!                                                                              *
-!*******************************************************************************
-
-module wind_mod
-
-  implicit none
-
-  !*********************************************
-  ! Maximum dimensions of the input mother grids
-  !*********************************************
-
-!  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF new 
-   integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 
-!  integer,parameter :: nxmax=721,nymax=361,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 0.5 
-   integer,parameter :: nxshift=359
-!  integer,parameter :: nxshift=718
-!  integer,parameter :: nxshift=0
-
-  !*********************************************
-  ! Maximum dimensions of the nested input grids
-  !*********************************************
-
-  integer,parameter :: maxnests=0,nxmaxn=361,nymaxn=181
-
-  ! nxmax,nymax        maximum dimension of wind fields in x and y
-  !                    direction, respectively
-  ! nuvzmax,nwzmax     maximum dimension of (u,v) and (w) wind fields in z
-  !                    direction (for fields on eta levels)
-  ! nzmax              maximum dimension of wind fields in z direction
-  !                    for the transformed Cartesian coordinates
-  ! nxshift            for global grids (in x), the grid can be shifted by
-  !                    nxshift grid points, in order to accomodate nested
-  !                    grids, and output grids overlapping the domain "boundary"
-  !                    nxshift must not be negative; "normal" setting would be 0
-
-end module wind_mod
Index: c/gfs_mod.f90
===================================================================
--- src/gfs_mod.f90	(revision 4c64400c622119505b289ef83be130b0b81b9590)
+++ 	(revision )
@@ -1,64 +1,0 @@
-!**********************************************************************
-! 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/>.   *
-!**********************************************************************
-
-!*******************************************************************************
-!   Include file for calculation of particle trajectories (Program FLEXPART)   *
-!                                                                              *
-!       This file contains GFS specific parameters used in FLEXPART            *
-!       Note that module name differs from file name.                          *
-!       The makefile selects either this file, or ecmwf_mod.f90, depending     *
-!       on target.                                                             *
-!                                                                              *
-!       Author: ESO                                                            *
-!                                                                              *
-!       2015                                                                   *
-!                                                                              *
-!*******************************************************************************
-
-module wind_mod
-
-  implicit none
-
-  !*********************************************
-  ! Maximum dimensions of the input mother grids
-  !*********************************************
-
-!  integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
-  integer,parameter :: nxmax=361,nymax=181,nuvzmax=64,nwzmax=64,nzmax=64
-  integer,parameter :: nxshift=0     ! for GFS or FNL
-
-  !*********************************************
-  ! Maximum dimensions of the nested input grids
-  !*********************************************
-
-  integer,parameter :: maxnests=0,nxmaxn=361,nymaxn=181
-
-  ! nxmax,nymax        maximum dimension of wind fields in x and y
-  !                    direction, respectively
-  ! nuvzmax,nwzmax     maximum dimension of (u,v) and (w) wind fields in z
-  !                    direction (for fields on eta levels)
-  ! nzmax              maximum dimension of wind fields in z direction
-  !                    for the transformed Cartesian coordinates
-  ! nxshift            for global grids (in x), the grid can be shifted by
-  !                    nxshift grid points, in order to accomodate nested
-  !                    grids, and output grids overlapping the domain "boundary"
-  !                    nxshift must not be negative; "normal" setting would be 0
-end module wind_mod
Index: c/gridcheck.f90
===================================================================
--- src/gridcheck.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
+++ 	(revision )
@@ -1,571 +1,0 @@
-!**********************************************************************
-! 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 gridcheck
-
-  !**********************************************************************
-  !                                                                     *
-  !             FLEXPART MODEL SUBROUTINE GRIDCHECK                     *
-  !                                                                     *
-  !**********************************************************************
-  !                                                                     *
-  !             AUTHOR:      G. WOTAWA                                  *
-  !             DATE:        1997-08-06                                 *
-  !             LAST UPDATE: 1997-10-10                                 *
-  !                                                                     *
-  !             Update:      1999-02-08, global fields allowed, A. Stohl*
-  !             CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
-  !                                 ECMWF grib_api                      *
-  !             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
-  !                                 ECMWF grib_api                      *
-  !                                                                     *
-  !**********************************************************************
-  !                                                                     *
-  ! DESCRIPTION:                                                        *
-  !                                                                     *
-  ! THIS SUBROUTINE DETERMINES THE GRID SPECIFICATIONS (LOWER LEFT      *
-  ! LONGITUDE, LOWER LEFT LATITUDE, NUMBER OF GRID POINTS, GRID DIST-   *
-  ! ANCE AND VERTICAL DISCRETIZATION OF THE ECMWF MODEL) FROM THE       *
-  ! GRIB HEADER OF THE FIRST INPUT FILE. THE CONSISTANCY (NO CHANGES    *
-  ! WITHIN ONE FLEXPART RUN) IS CHECKED IN THE ROUTINE "READWIND" AT    *
-  ! ANY CALL.                                                           *
-  !                                                                     *
-  ! XLON0                geographical longitude of lower left gridpoint *
-  ! YLAT0                geographical latitude of lower left gridpoint  *
-  ! NX                   number of grid points x-direction              *
-  ! NY                   number of grid points y-direction              *
-  ! DX                   grid distance x-direction                      *
-  ! DY                   grid distance y-direction                      *
-  ! NUVZ                 number of grid points for horizontal wind      *
-  !                      components in z direction                      *
-  ! NWZ                  number of grid points for vertical wind        *
-  !                      component in z direction                       *
-  ! sizesouth, sizenorth give the map scale (i.e. number of virtual grid*
-  !                      points of the polar stereographic grid):       *
-  !                      used to check the CFL criterion                *
-  ! UVHEIGHT(1)-         heights of gridpoints where u and v are        *
-  ! UVHEIGHT(NUVZ)       given                                          *
-  ! WHEIGHT(1)-          heights of gridpoints where w is given         *
-  ! WHEIGHT(NWZ)                                                        *
-  !                                                                     *
-  !**********************************************************************
-
-  use grib_api
-  use par_mod
-  use com_mod
-  use conv_mod
-  use cmapf_mod, only: stlmbr,stcm2p
-
-  implicit none
-
-  !HSO  parameters for grib_api
-  integer :: ifile
-  integer :: iret
-  integer :: igrib
-  integer :: gotGrid
-  real(kind=4) :: xaux1,xaux2,yaux1,yaux2
-  real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
-  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
-  !HSO  end
-  integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
-  real :: sizesouth,sizenorth,xauxa,pint,conversion_factor
-
-  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
-
-  ! dimension of isec2 at least (22+n), where n is the number of parallels or
-  ! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
-
-  ! dimension of zsec2 at least (10+nn), where nn is the number of vertical
-  ! coordinate parameters
-
-  integer :: isec1(56),isec2(22+nxmax+nymax)
-  real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
-  character(len=1) :: opt
-
-  !HSO  grib api error messages
-  character(len=24) :: gribErrorMsg = 'Error reading grib file'
-  character(len=20) :: gribFunction = 'gridcheck'
-
-
-  iumax=0
-  iwmax=0
-
-  if(ideltas.gt.0) then
-    ifn=1
-  else
-    ifn=numbwf
-  endif
-  !
-  ! OPENING OF DATA FILE (GRIB CODE)
-  !
-5   call grib_open_file(ifile,path(3)(1:length(3)) &
-         //trim(wfname(ifn)),'r',iret)
-  if (iret.ne.GRIB_SUCCESS) then
-    goto 999   ! ERROR DETECTED
-  endif
-  !turn on support for multi fields messages
-  !call grib_multi_support_on
-
-  gotGrid=0
-  ifield=0
-10   ifield=ifield+1
-
-  !
-  ! GET NEXT FIELDS
-  !
-  call grib_new_from_file(ifile,igrib,iret)
-  if (iret.eq.GRIB_END_OF_FILE )  then
-    goto 30    ! EOF DETECTED
-  elseif (iret.ne.GRIB_SUCCESS) then
-    goto 999   ! ERROR DETECTED
-  endif
-
-  !first see if we read GRIB1 or GRIB2
-  call grib_get_int(igrib,'editionNumber',gribVer,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-
-  if (gribVer.eq.1) then ! GRIB Edition 1
-
-  !print*,'GRiB Edition 1'
-  !read the grib2 identifiers
-  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'level',isec1(8),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-
-  !change code for etadot to code for omega
-  if (isec1(6).eq.77) then
-    isec1(6)=135
-  endif
-
-  !print*,isec1(6),isec1(8)
-
-  else
-
-  !print*,'GRiB Edition 2'
-  !read the grib2 identifiers
-  call grib_get_int(igrib,'discipline',discipl,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'parameterCategory',parCat,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'parameterNumber',parNum,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'level',valSurf,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'paramId',parId,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-
-  !print*,discipl,parCat,parNum,typSurf,valSurf
-
-  !convert to grib1 identifiers
-  isec1(6)=-1
-  isec1(7)=-1
-  isec1(8)=-1
-  isec1(8)=valSurf     ! level
-  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
-    isec1(6)=130         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
-    isec1(6)=131         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
-    isec1(6)=132         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
-    isec1(6)=133         ! indicatorOfParameter
-!ZHG FOR CLOUDS FROM GRIB
-  elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
-    isec1(6)=246         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
-    isec1(6)=247         ! indicatorOfParameter
-!ZHG end
-! ESO qc(=clwc+ciwc)
-  elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
-    isec1(6)=201031      ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
-    isec1(6)=134         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
-    isec1(6)=135         ! indicatorOfParameter
-  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
-    isec1(6)=135         ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
-    isec1(6)=151         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
-    isec1(6)=165         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
-    isec1(6)=166         ! indicatorOfParameter
-  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
-    isec1(6)=167         ! indicatorOfParameter
-  elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
-    isec1(6)=168         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
-    isec1(6)=141         ! indicatorOfParameter
-  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
-    isec1(6)=164         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
-    isec1(6)=142         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
-    isec1(6)=143         ! indicatorOfParameter
-  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
-    isec1(6)=146         ! indicatorOfParameter
-  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
-    isec1(6)=176         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
-    isec1(6)=180         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
-    isec1(6)=181         ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
-    isec1(6)=129         ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
-    isec1(6)=160         ! indicatorOfParameter
-  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
-       (typSurf.eq.1)) then ! LSM
-    isec1(6)=172         ! indicatorOfParameter
-  else
-    print*,'***ERROR: undefined GRiB2 message found!',discipl, &
-         parCat,parNum,typSurf
-  endif
-  if(parId .ne. isec1(6) .and. parId .ne. 77) then
-    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
-!    stop
-  endif
-
-  endif
-
-  !get the size and data of the values array
-  if (isec1(6).ne.-1) then
-    call grib_get_real4_array(igrib,'values',zsec4,iret)
-    call grib_check(iret,gribFunction,gribErrorMsg)
-  endif
-
-  if (ifield.eq.1) then
-
-  !HSO  get the required fields from section 2 in a gribex compatible manner
-    call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
-         isec2(2),iret)
-    call grib_check(iret,gribFunction,gribErrorMsg)
-    call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
-         isec2(3),iret)
-    call grib_check(iret,gribFunction,gribErrorMsg)
-    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
-         xaux1in,iret)
-    call grib_check(iret,gribFunction,gribErrorMsg)
-    call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
-         isec2(12),iret)
-    call grib_check(iret,gribFunction,gribErrorMsg)
-
-  !  get the size and data of the vertical coordinate array
-    call grib_get_real4_array(igrib,'pv',zsec2,iret)
-    call grib_check(iret,gribFunction,gribErrorMsg)
-
-    nxfield=isec2(2)
-    ny=isec2(3)
-    nlev_ec=isec2(12)/2-1
-   endif
-
-  !HSO  get the second part of the grid dimensions only from GRiB1 messages
-  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
-    call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
-         xaux2in,iret)
-    call grib_check(iret,gribFunction,gribErrorMsg)
-    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
-         yaux1in,iret)
-    call grib_check(iret,gribFunction,gribErrorMsg)
-    call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
-         yaux2in,iret)
-    call grib_check(iret,gribFunction,gribErrorMsg)
-    xaux1=xaux1in
-    xaux2=xaux2in
-    yaux1=yaux1in
-    yaux2=yaux2in
-    if (xaux1.gt.180.) xaux1=xaux1-360.0
-    if (xaux2.gt.180.) xaux2=xaux2-360.0
-    if (xaux1.lt.-180.) xaux1=xaux1+360.0
-    if (xaux2.lt.-180.) xaux2=xaux2+360.0
-    if (xaux2.lt.xaux1) xaux2=xaux2+360.0
-    xlon0=xaux1
-    ylat0=yaux1
-    dx=(xaux2-xaux1)/real(nxfield-1)
-    dy=(yaux2-yaux1)/real(ny-1)
-    dxconst=180./(dx*r_earth*pi)
-    dyconst=180./(dy*r_earth*pi)
-    gotGrid=1
-  ! Check whether fields are global
-  ! If they contain the poles, specify polar stereographic map
-  ! projections using the stlmbr- and stcm2p-calls
-  !***********************************************************
-
-    xauxa=abs(xaux2+dx-360.-xaux1)
-    if (xauxa.lt.0.001) then
-      nx=nxfield+1                 ! field is cyclic
-      xglobal=.true.
-      if (abs(nxshift).ge.nx) &
-           stop 'nxshift in file par_mod is too large'
-      xlon0=xlon0+real(nxshift)*dx
-    else
-      nx=nxfield
-      xglobal=.false.
-      if (nxshift.ne.0) &
-           stop 'nxshift (par_mod) must be zero for non-global domain'
-    endif
-    nxmin1=nx-1
-    nymin1=ny-1
-    if (xlon0.gt.180.) xlon0=xlon0-360.
-    xauxa=abs(yaux1+90.)
-    if (xglobal.and.xauxa.lt.0.001) then
-      sglobal=.true.               ! field contains south pole
-  ! Enhance the map scale by factor 3 (*2=6) compared to north-south
-  ! map scale
-      sizesouth=6.*(switchsouth+90.)/dy
-      call stlmbr(southpolemap,-90.,0.)
-      call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
-           sizesouth,switchsouth,180.)
-      switchsouthg=(switchsouth-ylat0)/dy
-    else
-      sglobal=.false.
-      switchsouthg=999999.
-    endif
-    xauxa=abs(yaux2-90.)
-    if (xglobal.and.xauxa.lt.0.001) then
-      nglobal=.true.               ! field contains north pole
-  ! Enhance the map scale by factor 3 (*2=6) compared to north-south
-  ! map scale
-      sizenorth=6.*(90.-switchnorth)/dy
-      call stlmbr(northpolemap,90.,0.)
-      call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
-           sizenorth,switchnorth,180.)
-      switchnorthg=(switchnorth-ylat0)/dy
-    else
-      nglobal=.false.
-      switchnorthg=999999.
-    endif
-    if (nxshift.lt.0) &
-         stop 'nxshift (par_mod) must not be negative'
-    if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
-  endif ! gotGrid
-
-  if (nx.gt.nxmax) then
-   write(*,*) 'FLEXPART error: Too many grid points in x direction.'
-    write(*,*) 'Reduce resolution of wind fields.'
-    write(*,*) 'Or change parameter settings in file par_mod.'
-    write(*,*) nx,nxmax
-    stop
-  endif
-
-  if (ny.gt.nymax) then
-   write(*,*) 'FLEXPART error: Too many grid points in y direction.'
-    write(*,*) 'Reduce resolution of wind fields.'
-    write(*,*) 'Or change parameter settings in file par_mod.'
-    write(*,*) ny,nymax
-    stop
-  endif
-
-  k=isec1(8)
-  if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
-  if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
-
-  if(isec1(6).eq.129) then
-    do jy=0,ny-1
-      do ix=0,nxfield-1
-        oro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)/ga
-      end do
-    end do
-  endif
-  if(isec1(6).eq.172) then
-    do jy=0,ny-1
-      do ix=0,nxfield-1
-        lsm(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
-      end do
-    end do
-  endif
-  if(isec1(6).eq.160) then
-    do jy=0,ny-1
-      do ix=0,nxfield-1
-        excessoro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
-      end do
-    end do
-  endif
-
-  call grib_release(igrib)
-  goto 10                      !! READ NEXT LEVEL OR PARAMETER
-  !
-  ! CLOSING OF INPUT DATA FILE
-  !
-
-30   call grib_close_file(ifile)
-
-  !error message if no fields found with correct first longitude in it
-  if (gotGrid.eq.0) then
-    print*,'***ERROR: input file needs to contain GRiB1 formatted'// &
-         'messages'
-    stop
-  endif
-
-  nuvz=iumax
-  nwz =iwmax
-  if(nuvz.eq.nlev_ec) nwz=nlev_ec+1
-
-  if (nuvz+1.gt.nuvzmax) then
-    write(*,*) 'FLEXPART error: Too many u,v grid points in z '// &
-         'direction.'
-    write(*,*) 'Reduce resolution of wind fields.'
-    write(*,*) 'Or change parameter settings in file par_mod.'
-    write(*,*) nuvz+1,nuvzmax
-    stop
-  endif
-
-  if (nwz.gt.nwzmax) then
-    write(*,*) 'FLEXPART error: Too many w grid points in z '// &
-         'direction.'
-    write(*,*) 'Reduce resolution of wind fields.'
-    write(*,*) 'Or change parameter settings in file par_mod.'
-    write(*,*) nwz,nwzmax
-    stop
-  endif
-
-  ! If desired, shift all grids by nxshift grid cells
-  !**************************************************
-
-  if (xglobal) then
-    call shift_field_0(oro,nxfield,ny)
-    call shift_field_0(lsm,nxfield,ny)
-    call shift_field_0(excessoro,nxfield,ny)
-  endif
-
-  ! Output of grid info
-  !********************
-
-  if (lroot) then
-    write(*,'(a,2i7)') ' Vertical levels in ECMWF data: ', &
-         nuvz+1,nwz
-    write(*,*)
-    write(*,'(a)') ' Mother domain:'
-    write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
-         xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
-    write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range : ', &
-         ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
-    write(*,*)
-  end if
-
-  ! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
-  ! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
-
-  numskip=nlev_ec-nuvz  ! number of ecmwf model layers not used
-                        ! by trajectory model
-  !do 8940 i=1,244
-  !   write (*,*) 'zsec2:',i,ifield,zsec2(i),numskip
-  !940  continue
-  !   stop
-  ! SEC SEC SEC
-  ! for unknown reason zsec 1 to 10 is filled in this version
-  ! compared to the old one
-  ! SEC SEC SE
-  do i=1,nwz
-    j=numskip+i
-    k=nlev_ec+1+numskip+i
-    akm(nwz-i+1)=zsec2(j)
-  !   write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
-    bkm(nwz-i+1)=zsec2(k)
-  end do
-
-  !
-  ! CALCULATION OF AKZ, BKZ
-  ! AKZ,BKZ: model discretization parameters at the center of each model
-  !     layer
-  !
-  ! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0,
-  ! i.e. ground level
-  !*****************************************************************************
-
-  akz(1)=0.
-  bkz(1)=1.0
-  do i=1,nuvz
-     akz(i+1)=0.5*(akm(i+1)+akm(i))
-     bkz(i+1)=0.5*(bkm(i+1)+bkm(i))
-  end do
-  nuvz=nuvz+1
-
-  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
-  ! upon the transformation to z levels. In order to save computer memory, this is
-  ! not done anymore in the standard version. However, this option can still be
-  ! switched on by replacing the following lines with those below, that are
-  ! currently commented out. For this, similar changes are necessary in
-  ! verttransform.f and verttranform_nests.f
-  !*****************************************************************************
-
-  nz=nuvz
-  if (nz.gt.nzmax) stop 'nzmax too small'
-  do i=1,nuvz
-    aknew(i)=akz(i)
-    bknew(i)=bkz(i)
-  end do
-
-  ! Switch on following lines to use doubled vertical resolution
-  !*************************************************************
-  !nz=nuvz+nwz-1
-  !if (nz.gt.nzmax) stop 'nzmax too small'
-  !do 100 i=1,nwz
-  !  aknew(2*(i-1)+1)=akm(i)
-  !00     bknew(2*(i-1)+1)=bkm(i)
-  !do 110 i=2,nuvz
-  !  aknew(2*(i-1))=akz(i)
-  !10     bknew(2*(i-1))=bkz(i)
-  ! End doubled vertical resolution
-
-
-  ! Determine the uppermost level for which the convection scheme shall be applied
-  ! by assuming that there is no convection above 50 hPa (for standard SLP)
-  !*****************************************************************************
-
-  do i=1,nuvz-2
-    pint=akz(i)+bkz(i)*101325.
-    if (pint.lt.5000.) goto 96
-  end do
-96   nconvlev=i
-  if (nconvlev.gt.nconvlevmax-1) then
-    nconvlev=nconvlevmax-1
-    write(*,*) 'Attention, convection only calculated up to ', &
-         akz(nconvlev)+bkz(nconvlev)*1013.25,' hPa'
-  endif
-
-  return
-
-999   write(*,*)
-  write(*,*) ' ###########################################'// &
-       '###### '
-  write(*,*) '       TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
-  write(*,*) ' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
-  write(*,*) ' ###########################################'// &
-       '###### '
-  write(*,*)
-  write(*,'(a)') '!!! PLEASE INSERT A NEW CD-ROM AND   !!!'
-  write(*,'(a)') '!!! PRESS ANY KEY TO CONTINUE...     !!!'
-  write(*,'(a)') '!!! ...OR TERMINATE FLEXPART PRESSING!!!'
-  write(*,'(a)') '!!! THE "X" KEY...                   !!!'
-  write(*,*)
-  read(*,'(a)') opt
-  if(opt.eq.'X') then
-    stop
-  else
-    goto 5
-  endif
-
-end subroutine gridcheck
-
Index: c/obukhov_gfs.f90
===================================================================
--- src/obukhov_gfs.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ 	(revision )
@@ -1,76 +1,0 @@
-!**********************************************************************
-! 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/>.   *
-!**********************************************************************
-
-real function obukhov(ps,tsurf,tdsurf,tlev,ustar,hf,plev)
-
-  !********************************************************************
-  !                                                                   *
-  !                       Author: G. WOTAWA                           *
-  !                       Date:   1994-06-27                          *
-  !                                                                   *
-  !  Update: A. Stohl, 2000-09-25, avoid division by zero by          *
-  !  setting ustar to minimum value                                   *
-  !                                                                   *
-  !  CHANGE: 17/11/2005 Caroline Forster NCEP GFS version             *
-  !                                                                   *
-  !********************************************************************
-  !                                                                   *
-  !     This program calculates Obukhov scale height from surface     *
-  !     meteorological data and sensible heat flux.                   *
-  !                                                                   *
-  !********************************************************************
-  !                                                                   *
-  !     INPUT:                                                        *
-  !                                                                   *
-  !     ps      surface pressure [Pa]                                 *
-  !     tsurf   surface temperature [K]                               *
-  !     tdsurf  surface dew point [K]                                 *
-  !     tlev    temperature first model level [K]                     *
-  !     ustar   scale velocity [m/s]                                  *
-  !     hf      surface sensible heat flux [W/m2]                     *
-  !     akm     ECMWF vertical discretization parameter               *
-  !     bkm     ECMWF vertical discretization parameter               *
-  !                                                                   *
-  !********************************************************************
-
-  use par_mod
-
-  implicit none
-
-  real :: ps,tsurf,tdsurf,tlev,ustar,hf,e,ew,tv,rhoa,plev
-  real :: theta,thetastar
-
-
-  e=ew(tdsurf)                           ! vapor pressure
-  tv=tsurf*(1.+0.378*e/ps)               ! virtual temperature
-  rhoa=ps/(r_air*tv)                      ! air density
-  theta=tlev*(100000./plev)**(r_air/cpa) ! potential temperature
-  if (ustar.le.0.) ustar=1.e-8
-  thetastar=hf/(rhoa*cpa*ustar)           ! scale temperature
-  if(abs(thetastar).gt.1.e-10) then
-     obukhov=theta*ustar**2/(karman*ga*thetastar)
-  else
-     obukhov=9999                        ! zero heat flux
-  endif
-  if (obukhov.gt. 9999.) obukhov= 9999.
-  if (obukhov.lt.-9999.) obukhov=-9999.
-
-end function obukhov
Index: c/richardson_gfs.f90
===================================================================
--- src/richardson_gfs.f90	(revision c2162ce5e4e231a653816e16d96eacf839a90235)
+++ 	(revision )
@@ -1,196 +1,0 @@
-!**********************************************************************
-! 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 richardson(psurf,ust,ttlev,qvlev,ulev,vlev,nuvz, &
-       akz,bkz,hf,tt2,td2,h,wst,hmixplus)
-  !                        i    i    i     i    i    i    i
-  ! i   i  i   i   i  o  o     o
-  !****************************************************************************
-  !                                                                           *
-  !     Calculation of mixing height based on the critical Richardson number. *
-  !     Calculation of convective time scale.                                 *
-  !     For unstable conditions, one iteration is performed. An excess        *
-  !     temperature (dependent on hf and wst) is calculated, added to the     *
-  !     temperature at the lowest model level. Then the procedure is repeated.*
-  !                                                                           *
-  !     Author: A. Stohl                                                      *
-  !                                                                           *
-  !     22 August 1996                                                        *
-  !                                                                           *
-  !     Literature:                                                           *
-  !     Vogelezang DHP and Holtslag AAM (1996): Evaluation and model impacts  *
-  !     of alternative boundary-layer height formulations. Boundary-Layer     *
-  !     Meteor. 81, 245-269.                                                  *
-  !                                                                           *
-  !     Update: 1999-02-01 by G. Wotawa                                       *
-  !     CHANGE: 17/11/2005 Caroline Forster NCEP GFS version                  *
-  !                                                                           *
-  !     Two meter level (temperature, humidity) is taken as reference level   *
-  !     instead of first model level.                                         *
-  !     New input variables tt2, td2 introduced.                              *
-  !                                                                           *
-  !****************************************************************************
-  !                                                                           *
-  ! Variables:                                                                *
-  ! h                          mixing height [m]                              *
-  ! hf                         sensible heat flux                             *
-  ! psurf                      surface pressure at point (xt,yt) [Pa]         *
-  ! tv                         virtual temperature                            *
-  ! wst                        convective velocity scale                      *
-  !                                                                           *
-  ! Constants:                                                                *
-  ! ric                        critical Richardson number                     *
-  !                                                                           *
-  !****************************************************************************
-
-  use par_mod
-
-  implicit none
-
-  integer :: i,k,nuvz,iter,llev
-  real :: tv,tvold,zref,z,zold,pint,pold,theta,thetaref,ri
-  real :: akz(nuvz),bkz(nuvz),ulev(nuvz),vlev(nuvz),hf,wst,tt2,td2,ew
-  real :: psurf,ust,ttlev(nuvz),qvlev(nuvz),h,excess
-  real :: thetaold,zl,ul,vl,thetal,ril,hmixplus,wspeed,bvfsq,bvf
-  real :: f_qvsat,rh,rhold,rhl,theta1,theta2,zl1,zl2,thetam
-  real,parameter :: const=r_air/ga, ric=0.25, b=100., bs=8.5
-  integer,parameter :: itmax=3
-
-  excess=0.0
-  iter=0
-
-  ! NCEP version: find first model level above ground
-  !**************************************************
-
-      llev = 0
-      do i=1,nuvz
-       if (psurf.lt.akz(i)) llev=i
-      end do
-       llev = llev+1
-  ! sec llev should not be 1!
-       if (llev.eq.1) llev = 2
-       if (llev.gt.nuvz) llev = nuvz-1
-  ! NCEP version
-
-
-  ! Compute virtual temperature and virtual potential temperature at
-  ! reference level (2 m)
-  !*****************************************************************
-
-30   iter=iter+1
-
-  pold=psurf
-  tvold=tt2*(1.+0.378*ew(td2)/psurf)
-  zold=2.0
-  zref=zold
-  rhold=ew(td2)/ew(tt2)
-
-  thetaref=tvold*(100000./pold)**(r_air/cpa)+excess
-  thetaold=thetaref
-
-  ! Integrate z up to one level above zt
-  !*************************************
-
-  do k=llev,nuvz
-    pint=akz(k)+bkz(k)*psurf  ! pressure on model layers
-    tv=ttlev(k)*(1.+0.608*qvlev(k))
-
-    if (abs(tv-tvold).gt.0.2) then
-      z=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
-    else
-      z=zold+const*log(pold/pint)*tv
-    endif
-
-    theta=tv*(100000./pint)**(r_air/cpa)
-  ! Petra
-    rh = qvlev(k) / f_qvsat( pint, ttlev(k) )
-
-
-  !alculate Richardson number at each level
-  !****************************************
-
-    ri=ga/thetaref*(theta-thetaref)*(z-zref)/ &
-         max(((ulev(k)-ulev(2))**2+(vlev(k)-vlev(2))**2+b*ust**2),0.1)
-
-  !  addition of second condition: MH should not be placed in an
-  !  unstable layer (PS / Feb 2000)
-    if (ri.gt.ric .and. thetaold.lt.theta) goto 20
-
-    tvold=tv
-    pold=pint
-    rhold=rh
-    thetaold=theta
-    zold=z
-  end do
-  k=k-1 ! ESO: make sure k <= nuvz (ticket #139)
-20   continue
-
-  ! Determine Richardson number between the critical levels
-  !********************************************************
-
-  zl1=zold
-  theta1=thetaold
-  do i=1,20
-    zl=zold+real(i)/20.*(z-zold)
-    ul=ulev(k-1)+real(i)/20.*(ulev(k)-ulev(k-1))
-    vl=vlev(k-1)+real(i)/20.*(vlev(k)-vlev(k-1))
-    thetal=thetaold+real(i)/20.*(theta-thetaold)
-    rhl=rhold+real(i)/20.*(rh-rhold)
-    ril=ga/thetaref*(thetal-thetaref)*(zl-zref)/ &
-         max(((ul-ulev(2))**2+(vl-vlev(2))**2+b*ust**2),0.1)
-    zl2=zl
-    theta2=thetal
-    if (ril.gt.ric) goto 25
-    zl1=zl
-    theta1=thetal
-  end do
-
-25   continue
-  h=zl
-  thetam=0.5*(theta1+theta2)
-  wspeed=sqrt(ul**2+vl**2)                    ! Wind speed at z=hmix
-  bvfsq=(ga/thetam)*(theta2-theta1)/(zl2-zl1) ! Brunt-Vaisala frequency
-                                              ! at z=hmix
-
-  ! Under stable conditions, limit the maximum effect of the subgrid-scale topography
-  ! by the maximum lifting possible from the available kinetic energy
-  !*****************************************************************************
-
-  if(bvfsq.le.0.) then
-    hmixplus=9999.
-  else
-    bvf=sqrt(bvfsq)
-    hmixplus=wspeed/bvf*convke
-  endif
-
-
-  ! Calculate convective velocity scale
-  !************************************
-
-  if (hf.lt.0.) then
-    wst=(-h*ga/thetaref*hf/cpa)**0.333
-    excess=-bs*hf/cpa/wst
-    if (iter.lt.itmax) goto 30
-  else
-    wst=0.
-  endif
-
-end subroutine richardson
