Index: c/gridcheck_emos.f90
===================================================================
--- src/gridcheck_emos.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
+++ 	(revision )
@@ -1,396 +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*
-  !                                                                     *
-  !**********************************************************************
-  !                                                                     *
-  ! 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 par_mod
-  use com_mod
-  use conv_mod
-
-  implicit none
-
-  integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
-  real :: xaux1,xaux2,yaux1,yaux2,sizesouth,sizenorth,xauxa,pint
-
-  ! 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 :: isec0(2),isec1(56),isec2(22+nxmax+nymax),isec3(2)
-  integer :: isec4(64),inbuff(jpack),ilen,ierr,iword,lunit
-  !integer iswap
-  real :: zsec2(60+2*nuvzmax),zsec3(2),zsec4(jpunp)
-  character(len=1) :: opt, yoper = 'D'
-
-
-  iumax=0
-  iwmax=0
-
-  if(ideltas.gt.0) then
-    ifn=1
-  else
-    ifn=numbwf
-  endif
-  !
-  ! OPENING OF DATA FILE (GRIB CODE)
-  !
-5   call pbopen(lunit,path(3)(1:length(3))//wfname(ifn),'r',ierr)
-  if(ierr.lt.0) goto 999
-
-  ifield=0
-10   ifield=ifield+1
-  !
-  ! GET NEXT FIELDS
-  !
-  call pbgrib(lunit,inbuff,jpack,ilen,ierr)
-  if(ierr.eq.-1) goto 30    ! EOF DETECTED
-  if(ierr.lt.-1) goto 999   ! ERROR DETECTED
-
-  ierr=1
-
-  ! Check whether we are on a little endian or on a big endian computer
-  !********************************************************************
-
-  !if (inbuff(1).eq.1112101447) then         ! little endian, swap bytes
-  !  iswap=1+ilen/4
-  !  call swap32(inbuff,iswap)
-  !else if (inbuff(1).ne.1196575042) then    ! big endian
-  !  stop 'subroutine gridcheck: corrupt GRIB data'
-  !endif
-
-  call gribex(isec0,isec1,isec2,zsec2,isec3,zsec3,isec4, &
-       zsec4,jpunp,inbuff,jpack,iword,yoper,ierr)
-  if (ierr.ne.0) goto 999   ! ERROR DETECTED
-
-  if(ifield.eq.1) then
-    nxfield=isec2(2)
-    ny=isec2(3)
-    xaux1=real(isec2(5))/1000.
-    xaux2=real(isec2(8))/1000.
-    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.
-    yaux1=real(isec2(7))/1000.
-    yaux2=real(isec2(4))/1000.
-    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)
-    nlev_ec=isec2(12)/2-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
-  endif
-
-  if (nxshift.lt.0) stop 'nxshift (par_mod) must not be negative'
-  if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
-
-  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
-  !       write (*,*) 'ich stop!',nxfield,ny,jy,ix,ga
-  !       stop
-        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
-
-  goto 10                      !! READ NEXT LEVEL OR PARAMETER
-  !
-  ! CLOSING OF INPUT DATA FILE
-  !
-30   call pbclose(lunit,ierr)     !! FINISHED READING / CLOSING GRIB FILE
-
-  nuvz=iumax
-  nwz =iwmax
-  if(nuvz.eq.nlev_ec) nwz=nlev_ec+1
-
-  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
-
-  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
-  !********************
-
-  write(*,*)
-  write(*,*)
-  write(*,'(a,2i7)') '# of vertical levels in ECMWF data: ', &
-       nuvz+1,nwz
-  write(*,*)
-  write(*,'(a)') 'Mother domain:'
-  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
-       xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
-  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
-       ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
-  write(*,*)
-
-
-  ! 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 i=1,nwz
-    j=numskip+i
-    k=nlev_ec+1+numskip+i
-    akm(nwz-i+1)=zsec2(10+j)
-    bkm(nwz-i+1)=zsec2(10+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/gridcheck_fnl.f90
===================================================================
--- src/gridcheck_fnl.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ 	(revision )
@@ -1,545 +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: 17/11/2005, Caroline Forster, GFS data          *
-  !             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                *
-  ! VHEIGHT(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
-  real(kind=4) :: xaux1,xaux2,yaux1,yaux2
-  real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
-  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
-  !HSO  end
-  integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
-  real :: sizesouth,sizenorth,xauxa,pint
-  real :: akm_usort(nwzmax)
-  real,parameter :: eps=0.0001
-
-  ! NCEP GFS
-  real :: pres(nwzmax), help
-
-  integer :: i179,i180,i181
-
-  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
-
-  integer :: isec1(8),isec2(3)
-  real(kind=4) :: zsec4(jpunp)
-  character(len=1) :: opt
-
-  !HSO  grib api error messages
-  character(len=24) :: gribErrorMsg = 'Error reading grib file'
-  character(len=20) :: gribFunction = 'gridcheckwind_gfs'
-  !
-  if (numbnests.ge.1) then
-  write(*,*) ' ###########################################'
-  write(*,*) ' FLEXPART ERROR SUBROUTINE GRIDCHECK:'
-  write(*,*) ' NO NESTED WINDFIELDAS ALLOWED FOR GFS!      '
-  write(*,*) ' ###########################################'
-  stop
-  endif
-
-  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
-
-  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
-
-  !read the grib1 identifiers
-  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'level',isec1(8),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-
-  !get the size and data of the values array
-  call grib_get_real4_array(igrib,'values',zsec4,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-
-  else ! 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,'scaledValueOfFirstFixedSurface', &
-       valSurf,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-
-  !convert to grib1 identifiers
-  isec1(6)=-1
-  isec1(7)=-1
-  isec1(8)=-1
-  if ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
-    isec1(6)=33          ! indicatorOfParameter
-    isec1(7)=100         ! indicatorOfTypeOfLevel
-    isec1(8)=valSurf/100 ! level, convert to hPa
-  elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
-    isec1(6)=7           ! indicatorOfParameter
-    isec1(7)=1           ! indicatorOfTypeOfLevel
-    isec1(8)=0
-  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
-       .and.(discipl.eq.2)) then ! LSM
-    isec1(6)=81          ! indicatorOfParameter
-    isec1(7)=1           ! indicatorOfTypeOfLevel
-    isec1(8)=0
-  endif
-
-  if (isec1(6).ne.-1) then
-  !  get the size and data of the values array
-    call grib_get_real4_array(igrib,'values',zsec4,iret)
-!    call grib_check(iret,gribFunction,gribErrorMsg)
-  endif
-
-  endif ! gribVer
-
-  if(ifield.eq.1) then
-
-  !get the required fields from section 2
-  !store compatible to gribex input
-  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_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
-  xaux2=xaux2in+360.    !!! Transform FNL xaux2in from -1 to 359 (Xuekue Fang, 31 Jan 2013)
-  yaux1=yaux1in
-  yaux2=yaux2in
-  write(*,*) 'xaux1,xaux2,yaux1,yaux2',xaux1,xaux2,yaux1,yaux2
-  nxfield=isec2(2)
-  ny=isec2(3)
-
-  
-
-  if((abs(xaux1).lt.eps).and.(xaux2.ge.359)) then ! NCEP DATA FROM 0 TO
-    xaux1=-179.0                             ! 359 DEG EAST ->
-    xaux2=-179.0+360.-360./real(nxfield)    ! TRANSFORMED TO -179
-  endif                                      ! TO 180 DEG EAST
-  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.
-  xlon0=xaux1
-  ylat0=yaux1
-  write(*,*) 'xlon0,ylat0',xlon0,ylat0
-
-  dx=(xaux2-xaux1)/real(nxfield-1)
-  dy=(yaux2-yaux1)/real(ny-1)
-  dxconst=180./(dx*r_earth*pi)
-  dyconst=180./(dy*r_earth*pi)
-  write(*,*) 'dx,dy,dxconst,dyconst',dx,dy,dxconst,dyconst
-  !HSO end edits
-
-
-  ! 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
-  endif ! ifield.eq.1
-
-  if (nxshift.lt.0) stop 'nxshift (par_mod) must not be negative'
-  if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
-
-  ! NCEP ISOBARIC LEVELS
-  !*********************
-
-  if((isec1(6).eq.33).and.(isec1(7).eq.100)) then ! check for U wind
-    iumax=iumax+1
-    pres(iumax)=real(isec1(8))*100.0
-  endif
-
-
-  i179=nint(179./dx)
-  if (dx.lt.0.7) then
-    i180=nint(180./dx)+1    ! 0.5 deg data
-  else
-    i180=nint(179./dx)+1    ! 1 deg data
-  endif
-  i181=i180+1
-
-
-  ! NCEP TERRAIN
-  !*************
-
-  if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
-    do jy=0,ny-1
-      do ix=0,nxfield-1
-        help=zsec4(nxfield*(ny-jy-1)+ix+1)
-        if(ix.le.i180) then
-          oro(i179+ix,jy)=help
-          excessoro(i179+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
-        else
-          oro(ix-i181,jy)=help
-          excessoro(ix-i181,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
-        endif
-      end do
-    end do
-  endif
-
-  ! NCEP LAND SEA MASK
-  !*******************
-
-  if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
-    do jy=0,ny-1
-      do ix=0,nxfield-1
-        help=zsec4(nxfield*(ny-jy-1)+ix+1)
-        if(ix.le.i180) then
-          lsm(i179+ix,jy)=help
-        else
-          lsm(ix-i181,jy)=help
-        endif
-      end do
-    end do
-  endif
-
-  call grib_release(igrib)
-
-  goto 10                      !! READ NEXT LEVEL OR PARAMETER
-  !
-  ! CLOSING OF INPUT DATA FILE
-  !
-
-  ! HSO
-30   continue
-  call grib_close_file(ifile)
-  ! HSO end edits
-
-  nuvz=iumax
-  nwz =iumax
-  nlev_ec=iumax
-
-  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
-
-  if (nuvz.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,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
-  !********************
-
-  write(*,*)
-  write(*,*)
-  write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', &
-       nuvz,nwz
-  write(*,*)
-  write(*,'(a)') 'other domain:'
-  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
-       xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
-  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
-       ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
-  write(*,*)
-
-
-  ! 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 i=1,nwz
-    j=numskip+i
-    k=nlev_ec+1+numskip+i
-    akm_usort(nwz-i+1)=pres(nwz-i+1)
-    bkm(nwz-i+1)=0.0
-  end do
-
-  !******************************
-  ! change Sabine Eckhardt: akm should always be in descending order ... readwind adapted!
-  !******************************
-      do i=1,nwz
-         if (akm_usort(1).gt.akm_usort(2)) then
-            akm(i)=akm_usort(i)
-         else
-            akm(i)=akm_usort(nwz-i+1)
-         endif
-      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
-  !*****************************************************************************
-
-  do i=1,nuvz
-     akz(i)=akm(i)
-     bkz(i)=bkm(i)
-  end do
-
-  ! 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/gridcheck_gfs_emos.f90
===================================================================
--- src/gridcheck_gfs_emos.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
+++ 	(revision )
@@ -1,425 +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: 17/11/2005, Caroline Forster, GFS data          *
-  !                                                                     *
-  !**********************************************************************
-  !                                                                     *
-  ! 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 par_mod
-  use com_mod
-  use conv_mod
-
-  implicit none
-
-  integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
-  real :: xaux1,xaux2,yaux1,yaux2,sizesouth,sizenorth,xauxa,pint
-  real :: akm_usort(nwzmax)
-
-  ! NCEP GFS
-  real :: pres(nwzmax), help
-
-  ! 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 :: isec0(2),isec1(56),isec2(22+nxmax+nymax),isec3(2)
-  integer :: isec4(64),inbuff(jpack),ilen,iswap,ierr,iword,lunit
-  real :: zsec2(60+2*nuvzmax),zsec3(2),zsec4(jpunp)
-  character(len=1) :: opt, yoper = 'D'
-
-
-  iumax=0
-  iwmax=0
-
-  if(ideltas.gt.0) then
-    ifn=1
-  else
-    ifn=numbwf
-  endif
-  !
-  ! OPENING OF DATA FILE (GRIB CODE)
-  !
-5   call pbopen(lunit,path(3)(1:length(3))//wfname(ifn),'r',ierr)
-  if(ierr.lt.0) goto 999
-
-  ifield=0
-10   ifield=ifield+1
-  !
-  ! GET NEXT FIELDS
-  !
-  call pbgrib(lunit,inbuff,jpack,ilen,ierr)
-  if(ierr.eq.-1) goto 30    ! EOF DETECTED
-  if(ierr.lt.-1) goto 999   ! ERROR DETECTED
-
-  ierr=1
-
-  ! Check whether we are on a little endian or on a big endian computer
-  !********************************************************************
-
-  !if (inbuff(1).eq.1112101447) then         ! little endian, swap bytes
-  !  iswap=1+ilen/4
-  !  call swap32(inbuff,iswap)
-  !else if (inbuff(1).ne.1196575042) then    ! big endian
-  !  stop 'subroutine gridcheck: corrupt GRIB data'
-  !endif
-
-  call gribex(isec0,isec1,isec2,zsec2,isec3,zsec3,isec4, &
-       zsec4,jpunp,inbuff,jpack,iword,yoper,ierr)
-  if (ierr.ne.0) goto 10   ! ERROR DETECTED
-
-  if(ifield.eq.1) then
-    nxfield=isec2(2)
-    ny=isec2(3)
-    xaux1=real(isec2(5))/1000.
-    xaux2=real(isec2(8))/1000.
-    if((xaux1.eq.0.).and.(xaux2.eq.-1.0)) then ! NCEP DATA FROM 0 TO
-      xaux1=-179.0                             ! 359 DEG EAST ->
-      xaux2=180.0                              ! TRANSFORMED TO -179
-    endif                                      ! TO 180 DEG EAST
-    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.
-    yaux1=real(isec2(7))/1000.
-    yaux2=real(isec2(4))/1000.
-    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)
-
-
-  ! 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
-  endif
-
-  if (nxshift.lt.0) stop 'nxshift (par_mod) must not be negative'
-  if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
-
-  ! NCEP ISOBARIC LEVELS
-  !*********************
-
-  if((isec1(6).eq.33).and.(isec1(7).eq.100)) then
-    iumax=iumax+1
-    pres(iumax)=real(isec1(8))*100.0
-  endif
-
-  ! NCEP TERRAIN
-  !*************
-
-  if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
-    do jy=0,ny-1
-      do ix=0,nxfield-1
-        help=zsec4(nxfield*(ny-jy-1)+ix+1)
-        if(ix.le.180) then
-          oro(179+ix,jy)=help
-          excessoro(179+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
-        else
-          oro(ix-181,jy)=help
-          excessoro(ix-181,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
-        endif
-      end do
-    end do
-  endif
-
-  ! NCEP LAND SEA MASK
-  !*******************
-
-  if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
-    do jy=0,ny-1
-      do ix=0,nxfield-1
-        help=zsec4(nxfield*(ny-jy-1)+ix+1)
-        if(ix.le.180) then
-          lsm(179+ix,jy)=help
-        else
-          lsm(ix-181,jy)=help
-        endif
-      end do
-    end do
-  endif
-
-  goto 10                      !! READ NEXT LEVEL OR PARAMETER
-  !
-  ! CLOSING OF INPUT DATA FILE
-  !
-30   call pbclose(lunit,ierr)     !! FINISHED READING / CLOSING GRIB FILE
-
-  nuvz=iumax
-  nwz =iumax
-  nlev_ec=iumax
-
-  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
-
-  if (nuvz.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,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
-  !********************
-
-  write(*,*)
-  write(*,*)
-  write(*,'(a,2i7)') '# of vertical levels in ECMWF data: ', &
-       nuvz,nwz
-  write(*,*)
-  write(*,'(a)') 'Mother domain:'
-  write(*,'(a,f10.1,a1,f10.1,a,f10.1)') '  Longitude range: ', &
-       xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
-  write(*,'(a,f10.1,a1,f10.1,a,f10.1)') '  Latitude range: ', &
-       ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
-  write(*,*)
-
-
-  ! 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 i=1,nwz
-    j=numskip+i
-    k=nlev_ec+1+numskip+i
-    akm_usort(nwz-i+1)=pres(nwz-i+1)
-    bkm(nwz-i+1)=0.0
-  end do
-
-  !******************************
-  ! change Sabine Eckhardt: akm should always be in descending order ... readwind adapted!
-  !******************************
-      do i=1,nwz
-         if (akm_usort(1).gt.akm_usort(2)) then
-            akm(i)=akm_usort(i)
-         else
-            akm(i)=akm_usort(nwz-i+1)
-         endif
-      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
-  !*****************************************************************************
-
-  do i=1,nuvz
-     akz(i)=akm(i)
-     bkz(i)=bkm(i)
-  end do
-
-  ! 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/gridcheck_nests_emos.f90
===================================================================
--- src/gridcheck_nests_emos.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
+++ 	(revision )
@@ -1,301 +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_nests
-
-  !*****************************************************************************
-  !                                                                            *
-  !     This routine checks the grid specification for the nested model        *
-  !     domains. It is similar to subroutine gridcheck, which checks the       *
-  !     mother domain.                                                         *
-  !                                                                            *
-  !     Authors: A. Stohl, G. Wotawa                                           *
-  !                                                                            *
-  !     8 February 1999                                                        *
-  !                                                                            *
-  !*****************************************************************************
-
-  use par_mod
-  use com_mod
-
-  implicit none
-
-  integer :: i,j,k,l,ifn,ifield,iumax,iwmax,numskip,nlev_ecn
-  integer :: nuvzn,nwzn
-  real :: akmn(nwzmax),bkmn(nwzmax),akzn(nuvzmax),bkzn(nuvzmax)
-  real :: xaux1,xaux2,yaux1,yaux2
-
-  ! 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 :: isec0(2),isec1(56),isec2(22+nxmaxn+nymaxn),isec3(2)
-  integer :: isec4(64),inbuff(jpack),ilen,ierr,iword,lunit
-  !integer iswap
-  real :: zsec2(60+2*nuvzmax),zsec3(2),zsec4(jpunp)
-  character(len=1) :: yoper = 'D'
-
-  xresoln(0)=1.       ! resolution enhancement for mother grid
-  yresoln(0)=1.       ! resolution enhancement for mother grid
-
-  ! Loop about all nesting levels
-  !******************************
-
-  do l=1,numbnests
-
-    iumax=0
-    iwmax=0
-
-    if(ideltas.gt.0) then
-      ifn=1
-    else
-      ifn=numbwf
-    endif
-  !
-  ! OPENING OF DATA FILE (GRIB CODE)
-  !
-5   call pbopen(lunit,path(numpath+2*(l-1)+1) &
-         (1:length(numpath+2*(l-1)+1))//wfnamen(l,ifn),'r',ierr)
-    if(ierr.lt.0) goto 999
-
-    ifield=0
-10    ifield=ifield+1
-  !
-  ! GET NEXT FIELDS
-  !
-      call pbgrib(lunit,inbuff,jpack,ilen,ierr)
-      if(ierr.eq.-1) goto 30    ! EOF DETECTED
-      if(ierr.lt.-1) goto 999   ! ERROR DETECTED
-
-      ierr=1
-
-  ! Check whether we are on a little endian or on a big endian computer
-  !********************************************************************
-
-  !    if (inbuff(1).eq.1112101447) then         ! little endian, swap bytes
-  !      iswap=1+ilen/4
-  !      call swap32(inbuff,iswap)
-  !    else if (inbuff(1).ne.1196575042) then    ! big endian
-  !      stop 'subroutine gridcheck: corrupt GRIB data'
-  !    endif
-
-      call gribex(isec0,isec1,isec2,zsec2,isec3,zsec3,isec4, &
-           zsec4,jpunp,inbuff,jpack,iword,yoper,ierr)
-      if (ierr.ne.0) goto 999   ! ERROR DETECTED
-
-      if(ifield.eq.1) then
-        nxn(l)=isec2(2)
-        nyn(l)=isec2(3)
-        xaux1=real(isec2(5))/1000.
-        xaux2=real(isec2(8))/1000.
-        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.
-        yaux1=real(isec2(7))/1000.
-        yaux2=real(isec2(4))/1000.
-        xlon0n(l)=xaux1
-        ylat0n(l)=yaux1
-        dxn(l)=(xaux2-xaux1)/real(nxn(l)-1)
-        dyn(l)=(yaux2-yaux1)/real(nyn(l)-1)
-        nlev_ecn=isec2(12)/2-1
-      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 j=0,nyn(l)-1
-        do i=0,nxn(l)-1
-          oron(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
-        end do
-      end do
-      endif
-      if(isec1(6).eq.172) then
-      do j=0,nyn(l)-1
-        do i=0,nxn(l)-1
-          lsmn(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
-        end do
-      end do
-      endif
-      if(isec1(6).eq.160) then
-      do j=0,nyn(l)-1
-        do i=0,nxn(l)-1
-          excessoron(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
-        end do
-      end do
-      endif
-
-      goto 10                 !! READ NEXT LEVEL OR PARAMETER
-  !
-  ! CLOSING OF INPUT DATA FILE
-  !
-30   call pbclose(lunit,ierr)     !! FINISHED READING / CLOSING GRIB FILE
-
-    nuvzn=iumax
-    nwzn=iwmax
-    if(nuvzn.eq.nlev_ec) nwzn=nlev_ecn+1
-
-    if (nxn(l).gt.nxmaxn) then
-  write(*,*) 'FLEXPART error: Too many grid points in x direction.'
-  write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'
-  write(*,*) 'for nesting level ',l
-  write(*,*) 'Or change parameter settings in file par_mod.'
-  write(*,*) nxn(l),nxmaxn
-      stop
-    endif
-
-    if (nyn(l).gt.nymaxn) then
-  write(*,*) 'FLEXPART error: Too many grid points in y direction.'
-  write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'
-  write(*,*) 'for nesting level ',l
-  write(*,*) 'Or change parameter settings in file par_mod.'
-  write(*,*) nyn(l),nymaxn
-      stop
-    endif
-
-    if ((nuvzn.gt.nuvzmax).or.(nwzn.gt.nwzmax)) then
-  write(*,*) 'FLEXPART error: Nested wind fields have too many'// &
-       'vertical levels.'
-  write(*,*) 'Problem was encountered for nesting level ',l
-      stop
-    endif
-
-
-  ! Output of grid info
-  !********************
-
-  write(*,'(a,i2)') 'Nested domain #: ',l
-  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
-       xlon0n(l),' to ',xlon0n(l)+(nxn(l)-1)*dxn(l), &
-       '   Grid distance: ',dxn(l)
-  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
-       ylat0n(l),' to ',ylat0n(l)+(nyn(l)-1)*dyn(l), &
-       '   Grid distance: ',dyn(l)
-  write(*,*)
-
-  ! Determine, how much the resolutions in the nests are enhanced as
-  ! compared to the mother grid
-  !*****************************************************************
-
-    xresoln(l)=dx/dxn(l)
-    yresoln(l)=dy/dyn(l)
-
-  ! Determine the mother grid coordinates of the corner points of the
-  ! nested grids
-  ! Convert first to geographical coordinates, then to grid coordinates
-  !********************************************************************
-
-    xaux1=xlon0n(l)
-    xaux2=xlon0n(l)+real(nxn(l)-1)*dxn(l)
-    yaux1=ylat0n(l)
-    yaux2=ylat0n(l)+real(nyn(l)-1)*dyn(l)
-
-    xln(l)=(xaux1-xlon0)/dx
-    xrn(l)=(xaux2-xlon0)/dx
-    yln(l)=(yaux1-ylat0)/dy
-    yrn(l)=(yaux2-ylat0)/dy
-
-
-    if ((xln(l).lt.0.).or.(yln(l).lt.0.).or. &
-         (xrn(l).gt.real(nxmin1)).or.(yrn(l).gt.real(nymin1))) then
-      write(*,*) 'Nested domain does not fit into mother domain'
-      write(*,*) 'For global mother domain fields, you can shift'
-      write(*,*) 'shift the mother domain into x-direction'
-      write(*,*) 'by setting nxshift (file par_mod) to a'
-      write(*,*) 'positive value. Execution is terminated.'
-      stop
-    endif
-
-
-  ! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
-  ! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
-
-    numskip=nlev_ecn-nuvzn ! number of ecmwf model layers not used by FLEXPART
-    do i=1,nwzn
-      j=numskip+i
-      k=nlev_ecn+1+numskip+i
-      akmn(nwzn-i+1)=zsec2(10+j)
-      bkmn(nwzn-i+1)=zsec2(10+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
-  !*****************************************************************************
-
-    akzn(1)=0.
-    bkzn(1)=1.0
-    do i=1,nuvzn
-      akzn(i+1)=0.5*(akmn(i+1)+akmn(i))
-      bkzn(i+1)=0.5*(bkmn(i+1)+bkmn(i))
-    end do
-    nuvzn=nuvzn+1
-
-  ! Check, whether the heights of the model levels of the nested
-  ! wind fields are consistent with those of the mother domain.
-  ! If not, terminate model run.
-  !*************************************************************
-
-    do i=1,nuvz
-      if ((akzn(i).ne.akz(i)).or.(bkzn(i).ne.bkz(i))) then
-  write(*,*) 'FLEXPART error: The wind fields of nesting level',l
-  write(*,*) 'are not consistent with the mother domain:'
-  write(*,*) 'Differences in vertical levels detected.'
-        stop
-      endif
-    end do
-
-    do i=1,nwz
-      if ((akmn(i).ne.akm(i)).or.(bkmn(i).ne.bkm(i))) then
-  write(*,*) 'FLEXPART error: The wind fields of nesting level',l
-  write(*,*) 'are not consistent with the mother domain:'
-  write(*,*) 'Differences in vertical levels detected.'
-        stop
-      endif
-    end do
-
-  end do
-
-  return
-
-999   write(*,*)
-  write(*,*) ' ###########################################'// &
-       '###### '
-  write(*,*) '                FLEXPART SUBROUTINE GRIDCHECK:'
-  write(*,*) ' CAN NOT OPEN INPUT DATA FILE '//wfnamen(l,ifn)
-  write(*,*) ' FOR NESTING LEVEL ',k
-  write(*,*) ' ###########################################'// &
-       '###### '
-  stop
-
-end subroutine gridcheck_nests
Index: c/gridcheck_orig_ecmwf.f90
===================================================================
--- src/gridcheck_orig_ecmwf.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
+++ 	(revision )
@@ -1,560 +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
-  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
-
-  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 (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
-
-  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
-  !********************
-
-  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(*,*)
-
-
-  ! 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: src/makefile
===================================================================
--- src/makefile	(revision d8eed026ebec594e9a91de6c6fbc45b354a0e4b2)
+++ src/makefile	(revision de4c5e9885c62eece05be19e1c0165ade923e201)
@@ -64,7 +64,7 @@
 	MPIF90    = /usr/bin/mpif90.openmpi
 
-	INCPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/include
+	INCPATH1 = /xnilu_wrk/projects/FLEXPART/flex_wrk/bin64/grib_api/include
 	INCPATH2 = /usr/include
-	LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
+	LIBPATH1 = /xnilu_wrk/projects/FLEXPART/flex_wrk/bin64/grib_api/lib
 endif
 
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision 08a38b55d3749abb3d1df461e12c98d18b329a51)
+++ src/par_mod.f90	(revision de4c5e9885c62eece05be19e1c0165ade923e201)
@@ -58,5 +58,5 @@
   !****************************************************************
 
-  logical, parameter :: lnokernel=.false.
+  logical, parameter :: lnokernel=.true.
 
   !*********************************************************************
@@ -218,6 +218,6 @@
   !**************************************************
 
-  integer,parameter :: maxpart=3000000
-  integer,parameter :: maxspec=4
+  integer,parameter :: maxpart=300000000
+  integer,parameter :: maxspec=1
 
   real,parameter :: minmass=0.0001
Index: c/readlanduse_int1.f90
===================================================================
--- src/readlanduse_int1.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ 	(revision )
@@ -1,161 +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 readlanduse
-
-  !*****************************************************************************
-  !                                                                            *
-  !      Reads the landuse inventory into memory and relates it to Leaf Area   *
-  !      Index and roughness length.                                           *
-  !                                                                            *
-  !      AUTHOR: Andreas Stohl, 10 January 1994                                *
-  !                                                                            *
-  !*****************************************************************************
-  !                                                                            *
-  ! Variables:                                                                 *
-  ! i                       loop indices                                       *
-  ! landinvent(1200,600,13) area fractions of 13 landuse categories            *
-  ! LENGTH(numpath)         length of the path names                           *
-  ! PATH(numpath)           contains the path names                            *
-  ! unitland                unit connected with landuse inventory              *
-  !                                                                            *
-  ! -----                                                                      *
-  ! Sabine Eckhardt, Dec 06 - new landuse inventary                            *
-  ! after                                                                      *
-  ! Belward, A.S., Estes, J.E., and Kline, K.D., 1999,                         *
-  ! The IGBP-DIS 1-Km Land-Cover Data Set DISCover:                            *
-  ! A Project Overview: Photogrammetric Engineering and Remote Sensing,        *
-  ! v. 65, no. 9, p. 1013-1020                                                 *
-  !                                                                            *
-  ! LANDUSE CATEGORIES:                                                        *
-  !                                                                            *
-  ! 1	Urban land                                                             *
-  ! 2	Agricultural land					               *
-  ! 3	Range land					                       *
-  ! 4	Deciduous forest			                               *
-  ! 5	Coniferous forest		                                       *
-  ! 6	Mixed forest including wetland                                         *
-  ! 7	water, both salt and fresh                                             *
-  ! 8	barren land mostly desert                                              *
-  ! 9	nonforested wetland	                                               *
-  ! 10 	mixed agricultural and range land		                       *
-  ! 11 	rocky open areas with low growing shrubs                               *
-  ! 12  ice                                                                    *
-  ! 13  rainforest                                                             *
-  !                                                                            *
-  !*****************************************************************************
-
-  use par_mod
-  use com_mod
-
-  implicit none
-
-  integer :: ix,jy,i,k,lu_cat,lu_perc
-  integer(kind=1) :: ilr
-  integer(kind=1) :: ilr_buffer(2160000)
-  integer :: il,irecread
-  real :: rlr, r2lr
-
-
-  ! Read landuse inventory
-  !***********************
-  ! The landuse information is saved in a compressed format and written
-  ! out by records of the length of 1 BYTE. Each grid cell consists of 3
-  ! Bytes, which include 3 landuse categories (val 1-13 and 16 percentage
-  ! categories) So one half byte is used to store the Landusecat the other
-  ! for the percentageclass in 6.25 steps (100/6.25=16)
-  ! e.g.
-  ! 4 3  percentage 4 = 4*6.25 => 25% landuse class 3
-  ! 2 1  percentage 2 = 2*6.25 => 13% landuse class 1
-  ! 1 12 percentage 1 = 1*6.26 => 6.25% landuse class 12
-
-
-  write (*,*) 'reading: ',path(1)(1:length(1))
-  open(unitland,file=path(1)(1:length(1)) &
-       //'IGBP_int1.dat',status='old', &
-       form='UNFORMATTED', err=998)
-  read (unitland) (ilr_buffer(i),i=1,2160000)
-  close(unitland)
-  write (*,*) 'reading: '
-
-  irecread=1
-  do ix=1,1200
-    do jy=1,600
-  ! the 3 most abundant landuse categories in the inventory
-  ! first half byte contains the landuse class
-  ! second half byte contains the respective percentage
-      do k=1,3
-  ! 1 byte is read
-        ilr=ilr_buffer(irecread)
-        irecread=irecread+1
-  ! as only signed integer values exist an unsigned value is constructed
-        if (ilr.lt.0) then
-           il=ilr+256
-        else
-           il=ilr
-        endif
-  ! dividing by 16 has the effect to get rid of the right half of the byte
-  ! so just the left half remains, this corresponds to a shift right of 4
-  ! bits
-        rlr=real(il)/16.
-        lu_cat=int(rlr)
-  ! the left half of the byte is substracted from the whole in order to
-  ! get only the right half of the byte
-        r2lr=rlr-int(rlr)
-  ! shift left by 4
-        lu_perc=r2lr*16.
-        landinvent(ix,jy,k)=lu_cat
-        landinvent(ix,jy,k+3)=lu_perc
-        if ((jy.lt.10).and.(ix.lt.10)) then
-            write(*,*) 'reading: ', ix, jy, lu_cat, lu_perc
-        endif
-      end do
-    end do
-  end do
-
-  ! Read relation landuse,z0
-  !*****************************
-
-  open(unitsurfdata,file=path(1)(1:length(1))//'surfdata.t', &
-       status='old',err=999)
-
-  do i=1,4
-    read(unitsurfdata,*)
-  end do
-  do i=1,numclass
-    read(unitsurfdata,'(45x,f15.3)') z0(i)
-  end do
-  close(unitsurfdata)
-
-  return
-
-  ! Issue error messages
-  !*********************
-
-998   write(*,*) ' #### FLEXPART ERROR! FILE CONTAINING          ####'
-  write(*,*) ' #### LANDUSE INVENTORY DOES NOT EXIST         ####'
-  stop
-
-999   write(*,*) ' #### FLEXPART ERROR! FILE CONTAINING          ####'
-  write(*,*) ' #### RELATION LANDUSE,z0 DOES NOT EXIST       ####'
-  stop
-
-end subroutine readlanduse
Index: c/readwind.f90
===================================================================
--- src/readwind.f90	(revision b0434e1968915cd3d4a8afcbf2cc275acecb9ff7)
+++ 	(revision )
@@ -1,503 +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 readwind(indj,n,uuh,vvh,wwh)
-
-!**********************************************************************
-!                                                                     *
-!             TRAJECTORY MODEL SUBROUTINE READWIND                    *
-!                                                                     *
-!**********************************************************************
-!                                                                     *
-!             AUTHOR:      G. WOTAWA                                  *
-!             DATE:        1997-08-05                                 *
-!             LAST UPDATE: 2000-10-17, Andreas 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                      *
-!                                                                     *
-!**********************************************************************
-!  Changes, Bernd C. Krueger, Feb. 2001:
-!   Variables tth and qvh (on eta coordinates) in common block
-!**********************************************************************
-!                                                                     *
-! DESCRIPTION:                                                        *
-!                                                                     *
-! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE   *
-! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE          *
-!                                                                     *
-! INPUT:                                                              *
-! indj               indicates number of the wind field to be read in *
-! n                  temporal index for meteorological fields (1 to 3)*
-!                                                                     *
-! IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
-!                                                                     *
-! wfname             File name of data to be read in                  *
-! nx,ny,nuvz,nwz     expected field dimensions                        *
-! nlev_ec            number of vertical levels ecmwf model            *
-! uu,vv,ww           wind fields                                      *
-! tt,qv              temperature and specific humidity                *
-! ps                 surface pressure                                 *
-!                                                                     *
-!**********************************************************************
-
-  use grib_api
-  use par_mod
-  use com_mod
-
-  implicit none
-
-!  include 'grib_api.h'
-
-!HSO  parameters for grib_api
-  integer :: ifile
-  integer :: iret
-  integer :: igrib
-  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
-  integer :: gotGrid
-!HSO  end
-
-  real(sp) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
-  real(sp) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
-  real(sp) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
-  integer :: indj,i,j,k,n,levdiff2,ifield,iumax,iwmax
-
-! 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(sp) :: zsec4(jpunp)
-  real(sp) :: xaux,yaux
-  real(dp) :: xauxin,yauxin
-  real(sp),parameter :: eps=1.e-4
-  real(sp) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
-  real(sp) :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
-
-  logical :: hflswitch,strswitch!,readcloud
-
-!HSO  grib api error messages
-  character(len=24) :: gribErrorMsg = 'Error reading grib file'
-  character(len=20) :: gribFunction = 'readwind'
-
-  hflswitch=.false.
-  strswitch=.false.
-!ZHG test the grib fields that have lcwc without using them
-!  readcloud=.false.
-
-  levdiff2=nlev_ec-nwz+1
-  iumax=0
-  iwmax=0
-
-!
-! OPENING OF DATA FILE (GRIB CODE)
-!
-5 call grib_open_file(ifile,path(3)(1:length(3)) &
-       //trim(wfname(indj)),'r',iret)
-  if (iret.ne.GRIB_SUCCESS) then
-    goto 888   ! 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 50    ! EOF DETECTED
-  elseif (iret.ne.GRIB_SUCCESS) then
-    goto 888   ! 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
-
-    conversion_factor=1.
-
-  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
-    conversion_factor=1.
-    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
-! ESO Cloud water is in a) fields CLWC and CIWC, *or* b) field QC 
-    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
-! 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
-      conversion_factor=1000.
-    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
-      conversion_factor=1000.
-    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 --wrong
-    elseif ((parCat.eq.2).and.(parNum.eq.38) .or. parId .eq. 180) then ! EWSS --correct
-      isec1(6)=180         ! indicatorOfParameter
-!    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS --wrong
-    elseif ((parCat.eq.2).and.(parNum.eq.37) .or. parId .eq. 181) then ! NSSS --correct
-      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*,'***WARNING: 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
-
-!HSO  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
-
-!HSO  get the required fields from section 2 in a gribex compatible manner
-  if (ifield.eq.1) then
-    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_int(igrib,'numberOfVerticalCoordinateValues',isec2(12))
-    call grib_check(iret,gribFunction,gribErrorMsg)
-! CHECK GRID SPECIFICATIONS
-    if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
-    if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
-    if(isec2(12)/2-1.ne.nlev_ec) &
-         stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
-  endif ! ifield
-
-!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,'longitudeOfFirstGridPointInDegrees', &
-         xauxin,iret)
-    call grib_check(iret,gribFunction,gribErrorMsg)
-    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
-         yauxin,iret)
-    call grib_check(iret,gribFunction,gribErrorMsg)
-    if (xauxin.gt.180.) xauxin=xauxin-360.0
-    if (xauxin.lt.-180.) xauxin=xauxin+360.0
-
-    xaux=xauxin+real(nxshift)*dx
-    yaux=yauxin
-    if (xaux.gt.180.) xaux=xaux-360.0
-    if(abs(xaux-xlon0).gt.eps) &
-         stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
-    if(abs(yaux-ylat0).gt.eps) &
-         stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
-    gotGrid=1
-  endif ! gotGrid
-
-  do j=0,nymin1
-    do i=0,nxfield-1
-      k=isec1(8)
-      if(isec1(6).eq.130) tth(i,j,nlev_ec-k+2,n)= &!! TEMPERATURE
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.131) uuh(i,j,nlev_ec-k+2)= &!! U VELOCITY
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.132) vvh(i,j,nlev_ec-k+2)= &!! V VELOCITY
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.133) then                      !! SPEC. HUMIDITY
-        qvh(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
-        if (qvh(i,j,nlev_ec-k+2,n) .lt. 0.) &
-             qvh(i,j,nlev_ec-k+2,n) = 0.
-!        this is necessary because the gridded data may contain
-!        spurious negative values
-      endif
-      if(isec1(6).eq.134) ps(i,j,1,n)= &!! SURF. PRESS.
-           zsec4(nxfield*(ny-j-1)+i+1)
-
-      if(isec1(6).eq.135) wwh(i,j,nlev_ec-k+1)= &!! W VELOCITY
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.141) sd(i,j,1,n)= &!! SNOW DEPTH
-           zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
-      if(isec1(6).eq.151) msl(i,j,1,n)= &!! SEA LEVEL PRESS.
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.164) tcc(i,j,1,n)= &!! CLOUD COVER
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.165) u10(i,j,1,n)= &!! 10 M U VELOCITY
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.166) v10(i,j,1,n)= &!! 10 M V VELOCITY
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.167) tt2(i,j,1,n)= &!! 2 M TEMPERATURE
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.168) td2(i,j,1,n)= &!! 2 M DEW POINT
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.142) then                      !! LARGE SCALE PREC.
-        lsprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
-        if (lsprec(i,j,1,n).lt.0.) lsprec(i,j,1,n)=0.
-      endif
-      if(isec1(6).eq.143) then                      !! CONVECTIVE PREC.
-        convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
-        if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
-      endif
-      if(isec1(6).eq.146) sshf(i,j,1,n)= &!! SENS. HEAT FLUX
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if((isec1(6).eq.146).and.(zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) &
-           hflswitch=.true.    ! Heat flux available
-      if(isec1(6).eq.176) then                      !! SOLAR RADIATION
-        ssr(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
-        if (ssr(i,j,1,n).lt.0.) ssr(i,j,1,n)=0.
-      endif
-      if(isec1(6).eq.180) ewss(i,j)= &!! EW SURFACE STRESS
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.181) nsss(i,j)= &!! NS SURFACE STRESS
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(((isec1(6).eq.180).or.(isec1(6).eq.181)).and. &
-           (zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) strswitch=.true.    ! stress available
-!sec        strswitch=.true.
-      if(isec1(6).eq.129) oro(i,j)= &!! ECMWF OROGRAPHY
-           zsec4(nxfield*(ny-j-1)+i+1)/ga
-      if(isec1(6).eq.160) excessoro(i,j)= &!! STANDARD DEVIATION OF OROGRAPHY
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.172) lsm(i,j)= &!! ECMWF LAND SEA MASK
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
-      if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
-!ZHG READING CLOUD FIELDS ASWELL
-! ESO TODO: add check for if one of clwc/ciwc missing (error),
-! also if all 3 cw fields present, use qc and disregard the others
-      if(isec1(6).eq.246) then  !! CLWC  Cloud liquid water content [kg/kg]
-        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
-        readclouds=.true.
-        sumclouds=.false.
-      endif
-      if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
-        ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
-      endif
-!ZHG end
-!ESO read qc (=clwc+ciwc)
-      if(isec1(6).eq.201031) then  !! QC  Cloud liquid water content [kg/kg]
-        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
-        readclouds=.true.
-        sumclouds=.true.
-      endif
-
-    end do
-  end do
-
-  call grib_release(igrib)
-  goto 10                      !! READ NEXT LEVEL OR PARAMETER
-!
-! CLOSING OF INPUT DATA FILE
-!
-
-50 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
-
-  if(levdiff2.eq.0) then
-    iwmax=nlev_ec+1
-    do i=0,nxmin1
-      do j=0,nymin1
-        wwh(i,j,nlev_ec+1)=0.
-      end do
-    end do
-  endif
-
-! For global fields, assign the leftmost data column also to the rightmost
-! data column; if required, shift whole grid by nxshift grid points
-!*************************************************************************
-
-  if (xglobal) then
-    call shift_field_0(ewss,nxfield,ny)
-    call shift_field_0(nsss,nxfield,ny)
-    call shift_field_0(oro,nxfield,ny)
-    call shift_field_0(excessoro,nxfield,ny)
-    call shift_field_0(lsm,nxfield,ny)
-    call shift_field(ps,nxfield,ny,1,1,2,n)
-    call shift_field(sd,nxfield,ny,1,1,2,n)
-    call shift_field(msl,nxfield,ny,1,1,2,n)
-    call shift_field(tcc,nxfield,ny,1,1,2,n)
-    call shift_field(u10,nxfield,ny,1,1,2,n)
-    call shift_field(v10,nxfield,ny,1,1,2,n)
-    call shift_field(tt2,nxfield,ny,1,1,2,n)
-    call shift_field(td2,nxfield,ny,1,1,2,n)
-    call shift_field(lsprec,nxfield,ny,1,1,2,n)
-    call shift_field(convprec,nxfield,ny,1,1,2,n)
-    call shift_field(sshf,nxfield,ny,1,1,2,n)
-    call shift_field(ssr,nxfield,ny,1,1,2,n)
-    call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
-    call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
-    call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
-    call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
-    call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
-!ZHG
-    call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,2,n)
-    if (.not.sumclouds) call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,2,n)
-!ZHG end
-
-  endif
-
-  do i=0,nxmin1
-    do j=0,nymin1
-      surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
-    end do
-  end do
-
-  if ((.not.hflswitch).or.(.not.strswitch)) then
-    write(*,*) 'WARNING: No flux data contained in GRIB file ', &
-         wfname(indj)
-
-! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
-! As ECMWF has increased the model resolution, such that now the first model
-! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
-! (3rd model level in FLEXPART) for the profile method
-!***************************************************************************
-
-    do i=0,nxmin1
-      do j=0,nymin1
-        plev1=akz(3)+bkz(3)*ps(i,j,1,n)
-        pmean=0.5*(ps(i,j,1,n)+plev1)
-        tv=tth(i,j,3,n)*(1.+0.61*qvh(i,j,3,n))
-        fu=-r_air*tv/ga/pmean
-        hlev1=fu*(plev1-ps(i,j,1,n))   ! HEIGTH OF FIRST MODEL LAYER
-        ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
-        fflev1=sqrt(uuh(i,j,3)**2+vvh(i,j,3)**2)
-        call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
-             tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, &
-             surfstr(i,j,1,n),sshf(i,j,1,n))
-        if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
-        if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
-      end do
-    end do
-  endif
-
-
-! Assign 10 m wind to model level at eta=1.0 to have one additional model
-! level at the ground
-! Specific humidity is taken the same as at one level above
-! Temperature is taken as 2 m temperature
-!**************************************************************************
-
-  do i=0,nxmin1
-    do j=0,nymin1
-      uuh(i,j,1)=u10(i,j,1,n)
-      vvh(i,j,1)=v10(i,j,1,n)
-      qvh(i,j,1,n)=qvh(i,j,2,n)
-      tth(i,j,1,n)=tt2(i,j,1,n)
-    end do
-  end do
-
-  if(iumax.ne.nuvz-1) stop 'READWIND: NUVZ NOT CONSISTENT'
-  if(iwmax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'
-
-  return
-
-888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
-  write(*,*) ' #### ',wfname(indj),'                    #### '
-  write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
-  stop 'Execution terminated'
-
-end subroutine readwind
-
Index: c/readwind_gfs_emos.f90
===================================================================
--- src/readwind_gfs_emos.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
+++ 	(revision )
@@ -1,533 +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 readwind(indj,n,uuh,vvh,wwh)
-
-  !**********************************************************************
-  !                                                                     *
-  !             TRAJECTORY MODEL SUBROUTINE READWIND                    *
-  !                                                                     *
-  !**********************************************************************
-  !                                                                     *
-  !             AUTHOR:      G. WOTAWA                                  *
-  !             DATE:        1997-08-05                                 *
-  !             LAST UPDATE: 2000-10-17, Andreas Stohl                  *
-  !             CAHENGE: 16/11/2005, Caroline Forster, GFS data         *
-  !                                                                     *
-  !**********************************************************************
-  !  Changes, Bernd C. Krueger, Feb. 2001:
-  !   Variables tth and qvh (on eta coordinates) in common block
-  !**********************************************************************
-  !                                                                     *
-  ! DESCRIPTION:                                                        *
-  !                                                                     *
-  ! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE   *
-  ! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE          *
-  !                                                                     *
-  ! INPUT:                                                              *
-  ! indj               indicates number of the wind field to be read in *
-  ! n                  temporal index for meteorological fields (1 to 3)*
-  !                                                                     *
-  ! IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
-  !                                                                     *
-  ! wfname             File name of data to be read in                  *
-  ! nx,ny,nuvz,nwz     expected field dimensions                        *
-  ! nlev_ec            number of vertical levels ecmwf model            *
-  ! uu,vv,ww           wind fields                                      *
-  ! tt,qv              temperature and specific humidity                *
-  ! ps                 surface pressure                                 *
-  !                                                                     *
-  !**********************************************************************
-
-  use par_mod
-  use com_mod
-
-  implicit none
-
-  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
-  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
-  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
-  integer :: ii,indj,i,j,k,n,levdiff2,ifield,iumax,iwmax,lunit
-
-  ! NCEP
-
-  integer :: numpt,numpu,numpv,numpw,numprh
-  real :: help, temp, ew
-  real :: elev
-  real :: ulev1(0:nxmax-1,0:nymax-1),vlev1(0:nxmax-1,0:nymax-1)
-  real :: tlev1(0:nxmax-1,0:nymax-1)
-  real :: qvh2(0:nxmax-1,0:nymax-1)
-
-  ! 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 :: isec0(2),isec1(56),isec2(22+nxmax+nymax),isec3(2)
-  integer :: isec4(64),inbuff(jpack),ilen,iswap,ierr,iword
-  real :: zsec2(60+2*nuvzmax),zsec3(2),zsec4(jpunp)
-  real :: xaux,yaux,xaux0,yaux0
-  real,parameter :: eps=1.e-4
-  real :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
-  real :: plev1,hlev1,ff10m,fflev1
-
-  character(len=1) :: yoper = 'D'
-  logical :: hflswitch,strswitch
-
-  hflswitch=.false.
-  strswitch=.false.
-  levdiff2=nlev_ec-nwz+1
-  iumax=0
-  iwmax=0
-
-  !
-  ! OPENING OF DATA FILE (GRIB CODE)
-  !
-5   call pbopen(lunit,path(3)(1:length(3))//wfname(indj),'r',ierr)
-  if(ierr.lt.0) goto 999
-
-  numpt=0
-  numpu=0
-  numpv=0
-  numpw=0
-  numprh=0
-  ifield=0
-10   ifield=ifield+1
-  !
-  ! GET NEXT FIELDS
-  !
-  call pbgrib(lunit,inbuff,jpack,ilen,ierr)
-  if(ierr.eq.-1) goto 50    ! EOF DETECTED
-  if(ierr.lt.-1) goto 888   ! ERROR DETECTED
-
-  ierr=1
-
-  ! Check whether we are on a little endian or on a big endian computer
-  !********************************************************************
-
-  !if (inbuff(1).eq.1112101447) then         ! little endian, swap bytes
-  !  iswap=1+ilen/4
-  !  call swap32(inbuff,iswap)
-  !else if (inbuff(1).ne.1196575042) then    ! big endian
-  !  stop 'subroutine gridcheck: corrupt GRIB data'
-  !endif
-
-
-  call gribex(isec0,isec1,isec2,zsec2,isec3,zsec3,isec4, &
-       zsec4,jpunp,inbuff,jpack,iword,yoper,ierr)
-  if (ierr.ne.0) goto 10   ! ERROR DETECTED
-
-  if(ifield.eq.1) then
-
-  ! CHECK GRID SPECIFICATIONS
-
-    if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
-    if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
-    xaux=real(isec2(5))/1000.+real(nxshift)*dx
-    if(xaux.eq.0.) xaux=-179.0     ! NCEP DATA
-    yaux=real(isec2(7))/1000.
-    xaux0=xlon0
-    yaux0=ylat0
-    if(xaux.lt.0.) xaux=xaux+360.
-    if(yaux.lt.0.) yaux=yaux+360.
-    if(xaux0.lt.0.) xaux0=xaux0+360.
-    if(yaux0.lt.0.) yaux0=yaux0+360.
-    if(abs(xaux-xaux0).gt.eps) &
-         stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
-    if(abs(yaux-yaux0).gt.eps) &
-         stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
-  endif
-
-  do j=0,nymin1
-    do i=0,nxfield-1
-      if((isec1(6).eq.011).and.(isec1(7).eq.100)) then
-  ! TEMPERATURE
-         if((i.eq.0).and.(j.eq.0)) then
-            do ii=1,nuvz
-              if ((isec1(8)*100.0).eq.akz(ii)) numpt=ii
-            end do
-        endif
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          tth(179+i,j,numpt,n)=help
-        else
-          tth(i-181,j,numpt,n)=help
-        endif
-      endif
-      if((isec1(6).eq.033).and.(isec1(7).eq.100)) then
-  ! U VELOCITY
-         if((i.eq.0).and.(j.eq.0)) then
-            do ii=1,nuvz
-              if ((isec1(8)*100.0).eq.akz(ii)) numpu=ii
-            end do
-        endif
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          uuh(179+i,j,numpu)=help
-        else
-          uuh(i-181,j,numpu)=help
-        endif
-      endif
-      if((isec1(6).eq.034).and.(isec1(7).eq.100)) then
-  ! V VELOCITY
-         if((i.eq.0).and.(j.eq.0)) then
-            do ii=1,nuvz
-              if ((isec1(8)*100.0).eq.akz(ii)) numpv=ii
-            end do
-        endif
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          vvh(179+i,j,numpv)=help
-        else
-          vvh(i-181,j,numpv)=help
-        endif
-      endif
-      if((isec1(6).eq.052).and.(isec1(7).eq.100)) then
-  ! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
-         if((i.eq.0).and.(j.eq.0)) then
-            do ii=1,nuvz
-              if ((isec1(8)*100.0).eq.akz(ii)) numprh=ii
-            end do
-        endif
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          qvh(179+i,j,numprh,n)=help
-        else
-          qvh(i-181,j,numprh,n)=help
-        endif
-      endif
-      if((isec1(6).eq.001).and.(isec1(7).eq.001)) then
-  ! SURFACE PRESSURE
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          ps(179+i,j,1,n)=help
-        else
-          ps(i-181,j,1,n)=help
-        endif
-      endif
-      if((isec1(6).eq.039).and.(isec1(7).eq.100)) then
-  ! W VELOCITY
-         if((i.eq.0).and.(j.eq.0)) then
-            do ii=1,nuvz
-              if ((isec1(8)*100.0).eq.akz(ii)) numpw=ii
-            end do
-        endif
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          wwh(179+i,j,numpw)=help
-        else
-          wwh(i-181,j,numpw)=help
-        endif
-      endif
-      if((isec1(6).eq.066).and.(isec1(7).eq.001)) then
-  ! SNOW DEPTH
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          sd(179+i,j,1,n)=help
-        else
-          sd(i-181,j,1,n)=help
-        endif
-      endif
-      if((isec1(6).eq.002).and.(isec1(7).eq.102)) then
-  ! MEAN SEA LEVEL PRESSURE
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          msl(179+i,j,1,n)=help
-        else
-          msl(i-181,j,1,n)=help
-        endif
-      endif
-      if((isec1(6).eq.071).and.(isec1(7).eq.244)) then
-  ! TOTAL CLOUD COVER
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          tcc(179+i,j,1,n)=help
-        else
-          tcc(i-181,j,1,n)=help
-        endif
-      endif
-      if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
-           (isec1(8).eq.10)) then
-  ! 10 M U VELOCITY
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          u10(179+i,j,1,n)=help
-        else
-          u10(i-181,j,1,n)=help
-        endif
-      endif
-      if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
-           (isec1(8).eq.10)) then
-  ! 10 M V VELOCITY
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          v10(179+i,j,1,n)=help
-        else
-          v10(i-181,j,1,n)=help
-        endif
-      endif
-      if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
-           (isec1(8).eq.02)) then
-  ! 2 M TEMPERATURE
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          tt2(179+i,j,1,n)=help
-        else
-          tt2(i-181,j,1,n)=help
-        endif
-      endif
-      if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
-           (isec1(8).eq.02)) then
-  ! 2 M DEW POINT TEMPERATURE
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          td2(179+i,j,1,n)=help
-        else
-          td2(i-181,j,1,n)=help
-        endif
-      endif
-      if((isec1(6).eq.062).and.(isec1(7).eq.001)) then
-  ! LARGE SCALE PREC.
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          lsprec(179+i,j,1,n)=help
-        else
-          lsprec(i-181,j,1,n)=help
-        endif
-      endif
-      if((isec1(6).eq.063).and.(isec1(7).eq.001)) then
-  ! CONVECTIVE PREC.
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          convprec(179+i,j,1,n)=help
-        else
-          convprec(i-181,j,1,n)=help
-        endif
-      endif
-  ! SENS. HEAT FLUX
-      sshf(i,j,1,n)=0.0     ! not available from gfs.tccz.pgrbfxx files
-      hflswitch=.false.    ! Heat flux not available
-  ! SOLAR RADIATIVE FLUXES
-      ssr(i,j,1,n)=0.0      ! not available from gfs.tccz.pgrbfxx files
-  ! EW SURFACE STRESS
-      ewss(i,j)=0.0         ! not available from gfs.tccz.pgrbfxx files
-  ! NS SURFACE STRESS
-      nsss(i,j)=0.0         ! not available from gfs.tccz.pgrbfxx files
-      strswitch=.false.    ! stress not available
-      if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
-  ! TOPOGRAPHY
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          oro(179+i,j)=help
-          excessoro(179+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
-        else
-          oro(i-181,j)=help
-          excessoro(i-181,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
-        endif
-      endif
-      if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
-  ! LAND SEA MASK
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          lsm(179+i,j)=help
-        else
-          lsm(i-181,j)=help
-        endif
-      endif
-      if((isec1(6).eq.221).and.(isec1(7).eq.001)) then
-  ! MIXING HEIGHT
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          hmix(179+i,j,1,n)=help
-        else
-          hmix(i-181,j,1,n)=help
-        endif
-      endif
-      if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
-           (isec1(8).eq.02)) then
-  ! 2 M RELATIVE HUMIDITY
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          qvh2(179+i,j)=help
-        else
-          qvh2(i-181,j)=help
-        endif
-      endif
-      if((isec1(6).eq.011).and.(isec1(7).eq.107)) then
-  ! TEMPERATURE LOWEST SIGMA LEVEL
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          tlev1(179+i,j)=help
-        else
-          tlev1(i-181,j)=help
-        endif
-      endif
-      if((isec1(6).eq.033).and.(isec1(7).eq.107)) then
-  ! U VELOCITY LOWEST SIGMA LEVEL
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          ulev1(179+i,j)=help
-        else
-          ulev1(i-181,j)=help
-        endif
-      endif
-      if((isec1(6).eq.034).and.(isec1(7).eq.107)) then
-  ! V VELOCITY LOWEST SIGMA LEVEL
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.180) then
-          vlev1(179+i,j)=help
-        else
-          vlev1(i-181,j)=help
-        endif
-      endif
-
-    end do
-  end do
-
-  if((isec1(6).eq.33).and.(isec1(7).eq.100)) then
-  ! NCEP ISOBARIC LEVELS
-    iumax=iumax+1
-  endif
-
-
-  goto 10                      !! READ NEXT LEVEL OR PARAMETER
-  !
-  ! CLOSING OF INPUT DATA FILE
-  !
-50   call pbclose(lunit,ierr)     !! FINNISHED READING / CLOSING GRIB FILE
-
-  ! TRANSFORM RH TO SPECIFIC HUMIDITY
-
-  do j=0,ny-1
-    do i=0,nxfield-1
-      do k=1,nuvz
-        help=qvh(i,j,k,n)
-        temp=tth(i,j,k,n)
-        plev1=akm(k)+bkm(k)*ps(i,j,1,n)
-        elev=ew(temp)*help/100.0
-        qvh(i,j,k,n)=xmwml*(elev/(plev1-((1.0-xmwml)*elev)))
-      end do
-    end do
-  end do
-
-  ! CALCULATE 2 M DEW POINT FROM 2 M RELATIVE HUMIDITY
-  ! USING BOLTON'S (1980) FORMULA
-  ! BECAUSE td2 IS NOT AVAILABLE FROM NCEP GFS DATA
-
-  do j=0,ny-1
-    do i=0,nxfield-1
-        help=qvh2(i,j)
-        temp=tt2(i,j,1,n)
-        elev=ew(temp)/100.*help/100.   !vapour pressure in hPa
-        td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273.
-        if (help.le.0.) td2(i,j,1,n)=tt2(i,j,1,n)
-    end do
-  end do
-
-  if(levdiff2.eq.0) then
-    iwmax=nlev_ec+1
-    do i=0,nxmin1
-      do j=0,nymin1
-        wwh(i,j,nlev_ec+1)=0.
-      end do
-    end do
-  endif
-
-
-  ! For global fields, assign the leftmost data column also to the rightmost
-  ! data column; if required, shift whole grid by nxshift grid points
-  !*************************************************************************
-
-  if (xglobal) then
-    call shift_field_0(ewss,nxfield,ny)
-    call shift_field_0(nsss,nxfield,ny)
-    call shift_field_0(oro,nxfield,ny)
-    call shift_field_0(excessoro,nxfield,ny)
-    call shift_field_0(lsm,nxfield,ny)
-    call shift_field_0(ulev1,nxfield,ny)
-    call shift_field_0(vlev1,nxfield,ny)
-    call shift_field_0(tlev1,nxfield,ny)
-    call shift_field_0(qvh2,nxfield,ny)
-    call shift_field(ps,nxfield,ny,1,1,2,n)
-    call shift_field(sd,nxfield,ny,1,1,2,n)
-    call shift_field(msl,nxfield,ny,1,1,2,n)
-    call shift_field(tcc,nxfield,ny,1,1,2,n)
-    call shift_field(u10,nxfield,ny,1,1,2,n)
-    call shift_field(v10,nxfield,ny,1,1,2,n)
-    call shift_field(tt2,nxfield,ny,1,1,2,n)
-    call shift_field(td2,nxfield,ny,1,1,2,n)
-    call shift_field(lsprec,nxfield,ny,1,1,2,n)
-    call shift_field(convprec,nxfield,ny,1,1,2,n)
-    call shift_field(sshf,nxfield,ny,1,1,2,n)
-    call shift_field(ssr,nxfield,ny,1,1,2,n)
-    call shift_field(hmix,nxfield,ny,1,1,2,n)
-    call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
-    call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
-    call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
-    call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
-    call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
-  endif
-
-  do i=0,nxmin1
-    do j=0,nymin1
-      surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
-    end do
-  end do
-
-  if ((.not.hflswitch).or.(.not.strswitch)) then
-  !  write(*,*) 'WARNING: No flux data contained in GRIB file ',
-  !    +  wfname(indj)
-
-  ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
-  !***************************************************************************
-
-    do i=0,nxmin1
-      do j=0,nymin1
-        hlev1=30.0                     ! HEIGHT OF FIRST MODEL SIGMA LAYER
-        ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
-        fflev1=sqrt(ulev1(i,j)**2+vlev1(i,j)**2)
-        call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
-             tt2(i,j,1,n),tlev1(i,j),ff10m,fflev1, &
-             surfstr(i,j,1,n),sshf(i,j,1,n))
-        if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
-        if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
-      end do
-    end do
-  endif
-
-
-  if(iumax.ne.nuvz) stop 'READWIND: NUVZ NOT CONSISTENT'
-  if(iumax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'
-
-  return
-888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
-  write(*,*) ' #### ',wfname(indj),'                    #### '
-  write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
-  stop 'Execution terminated'
-999   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
-  write(*,*) ' #### ',wfname(indj),'                    #### '
-  write(*,*) ' #### CANNOT BE OPENED !!!                    #### '
-  stop 'Execution terminated'
-
-end subroutine readwind
Index: c/readwind_mpi.f90
===================================================================
--- src/readwind_mpi.f90	(revision fdc0f03d0fd969385acd98d3a3206e9d802172bb)
+++ 	(revision )
@@ -1,528 +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 readwind(indj,n,uuh,vvh,wwh)
-
-!**********************************************************************
-!                                                                     *
-!             TRAJECTORY MODEL SUBROUTINE READWIND                    *
-!                                                                     *
-!**********************************************************************
-!                                                                     *
-!             AUTHOR:      G. WOTAWA                                  *
-!             DATE:        1997-08-05                                 *
-!             LAST UPDATE: 2000-10-17, Andreas 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                      *
-!                                                                     *
-!**********************************************************************
-!  Changes, Bernd C. Krueger, Feb. 2001:
-!   Variables tth and qvh (on eta coordinates) in common block
-!  eso 2014:
-!    This version for reading windfields in advance by a separate
-!    MPI process.
-!    TODO: can be merged with readwind.f90 if using numwfmem in
-!          shift_field
-!**********************************************************************
-!                                                                     *
-! DESCRIPTION:                                                        *
-!                                                                     *
-! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE   *
-! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE          *
-!                                                                     *
-! INPUT:                                                              *
-! indj               indicates number of the wind field to be read in *
-! n                  temporal index for meteorological fields (1 to 3)*
-!                                                                     *
-! IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
-!                                                                     *
-! wfname             File name of data to be read in                  *
-! nx,ny,nuvz,nwz     expected field dimensions                        *
-! nlev_ec            number of vertical levels ecmwf model            *
-! uu,vv,ww           wind fields                                      *
-! tt,qv              temperature and specific humidity                *
-! ps                 surface pressure                                 *
-!                                                                     *
-!**********************************************************************
-
-  use GRIB_API
-  use par_mod
-  use com_mod
-  use mpi_mod
-
-  implicit none
-
-!  include 'grib_api.h'
-
-!HSO  parameters for grib_api
-  integer :: ifile
-  integer :: iret
-  integer :: igrib
-  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
-  integer :: gotGrid
-!HSO  end
-
-  real(sp), intent(inout) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
-  real(sp), intent(inout) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
-  real(sp), intent(inout) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
-  integer, intent(in) :: indj,n
-  integer :: i,j,k,levdiff2,ifield,iumax,iwmax
-
-! 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) :: zsec4(jpunp)
-  real(kind=4) :: xaux,yaux
-  real(kind=8) :: xauxin,yauxin
-  real,parameter :: eps=1.e-4
-  real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
-  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
-
-  logical :: hflswitch,strswitch !,readcloud
-
-!HSO  grib api error messages
-  character(len=24) :: gribErrorMsg = 'Error reading grib file'
-  character(len=20) :: gribFunction = 'readwind'
-
-! measure time
-  if (mp_measure_time) call mpif_mtime('readwind',0)
-
-  hflswitch=.false.
-  strswitch=.false.
-!ZHG test the grib fields that have lcwc without using them
-!  readcloud=.false.
-!hg end
-  levdiff2=nlev_ec-nwz+1
-  iumax=0
-  iwmax=0
-
-!
-! OPENING OF DATA FILE (GRIB CODE)
-!
-5 call grib_open_file(ifile,path(3)(1:length(3)) &
-       //trim(wfname(indj)),'r',iret)
-  if (iret.ne.GRIB_SUCCESS) then
-    goto 888   ! 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 50    ! EOF DETECTED
-  elseif (iret.ne.GRIB_SUCCESS) then
-    goto 888   ! 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
-
-    conversion_factor=1.
-
-  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
-    conversion_factor=1.
-    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
-! ESO Cloud water is in a) fields CLWC and CIWC, *or* b) field qc 
-    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
-! 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
-      conversion_factor=1000.
-    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
-      conversion_factor=1000.
-    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 --wrong
-    elseif ((parCat.eq.2).and.(parNum.eq.38) .or. parId .eq. 180) then ! EWSS --correct
-      isec1(6)=180         ! indicatorOfParameter
-!    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS --wrong
-    elseif ((parCat.eq.2).and.(parNum.eq.37) .or. parId .eq. 181) then ! NSSS --correct
-      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*,'***WARNING: 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
-
-!HSO  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
-
-!HSO  get the required fields from section 2 in a gribex compatible manner
-  if (ifield.eq.1) then
-    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_int(igrib,'numberOfVerticalCoordinateValues',isec2(12))
-    call grib_check(iret,gribFunction,gribErrorMsg)
-! CHECK GRID SPECIFICATIONS
-    if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
-    if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
-    if(isec2(12)/2-1.ne.nlev_ec) &
-         stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
-  endif ! ifield
-
-!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,'longitudeOfFirstGridPointInDegrees', &
-         xauxin,iret)
-    call grib_check(iret,gribFunction,gribErrorMsg)
-    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
-         yauxin,iret)
-    call grib_check(iret,gribFunction,gribErrorMsg)
-    if (xauxin.gt.180.) xauxin=xauxin-360.0
-    if (xauxin.lt.-180.) xauxin=xauxin+360.0
-
-    xaux=xauxin+real(nxshift)*dx
-    yaux=yauxin
-    if (xaux.gt.180.) xaux=xaux-360.0
-    if(abs(xaux-xlon0).gt.eps) &
-         stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
-    if(abs(yaux-ylat0).gt.eps) &
-         stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
-    gotGrid=1
-  endif ! gotGrid
-
-  do j=0,nymin1
-    do i=0,nxfield-1
-      k=isec1(8)
-      if(isec1(6).eq.130) tth(i,j,nlev_ec-k+2,n)= &!! TEMPERATURE
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.131) uuh(i,j,nlev_ec-k+2)= &!! U VELOCITY
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.132) vvh(i,j,nlev_ec-k+2)= &!! V VELOCITY
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.133) then                      !! SPEC. HUMIDITY
-        qvh(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
-        if (qvh(i,j,nlev_ec-k+2,n) .lt. 0.) &
-             qvh(i,j,nlev_ec-k+2,n) = 0.
-!        this is necessary because the gridded data may contain
-!        spurious negative values
-      endif
-      if(isec1(6).eq.134) ps(i,j,1,n)= &!! SURF. PRESS.
-           zsec4(nxfield*(ny-j-1)+i+1)
-
-      if(isec1(6).eq.135) wwh(i,j,nlev_ec-k+1)= &!! W VELOCITY
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.141) sd(i,j,1,n)= &!! SNOW DEPTH
-           zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
-      if(isec1(6).eq.151) msl(i,j,1,n)= &!! SEA LEVEL PRESS.
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.164) tcc(i,j,1,n)= &!! CLOUD COVER
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.165) u10(i,j,1,n)= &!! 10 M U VELOCITY
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.166) v10(i,j,1,n)= &!! 10 M V VELOCITY
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.167) tt2(i,j,1,n)= &!! 2 M TEMPERATURE
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.168) td2(i,j,1,n)= &!! 2 M DEW POINT
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.142) then                      !! LARGE SCALE PREC.
-        lsprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
-        if (lsprec(i,j,1,n).lt.0.) lsprec(i,j,1,n)=0.
-      endif
-      if(isec1(6).eq.143) then                      !! CONVECTIVE PREC.
-        convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
-        if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
-      endif
-      if(isec1(6).eq.146) sshf(i,j,1,n)= &!! SENS. HEAT FLUX
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if((isec1(6).eq.146).and.(zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) &
-           hflswitch=.true.    ! Heat flux available
-      if(isec1(6).eq.176) then                      !! SOLAR RADIATION
-        ssr(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
-        if (ssr(i,j,1,n).lt.0.) ssr(i,j,1,n)=0.
-      endif
-      if(isec1(6).eq.180) ewss(i,j)= &!! EW SURFACE STRESS
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.181) nsss(i,j)= &!! NS SURFACE STRESS
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(((isec1(6).eq.180).or.(isec1(6).eq.181)).and. &
-           (zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) strswitch=.true.    ! stress available
-!sec        strswitch=.true.
-      if(isec1(6).eq.129) oro(i,j)= &!! ECMWF OROGRAPHY
-           zsec4(nxfield*(ny-j-1)+i+1)/ga
-      if(isec1(6).eq.160) excessoro(i,j)= &!! STANDARD DEVIATION OF OROGRAPHY
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.172) lsm(i,j)= &!! ECMWF LAND SEA MASK
-           zsec4(nxfield*(ny-j-1)+i+1)
-      if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
-      if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
-!ZHG READING CLOUD FIELDS ASWELL
-      if(isec1(6).eq.246) then  !! CLWC  Cloud liquid water content [kg/kg]
-        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
-        readclouds = .true.
-        sumclouds=.false.
-      endif
-      if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
-        ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
-      endif
-!ZHG end
-!ESO read qc (=clwc+ciwc)
-      if(isec1(6).eq.201031) then  !! QC  Cloud liquid water content [kg/kg]
-        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
-        readclouds=.true.
-        sumclouds=.true.
-      endif
-
-    end do
-  end do
-
-  call grib_release(igrib)
-  goto 10                      !! READ NEXT LEVEL OR PARAMETER
-!
-! CLOSING OF INPUT DATA FILE
-!
-
-50 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
-
-  if(levdiff2.eq.0) then
-    iwmax=nlev_ec+1
-    do i=0,nxmin1
-      do j=0,nymin1
-        wwh(i,j,nlev_ec+1)=0.
-      end do
-    end do
-  endif
-
-! For global fields, assign the leftmost data column also to the rightmost
-! data column; if required, shift whole grid by nxshift grid points
-!*************************************************************************
-
-  if (xglobal) then
-    call shift_field_0(ewss,nxfield,ny)
-    call shift_field_0(nsss,nxfield,ny)
-    call shift_field_0(oro,nxfield,ny)
-    call shift_field_0(excessoro,nxfield,ny)
-    call shift_field_0(lsm,nxfield,ny)
-    ! call shift_field(ps,nxfield,ny,1,1,2,n)
-    ! call shift_field(sd,nxfield,ny,1,1,2,n)
-    ! call shift_field(msl,nxfield,ny,1,1,2,n)
-    ! call shift_field(tcc,nxfield,ny,1,1,2,n)
-    ! call shift_field(u10,nxfield,ny,1,1,2,n)
-    ! call shift_field(v10,nxfield,ny,1,1,2,n)
-    ! call shift_field(tt2,nxfield,ny,1,1,2,n)
-    ! call shift_field(td2,nxfield,ny,1,1,2,n)
-    ! call shift_field(lsprec,nxfield,ny,1,1,2,n)
-    ! call shift_field(convprec,nxfield,ny,1,1,2,n)
-    ! call shift_field(sshf,nxfield,ny,1,1,2,n)
-    ! call shift_field(ssr,nxfield,ny,1,1,2,n)
-    ! call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
-    ! call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
-    call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
-    call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
-    call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
-! numwfmem
-    call shift_field(ps,nxfield,ny,1,1,numwfmem,n)
-    call shift_field(sd,nxfield,ny,1,1,numwfmem,n)
-    call shift_field(msl,nxfield,ny,1,1,numwfmem,n)
-    call shift_field(tcc,nxfield,ny,1,1,numwfmem,n)
-    call shift_field(u10,nxfield,ny,1,1,numwfmem,n)
-    call shift_field(v10,nxfield,ny,1,1,numwfmem,n)
-    call shift_field(tt2,nxfield,ny,1,1,numwfmem,n)
-    call shift_field(td2,nxfield,ny,1,1,numwfmem,n)
-    call shift_field(lsprec,nxfield,ny,1,1,numwfmem,n)
-    call shift_field(convprec,nxfield,ny,1,1,numwfmem,n)
-    call shift_field(sshf,nxfield,ny,1,1,numwfmem,n)
-    call shift_field(ssr,nxfield,ny,1,1,numwfmem,n)
-    call shift_field(tth,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
-    call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
-!hg
-    call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
-    if (.not.sumclouds) call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
-!hg end
-
-  endif
-
-  do i=0,nxmin1
-    do j=0,nymin1
-      surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
-    end do
-  end do
-
-  if (mp_measure_time) call mpif_mtime('readwind',1)
-
-  if ((.not.hflswitch).or.(.not.strswitch)) then
-    write(*,*) 'WARNING: No flux data contained in GRIB file ', &
-         wfname(indj)
-
-! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
-! As ECMWF has increased the model resolution, such that now the first model
-! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
-! (3rd model level in FLEXPART) for the profile method
-!***************************************************************************
-
-    do i=0,nxmin1
-      do j=0,nymin1
-        plev1=akz(3)+bkz(3)*ps(i,j,1,n)
-        pmean=0.5*(ps(i,j,1,n)+plev1)
-        tv=tth(i,j,3,n)*(1.+0.61*qvh(i,j,3,n))
-        fu=-r_air*tv/ga/pmean
-        hlev1=fu*(plev1-ps(i,j,1,n))   ! HEIGTH OF FIRST MODEL LAYER
-        ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
-        fflev1=sqrt(uuh(i,j,3)**2+vvh(i,j,3)**2)
-        call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
-             tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, &
-             surfstr(i,j,1,n),sshf(i,j,1,n))
-        if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
-        if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
-      end do
-    end do
-  endif
-
-
-! Assign 10 m wind to model level at eta=1.0 to have one additional model
-! level at the ground
-! Specific humidity is taken the same as at one level above
-! Temperature is taken as 2 m temperature
-!**************************************************************************
-
-  do i=0,nxmin1
-    do j=0,nymin1
-      uuh(i,j,1)=u10(i,j,1,n)
-      vvh(i,j,1)=v10(i,j,1,n)
-      qvh(i,j,1,n)=qvh(i,j,2,n)
-      tth(i,j,1,n)=tt2(i,j,1,n)
-    end do
-  end do
-
-  if(iumax.ne.nuvz-1) stop 'READWIND: NUVZ NOT CONSISTENT'
-  if(iwmax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'
-
-  return
-
-888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
-  write(*,*) ' #### ',wfname(indj),'                    #### '
-  write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
-  stop 'Execution terminated'
-
-end subroutine readwind
-
Index: c/readwind_nests_emos.f90
===================================================================
--- src/readwind_nests_emos.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
+++ 	(revision )
@@ -1,286 +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 readwind_nests(indj,n,uuhn,vvhn,wwhn)
-  !                           i   i  o    o    o
-  !*****************************************************************************
-  !                                                                            *
-  !     This routine reads the wind fields for the nested model domains.       *
-  !     It is similar to subroutine readwind, which reads the mother domain.   *
-  !                                                                            *
-  !     Authors: A. Stohl, G. Wotawa                                           *
-  !                                                                            *
-  !     8 February 1999                                                        *
-  !                                                                            *
-  !     Last update: 17 October 2000, A. Stohl                                 *
-  !                                                                            *
-  !*****************************************************************************
-  !  Changes, Bernd C. Krueger, Feb. 2001:
-  !   Variables tthn and qvhn (on eta coordinates) in common block
-  !*****************************************************************************
-
-  use par_mod
-  use com_mod
-
-  implicit none
-
-  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
-  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
-  real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
-  integer :: indj,i,j,k,n,levdiff2,ifield,iumax,iwmax,lunit,l
-
-  ! 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 :: isec0(2),isec1(56),isec2(22+nxmaxn+nymaxn),isec3(2)
-  integer :: isec4(64),inbuff(jpack),ilen,ierr,iword
-  !integer iswap
-  real :: zsec2(60+2*nuvzmax),zsec3(2),zsec4(jpunp)
-  real :: xaux,yaux,xaux0,yaux0
-  real :: ewss(0:nxmaxn-1,0:nymaxn-1),nsss(0:nxmaxn-1,0:nymaxn-1)
-  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
-
-  character(len=1) :: yoper = 'D'
-  logical :: hflswitch,strswitch
-
-
-  do l=1,numbnests
-    hflswitch=.false.
-    strswitch=.false.
-    levdiff2=nlev_ec-nwz+1
-    iumax=0
-    iwmax=0
-
-  !
-  ! OPENING OF DATA FILE (GRIB CODE)
-  !
-5   call pbopen(lunit,path(numpath+2*(l-1)+1) &
-         (1:length(numpath+2*(l-1)+1))//wfnamen(l,indj),'r',ierr)
-    if(ierr.lt.0) goto 999
-
-    ifield=0
-10    ifield=ifield+1
-  !
-  ! GET NEXT FIELDS
-  !
-    call pbgrib(lunit,inbuff,jpack,ilen,ierr)
-    if(ierr.eq.-1) goto 50    ! EOF DETECTED
-    if(ierr.lt.-1) goto 888   ! ERROR DETECTED
-
-    ierr=1
-
-  ! Check whether we are on a little endian or on a big endian computer
-  !********************************************************************
-
-  !  if (inbuff(1).eq.1112101447) then         ! little endian, swap bytes
-  !    iswap=1+ilen/4
-  !    call swap32(inbuff,iswap)
-  !  else if (inbuff(1).ne.1196575042) then    ! big endian
-  !    stop 'subroutine gridcheck: corrupt GRIB data'
-  !  endif
-
-    call gribex(isec0,isec1,isec2,zsec2,isec3,zsec3,isec4, &
-         zsec4,jpunp,inbuff,jpack,iword,yoper,ierr)
-    if (ierr.ne.0) goto 888   ! ERROR DETECTED
-
-    if(ifield.eq.1) then
-
-  ! CHECK GRID SPECIFICATIONS
-
-      if(isec2(2).ne.nxn(l)) stop &
-           'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
-      if(isec2(3).ne.nyn(l)) stop &
-           'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
-      if(isec2(12)/2-1.ne.nlev_ec) stop 'READWIND: VERTICAL DISCRET&
-           &IZATION NOT CONSISTENT FOR A NESTING LEVEL'
-      xaux=real(isec2(5))/1000.
-      yaux=real(isec2(7))/1000.
-      xaux0=xlon0n(l)
-      yaux0=ylat0n(l)
-      if(xaux.lt.0.) xaux=xaux+360.
-      if(yaux.lt.0.) yaux=yaux+360.
-      if(xaux0.lt.0.) xaux0=xaux0+360.
-      if(yaux0.lt.0.) yaux0=yaux0+360.
-      if(xaux.ne.xaux0) &
-           stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT FOR A NES&
-           &TING LEVEL'
-      if(yaux.ne.yaux0) &
-           stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT FOR A NEST&
-           &ING LEVEL'
-    endif
-
-
-    do j=0,nyn(l)-1
-      do i=0,nxn(l)-1
-        k=isec1(8)
-        if(isec1(6).eq.130) tthn(i,j,nlev_ec-k+2,n,l)= &!! TEMPERATURE
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.131) uuhn(i,j,nlev_ec-k+2,l)= &!! U VELOCITY
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.132) vvhn(i,j,nlev_ec-k+2,l)= &!! V VELOCITY
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.133) then                         !! SPEC. HUMIDITY
-          qvhn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-          if (qvhn(i,j,nlev_ec-k+2,n,l) .lt. 0.) &
-               qvhn(i,j,nlev_ec-k+2,n,l) = 0.
-  !          this is necessary because the gridded data may contain
-  !          spurious negative values
-        endif
-        if(isec1(6).eq.134) psn(i,j,1,n,l)= &!! SURF. PRESS.
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.135) wwhn(i,j,nlev_ec-k+1,l)= &!! W VELOCITY
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.141) sdn(i,j,1,n,l)= &!! SNOW DEPTH
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS.
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.164) tccn(i,j,1,n,l)= &!! CLOUD COVER
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.165) u10n(i,j,1,n,l)= &!! 10 M U VELOCITY
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.166) v10n(i,j,1,n,l)= &!! 10 M V VELOCITY
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.167) tt2n(i,j,1,n,l)= &!! 2 M TEMPERATURE
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.168) td2n(i,j,1,n,l)= &!! 2 M DEW POINT
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.142) then                         !! LARGE SCALE PREC.
-          lsprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-          if (lsprecn(i,j,1,n,l).lt.0.) lsprecn(i,j,1,n,l)=0.
-        endif
-        if(isec1(6).eq.143) then                         !! CONVECTIVE PREC.
-          convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-          if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0.
-        endif
-        if(isec1(6).eq.146) sshfn(i,j,1,n,l)= &!! SENS. HEAT FLUX
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if((isec1(6).eq.146).and. &
-             (zsec4(nxn(l)*(nyn(l)-j-1)+i+1).ne.0.)) hflswitch=.true.    ! Heat flux available
-        if(isec1(6).eq.176) then                         !! SOLAR RADIATION
-          ssrn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-          if (ssrn(i,j,1,n,l).lt.0.) ssrn(i,j,1,n,l)=0.
-        endif
-        if(isec1(6).eq.180) ewss(i,j)= &!! EW SURFACE STRESS
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.181) nsss(i,j)= &!! NS SURFACE STRESS
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(((isec1(6).eq.180).or.(isec1(6).eq.181)).and. &
-             (zsec4(nxn(l)*(nyn(l)-j-1)+i+1).ne.0.)) strswitch=.true.    ! stress available
-        if(isec1(6).eq.129) oron(i,j,l)= &!! ECMWF OROGRAPHY
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
-        if(isec1(6).eq.160) excessoron(i,j,l)= &!! STANDARD DEVIATION OF OROGRAPHY
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.172) lsmn(i,j,l)= &!! ECMWF LAND SEA MASK
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
-        if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
-        if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
-      end do
-    end do
-
-    goto 10                      !! READ NEXT LEVEL OR PARAMETER
-  !
-  ! CLOSING OF INPUT DATA FILE
-  !
-50   call pbclose(lunit,ierr)     !! FINISHED READING / CLOSING GRIB FILE
-
-    if(levdiff2.eq.0) then
-      iwmax=nlev_ec+1
-      do i=0,nxn(l)-1
-        do j=0,nyn(l)-1
-          wwhn(i,j,nlev_ec+1,l)=0.
-        end do
-      end do
-    endif
-
-    do i=0,nxn(l)-1
-      do j=0,nyn(l)-1
-        surfstrn(i,j,1,n,l)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
-      end do
-    end do
-
-    if ((.not.hflswitch).or.(.not.strswitch)) then
-      write(*,*) 'WARNING: No flux data contained in GRIB file ', &
-           wfnamen(l,indj)
-
-  ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
-  ! As ECMWF has increased the model resolution, such that now the first model
-  ! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
-  ! (3rd model level in FLEXPART) for the profile method
-  !***************************************************************************
-
-      do i=0,nxn(l)-1
-        do j=0,nyn(l)-1
-        plev1=akz(3)+bkz(3)*psn(i,j,1,n,l)
-        pmean=0.5*(psn(i,j,1,n,l)+plev1)
-        tv=tthn(i,j,3,n,l)*(1.+0.61*qvhn(i,j,3,n,l))
-        fu=-r_air*tv/ga/pmean
-        hlev1=fu*(plev1-psn(i,j,1,n,l))   ! HEIGTH OF FIRST MODEL LAYER
-        ff10m= sqrt(u10n(i,j,1,n,l)**2+v10n(i,j,1,n,l)**2)
-        fflev1=sqrt(uuhn(i,j,3,l)**2+vvhn(i,j,3,l)**2)
-        call pbl_profile(psn(i,j,1,n,l),td2n(i,j,1,n,l),hlev1, &
-             tt2n(i,j,1,n,l),tthn(i,j,3,n,l),ff10m,fflev1, &
-             surfstrn(i,j,1,n,l),sshfn(i,j,1,n,l))
-        if(sshfn(i,j,1,n,l).gt.200.) sshfn(i,j,1,n,l)=200.
-        if(sshfn(i,j,1,n,l).lt.-400.) sshfn(i,j,1,n,l)=-400.
-        end do
-      end do
-    endif
-
-
-  ! Assign 10 m wind to model level at eta=1.0 to have one additional model
-  ! level at the ground
-  ! Specific humidity is taken the same as at one level above
-  ! Temperature is taken as 2 m temperature
-  !**************************************************************************
-
-    do i=0,nxn(l)-1
-      do j=0,nyn(l)-1
-        uuhn(i,j,1,l)=u10n(i,j,1,n,l)
-        vvhn(i,j,1,l)=v10n(i,j,1,n,l)
-        qvhn(i,j,1,n,l)=qvhn(i,j,2,n,l)
-        tthn(i,j,1,n,l)=tt2n(i,j,1,n,l)
-      end do
-    end do
-
-    if(iumax.ne.nuvz-1) stop &
-         'READWIND: NUVZ NOT CONSISTENT FOR A NESTING LEVEL'
-    if(iwmax.ne.nwz) stop &
-         'READWIND: NWZ NOT CONSISTENT FOR A NESTING LEVEL'
-
-  end do
-
-  return
-888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
-  write(*,*) ' #### ',wfnamen(l,indj),' FOR NESTING LEVEL  #### '
-  write(*,*) ' #### ',l,' IS NOT GRIB FORMAT !!!           #### '
-  stop 'Execution terminated'
-
-
-999   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
-  write(*,*) ' #### ',wfnamen(l,indj),'                    #### '
-  write(*,*) ' #### CANNOT BE OPENED FOR NESTING LEVEL ',l,'####'
-
-end subroutine readwind_nests
Index: c/verttransform.f90
===================================================================
--- src/verttransform.f90	(revision 72d3a5a739d68f5d3970051fd8ec853fea7376e2)
+++ 	(revision )
@@ -1,889 +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 verttransform(n,uuh,vvh,wwh,pvh)
-!                         i  i   i   i   i
-!*****************************************************************************
-!                                                                            *
-!     This subroutine transforms temperature, dew point temperature and      *
-!     wind components from eta to meter coordinates.                         *
-!     The vertical wind component is transformed from Pa/s to m/s using      *
-!     the conversion factor pinmconv.                                        *
-!     In addition, this routine calculates vertical density gradients        *
-!     needed for the parameterization of the turbulent velocities.           *
-!                                                                            *
-!     Author: A. Stohl, G. Wotawa                                            *
-!                                                                            *
-!     12 August 1996                                                         *
-!     Update: 16 January 1998                                                *
-!                                                                            *
-!     Major update: 17 February 1999                                         *
-!     by G. Wotawa                                                           *
-!                                                                            *
-!     - Vertical levels for u, v and w are put together                      *
-!     - Slope correction for vertical velocity: Modification of calculation  *
-!       procedure                                                            *
-!                                                                            *
-!*****************************************************************************
-!  Changes, Bernd C. Krueger, Feb. 2001:
-!   Variables tth and qvh (on eta coordinates) from common block
-!*****************************************************************************
-! Sabine Eckhardt, March 2007
-! added the variable cloud for use with scavenging - descr. in com_mod
-!*****************************************************************************
-! Date: 2017-05-30 modification of a bug in ew. Don Morton (CTBTO project)   *
-!*****************************************************************************
-!                                                                            *
-! Variables:                                                                 *
-! nx,ny,nz                        field dimensions in x,y and z direction    *
-! clouds(0:nxmax,0:nymax,0:nzmax,numwfmem) cloud field for wet deposition    *
-! uu(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in x-direction [m/s]*
-! vv(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in y-direction [m/s]*
-! ww(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in z-direction      *
-!                                          [deltaeta/s]                      *
-! tt(0:nxmax,0:nymax,nzmax,numwfmem)     temperature [K]                     *
-! pv(0:nxmax,0:nymax,nzmax,numwfmem)     potential voriticity (pvu)          *
-! ps(0:nxmax,0:nymax,numwfmem)           surface pressure [Pa]               *
-!                                                                            *
-!*****************************************************************************
-
-  use par_mod
-  use com_mod
-  use cmapf_mod, only: cc2gll
-!  use mpi_mod
-
-  implicit none
-
-  real,intent(in),dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: uuh,vvh,pvh
-  real,intent(in),dimension(0:nxmax-1,0:nymax-1,nwzmax) :: wwh
-
-  real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: rhoh,uvzlev,wzlev
-  real,dimension(0:nxmax-1,0:nymax-1,nzmax) :: pinmconv
-  real,dimension(0:nxmax-1,0:nymax-1) ::  tvold,pold,pint,tv
-  real,dimension(0:nymax-1) :: cosf
-
-  integer,dimension(0:nxmax-1,0:nymax-1) :: rain_cloud_above,idx
-
-  integer :: ix,jy,kz,iz,n,kmin,ix1,jy1,ixp,jyp,ixm,jym,kz_inv
-  real :: f_qvsat,pressure,rh,lsp,convp,cloudh_min,prec
-  real :: ew,dz1,dz2,dz
-  real :: xlon,ylat,xlonr,dzdx,dzdy
-  real :: dzdx1,dzdx2,dzdy1,dzdy2
-  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
-  real,parameter :: const=r_air/ga
-  real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
-
-  logical :: init = .true.
-  logical :: init_w = .false.
-  logical :: init_r = .false.
-
-
-  !ZHG SEP 2014 tests  
-  ! integer :: cloud_ver,cloud_min, cloud_max 
-  ! integer ::teller(5), convpteller=0, lspteller=0
-  ! real :: cloud_col_wat, cloud_water
-  !ZHG 2015 temporary variables for testing
-  ! real :: rcw(0:nxmax-1,0:nymax-1)
-  ! real :: rpc(0:nxmax-1,0:nymax-1)
-  ! character(len=60) :: zhgpath='/xnilu_wrk/flex_wrk/zhg/'
-  ! character(len=60) :: fnameA,fnameB,fnameC,fnameD,fnameE,fnameF,fnameG,fnameH
-  ! CHARACTER(LEN=3)  :: aspec
-  ! integer :: virr=0
-  !real :: tot_cloud_h
-  !real :: dbg_height(nzmax) 
-!ZHG
-
-!*************************************************************************
-! If verttransform is called the first time, initialize heights of the   *
-! z levels in meter. The heights are the heights of model levels, where  *
-! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
-! the vertical resolution in the z system is doubled. As reference point,*
-! the lower left corner of the grid is used.                             *
-! Unlike in the eta system, no difference between heights for u,v and    *
-! heights for w exists.                                                  *
-!*************************************************************************
-
-
-!eso measure CPU time
-!  call mpif_mtime('verttransform',0)
-
-  if (init) then
-
-
-    if (init_r) then
-
-        open(333,file='heights.txt', &
-          form='formatted')
-        do kz=1,nuvz
-            read(333,*) height(kz)
-        end do
-        close(333)
-        write(*,*) 'height read'
-    else
-
-
-! Search for a point with high surface pressure (i.e. not above significant topography)
-! Then, use this point to construct a reference z profile, to be used at all times
-!*****************************************************************************
-
-    do jy=0,nymin1
-      do ix=0,nxmin1
-        if (ps(ix,jy,1,n).gt.100000.) then
-          ixm=ix
-          jym=jy
-          goto 3
-        endif
-      end do
-    end do
-3   continue
-
-
-    tvold(ixm,jym)=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
-         ps(ixm,jym,1,n))
-    pold(ixm,jym)=ps(ixm,jym,1,n)
-    height(1)=0.
-
-    do kz=2,nuvz
-      pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
-      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
-
-      if (abs(tv(ixm,jym)-tvold(ixm,jym)).gt.0.2) then
-        height(kz)= &
-             height(kz-1)+const*log(pold(ixm,jym)/pint(ixm,jym))* &
-             (tv(ixm,jym)-tvold(ixm,jym))/log(tv(ixm,jym)/tvold(ixm,jym))
-      else
-        height(kz)=height(kz-1)+ &
-             const*log(pold(ixm,jym)/pint(ixm,jym))*tv(ixm,jym)
-      endif
-
-      tvold(ixm,jym)=tv(ixm,jym)
-      pold(ixm,jym)=pint(ixm,jym)
-    end do
-
-    if (init_w) then
-        open(333,file='heights.txt', &
-          form='formatted')
-        do kz=1,nuvz
-              write(333,*) height(kz)
-        end do
-        close(333)
-    endif
-
-    endif ! init
-
-! Determine highest levels that can be within PBL
-!************************************************
-
-    do kz=1,nz
-      if (height(kz).gt.hmixmax) then
-        nmixz=kz
-        goto 9
-      endif
-    end do
-9   continue
-
-! Do not repeat initialization of the Cartesian z grid
-!*****************************************************
-
-    init=.false.
-
-!    dbg_height = height
-
-  endif
-
-
-! Loop over the whole grid
-!*************************
-
-
-  do jy=0,nymin1
-    do ix=0,nxmin1
-      tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
-           ps(ix,jy,1,n))
-    enddo
-  enddo
-  pold=ps(:,:,1,n)
-  uvzlev(:,:,1)=0.
-  wzlev(:,:,1)=0.
-  rhoh(:,:,1)=pold/(r_air*tvold)
-
-
-! Compute heights of eta levels
-!******************************
-
-  do kz=2,nuvz
-    pint=akz(kz)+bkz(kz)*ps(:,:,1,n)
-    tv=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n))
-    rhoh(:,:,kz)=pint(:,:)/(r_air*tv)
-
-    where (abs(tv-tvold).gt.0.2) 
-      uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)* &
-           (tv-tvold)/log(tv/tvold)
-    elsewhere
-      uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)*tv
-    endwhere
-
-    tvold=tv
-    pold=pint
-  end do
-
-
-  do kz=2,nwz-1
-    wzlev(:,:,kz)=(uvzlev(:,:,kz+1)+uvzlev(:,:,kz))/2.
-  end do
-  wzlev(:,:,nwz)=wzlev(:,:,nwz-1)+ &
-       uvzlev(:,:,nuvz)-uvzlev(:,:,nuvz-1)
-
-! pinmconv=(h2-h1)/(p2-p1)
-
-  pinmconv(:,:,1)=(uvzlev(:,:,2))/ &
-       ((aknew(2)+bknew(2)*ps(:,:,1,n))- &
-       (aknew(1)+bknew(1)*ps(:,:,1,n)))
-  do kz=2,nz-1
-    pinmconv(:,:,kz)=(uvzlev(:,:,kz+1)-uvzlev(:,:,kz-1))/ &
-         ((aknew(kz+1)+bknew(kz+1)*ps(:,:,1,n))- &
-         (aknew(kz-1)+bknew(kz-1)*ps(:,:,1,n)))
-  end do
-  pinmconv(:,:,nz)=(uvzlev(:,:,nz)-uvzlev(:,:,nz-1))/ &
-       ((aknew(nz)+bknew(nz)*ps(:,:,1,n))- &
-       (aknew(nz-1)+bknew(nz-1)*ps(:,:,1,n)))
-
-! Levels, where u,v,t and q are given
-!************************************
-
-
-  uu(:,:,1,n)=uuh(:,:,1)
-  vv(:,:,1,n)=vvh(:,:,1)
-  tt(:,:,1,n)=tth(:,:,1,n)
-  qv(:,:,1,n)=qvh(:,:,1,n)
-!hg adding the cloud water 
-  if (readclouds) then
-    clwc(:,:,1,n)=clwch(:,:,1,n)
-    if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n)
-  end if
-!hg 
-  pv(:,:,1,n)=pvh(:,:,1)
-  rho(:,:,1,n)=rhoh(:,:,1)
-
-  uu(:,:,nz,n)=uuh(:,:,nuvz)
-  vv(:,:,nz,n)=vvh(:,:,nuvz)
-  tt(:,:,nz,n)=tth(:,:,nuvz,n)
-  qv(:,:,nz,n)=qvh(:,:,nuvz,n)
-!hg adding the cloud water
-  if (readclouds) then
-    clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
-    if (.not.sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
-  end if
-!hg
-  pv(:,:,nz,n)=pvh(:,:,nuvz)
-  rho(:,:,nz,n)=rhoh(:,:,nuvz)
-
-
-  kmin=2
-  idx=kmin
-  do iz=2,nz-1
-    do jy=0,nymin1
-      do ix=0,nxmin1
-        if(height(iz).gt.uvzlev(ix,jy,nuvz)) then
-          uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
-          vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
-          tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
-          qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
-!hg adding the cloud water
-          if (readclouds) then
-            clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
-            if (.not.sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
-          end if
-!hg
-          pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
-          rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
-        else
-          innuvz: do kz=idx(ix,jy),nuvz
-            if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
-                 (height(iz).le.uvzlev(ix,jy,kz))) then
-              idx(ix,jy)=kz
-              exit innuvz
-            endif
-          enddo innuvz
-        endif
-      enddo
-    enddo
-    do jy=0,nymin1
-      do ix=0,nxmin1
-        if(height(iz).le.uvzlev(ix,jy,nuvz)) then
-          kz=idx(ix,jy)
-          dz1=height(iz)-uvzlev(ix,jy,kz-1)
-          dz2=uvzlev(ix,jy,kz)-height(iz)
-          dz=dz1+dz2
-          uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
-          vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
-          tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
-               +tth(ix,jy,kz,n)*dz1)/dz
-          qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
-               +qvh(ix,jy,kz,n)*dz1)/dz
-!hg adding the cloud water
-          if (readclouds) then
-            clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
-            if (.not.sumclouds) &
-                 &ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
-          end if
-!hg
-          pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
-          rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)*dz2+rhoh(ix,jy,kz)*dz1)/dz
-        endif
-      enddo
-    enddo
-  enddo
-
-
-! Levels, where w is given
-!*************************
-
-  ww(:,:,1,n)=wwh(:,:,1)*pinmconv(:,:,1)
-  ww(:,:,nz,n)=wwh(:,:,nwz)*pinmconv(:,:,nz)
-  kmin=2
-  idx=kmin
-  do iz=2,nz
-    do jy=0,nymin1
-      do ix=0,nxmin1
-        inn:         do kz=idx(ix,jy),nwz
-          if(idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1).and. &
-               height(iz).le.wzlev(ix,jy,kz)) then
-            idx(ix,jy)=kz
-            exit inn
-          endif
-        enddo inn
-      enddo
-    enddo
-    do jy=0,nymin1
-      do ix=0,nxmin1
-        kz=idx(ix,jy)
-        dz1=height(iz)-wzlev(ix,jy,kz-1)
-        dz2=wzlev(ix,jy,kz)-height(iz)
-        dz=dz1+dz2
-        ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(ix,jy,kz-1)*dz2 &
-             +wwh(ix,jy,kz)*pinmconv(ix,jy,kz)*dz1)/dz
-      enddo
-    enddo
-  enddo
-
-! Compute density gradients at intermediate levels
-!*************************************************
-
-  drhodz(:,:,1,n)=(rho(:,:,2,n)-rho(:,:,1,n))/ &
-       (height(2)-height(1))
-  do kz=2,nz-1
-    drhodz(:,:,kz,n)=(rho(:,:,kz+1,n)-rho(:,:,kz-1,n))/ &
-         (height(kz+1)-height(kz-1))
-  end do
-  drhodz(:,:,nz,n)=drhodz(:,:,nz-1,n)
-
-!    end do
-!  end do
-
-
-!****************************************************************
-! Compute slope of eta levels in windward direction and resulting
-! vertical wind correction
-!****************************************************************
-
-  do jy=1,ny-2
-    cosf(jy)=1./cos((real(jy)*dy+ylat0)*pi180)
-  enddo
-
-  kmin=2
-  idx=kmin
-  do iz=2,nz-1
-    do jy=1,ny-2
-      do ix=1,nx-2
-
-        inneta: do kz=idx(ix,jy),nz
-          if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
-               (height(iz).le.uvzlev(ix,jy,kz))) then
-            idx(ix,jy)=kz
-            exit inneta
-          endif
-        enddo inneta
-      enddo
-    enddo
-
-    do jy=1,ny-2
-      do ix=1,nx-2
-        kz=idx(ix,jy)
-        dz1=height(iz)-uvzlev(ix,jy,kz-1)
-        dz2=uvzlev(ix,jy,kz)-height(iz)
-        dz=dz1+dz2
-        ix1=ix-1
-        jy1=jy-1
-        ixp=ix+1
-        jyp=jy+1
-
-        dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2.
-        dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2.
-        dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
-
-        dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2.
-        dzdy2=(uvzlev(ix,jyp,kz)-uvzlev(ix,jy1,kz))/2.
-        dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
-
-        ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*uu(ix,jy,iz,n)*dxconst*cosf(jy)+dzdy*vv(ix,jy,iz,n)*dyconst)
-
-      end do
-
-    end do
-  end do
-
-! If north pole is in the domain, calculate wind velocities in polar
-! stereographic coordinates
-!*******************************************************************
-
-  if (nglobal) then
-    do iz=1,nz
-      do jy=int(switchnorthg)-2,nymin1
-        ylat=ylat0+real(jy)*dy
-        do ix=0,nxmin1
-          xlon=xlon0+real(ix)*dx
-          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
-               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
-               vvpol(ix,jy,iz,n))
-        end do
-      end do
-    end do
-
-
-    do iz=1,nz
-
-! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
-!
-!   AMSnauffer Nov 18 2004 Added check for case vv=0
-!
-      xlon=xlon0+real(nx/2-1)*dx
-      xlonr=xlon*pi/180.
-      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
-           vv(nx/2-1,nymin1,iz,n)**2)
-      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
-        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
-             vv(nx/2-1,nymin1,iz,n))-xlonr
-      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
-        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
-             vv(nx/2-1,nymin1,iz,n))-xlonr
-      else
-        ddpol=pi/2-xlonr
-      endif
-      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
-      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
-
-! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
-      xlon=180.0
-      xlonr=xlon*pi/180.
-      ylat=90.0
-      uuaux=-ffpol*sin(xlonr+ddpol)
-      vvaux=-ffpol*cos(xlonr+ddpol)
-      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
-           vvpolaux)
-
-      jy=nymin1
-      do ix=0,nxmin1
-        uupol(ix,jy,iz,n)=uupolaux
-        vvpol(ix,jy,iz,n)=vvpolaux
-      end do
-    end do
-
-
-! Fix: Set W at pole to the zonally averaged W of the next equator-
-! ward parallel of latitude
-
-    do iz=1,nz
-      wdummy=0.
-      jy=ny-2
-      do ix=0,nxmin1
-        wdummy=wdummy+ww(ix,jy,iz,n)
-      end do
-      wdummy=wdummy/real(nx)
-      jy=nymin1
-      do ix=0,nxmin1
-        ww(ix,jy,iz,n)=wdummy
-      end do
-    end do
-
-  endif
-
-
-! If south pole is in the domain, calculate wind velocities in polar
-! stereographic coordinates
-!*******************************************************************
-
-  if (sglobal) then
-    do iz=1,nz
-      do jy=0,int(switchsouthg)+3
-        ylat=ylat0+real(jy)*dy
-        do ix=0,nxmin1
-          xlon=xlon0+real(ix)*dx
-          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
-               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
-               vvpol(ix,jy,iz,n))
-        end do
-      end do
-    end do
-
-    do iz=1,nz
-
-! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
-!
-!   AMSnauffer Nov 18 2004 Added check for case vv=0
-!
-      xlon=xlon0+real(nx/2-1)*dx
-      xlonr=xlon*pi/180.
-      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
-           vv(nx/2-1,0,iz,n)**2)
-      if (vv(nx/2-1,0,iz,n).lt.0.) then
-        ddpol=atan(uu(nx/2-1,0,iz,n)/ &
-             vv(nx/2-1,0,iz,n))+xlonr
-      else if (vv(nx/2-1,0,iz,n).gt.0.) then
-        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
-             vv(nx/2-1,0,iz,n))+xlonr
-      else
-        ddpol=pi/2-xlonr
-      endif
-      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
-      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
-
-! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
-      xlon=180.0
-      xlonr=xlon*pi/180.
-      ylat=-90.0
-      uuaux=+ffpol*sin(xlonr-ddpol)
-      vvaux=-ffpol*cos(xlonr-ddpol)
-      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
-           vvpolaux)
-
-      jy=0
-      do ix=0,nxmin1
-        uupol(ix,jy,iz,n)=uupolaux
-        vvpol(ix,jy,iz,n)=vvpolaux
-      end do
-    end do
-
-
-! Fix: Set W at pole to the zonally averaged W of the next equator-
-! ward parallel of latitude
-
-    do iz=1,nz
-      wdummy=0.
-      jy=1
-      do ix=0,nxmin1
-        wdummy=wdummy+ww(ix,jy,iz,n)
-      end do
-      wdummy=wdummy/real(nx)
-      jy=0
-      do ix=0,nxmin1
-        ww(ix,jy,iz,n)=wdummy
-      end do
-    end do
-  endif
-
-
-!***********************************************************************************  
-  if (readclouds) then !HG METHOD
-! The method is loops all grids vertically and constructs the 3D matrix for clouds
-! Cloud top and cloud bottom gid cells are assigned as well as the total column 
-! cloud water. For precipitating grids, the type and whether it is in or below 
-! cloud scavenging are assigned with numbers 2-5 (following the old metod).
-! Distinction is done for lsp and convp though they are treated the same in regards
-! to scavenging. Also clouds that are not precipitating are defined which may be 
-! to include future cloud processing by non-precipitating-clouds. 
-!***********************************************************************************
-    write(*,*) 'Global ECMWF fields: using cloud water'
-    clw(:,:,:,n)=0.0
-!    icloud_stats(:,:,:,n)=0.0
-    ctwc(:,:,n)=0.0
-    clouds(:,:,:,n)=0
-! If water/ice are read separately into clwc and ciwc, store sum in clwc
-    if (.not.sumclouds) then 
-      clwc(:,:,:,n) = clwc(:,:,:,n) + ciwc(:,:,:,n)
-    end if
-    do jy=0,nymin1
-      do ix=0,nxmin1
-        lsp=lsprec(ix,jy,1,n)
-        convp=convprec(ix,jy,1,n)
-        prec=lsp+convp
-!        tot_cloud_h=0
-! Find clouds in the vertical
-        do kz=1, nz-1 !go from top to bottom
-          if (clwc(ix,jy,kz,n).gt.0) then      
-! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3 
-            clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
-!            tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz)) 
-            
-!            icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]
-            ctwc(ix,jy,n) = ctwc(ix,jy,n)+clw(ix,jy,kz,n)
-!            icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz))                     ! Cloud BOT height stats      [m]
-            cloudh_min=min(height(kz+1),height(kz))
-!ZHG 2015 extra for testing
-!            clh(ix,jy,kz,n)=height(kz+1)-height(kz)
-!            icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
-!            icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz))          ! Cloud TOP height            [m]
-!ZHG
-          endif
-        end do
-
-! If Precipitation. Define removal type in the vertical
-        if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
-
-          do kz=nz,1,-1 !go Bottom up!
-            if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
-              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1) 
-              clouds(ix,jy,kz,n)=1                               ! is a cloud
-              if (lsp.ge.convp) then
-                clouds(ix,jy,kz,n)=3                            ! lsp in-cloud
-              else
-                clouds(ix,jy,kz,n)=2                             ! convp in-cloud
-              endif                                              ! convective or large scale
-            elseif((clw(ix,jy,kz,n).le.0) .and. (cloudh_min.ge.height(kz))) then ! is below cloud
-              if (lsp.ge.convp) then
-                clouds(ix,jy,kz,n)=5                             ! lsp dominated washout
-              else
-                clouds(ix,jy,kz,n)=4                             ! convp dominated washout
-              endif                                              ! convective or large scale 
-            endif
-
-            if (height(kz).ge. 19000) then                        ! set a max height for removal
-              clouds(ix,jy,kz,n)=0
-            endif !clw>0
-          end do !nz
-        endif ! precipitation
-      end do
-    end do
-
-! eso: copy the relevant data to clw4 to reduce amount of communicated data for MPI
-!    ctwc(:,:,n) = icloud_stats(:,:,4,n)
-
-!**************************************************************************
-  else       ! use old definitions
-!**************************************************************************
-!   create a cloud and rainout/washout field, clouds occur where rh>80%
-!   total cloudheight is stored at level 0
-    write(*,*) 'Global fields: using cloud water from Parameterization'
-    do jy=0,nymin1
-      do ix=0,nxmin1
-! OLD METHOD
-        rain_cloud_above(ix,jy)=0
-        lsp=lsprec(ix,jy,1,n)
-        convp=convprec(ix,jy,1,n)
-        cloudsh(ix,jy,n)=0
-        do kz_inv=1,nz-1
-          kz=nz-kz_inv+1
-          pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
-          rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
-          clouds(ix,jy,kz,n)=0
-          if (rh.gt.0.8) then ! in cloud
-            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
-              rain_cloud_above(ix,jy)=1
-              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
-                   height(kz)-height(kz-1)
-              if (lsp.ge.convp) then
-                clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
-              else
-                clouds(ix,jy,kz,n)=2 ! convp dominated rainout
-              endif
-            else ! no precipitation
-              clouds(ix,jy,kz,n)=1 ! cloud
-            endif
-          else ! no cloud
-            if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
-              if (lsp.ge.convp) then
-                clouds(ix,jy,kz,n)=5 ! lsp dominated washout
-              else
-                clouds(ix,jy,kz,n)=4 ! convp dominated washout
-              endif
-            endif
-          endif
-        end do
-!END OLD METHOD
-      end do
-    end do
-  endif !readclouds
-
-
-     !********* TEST ***************
-     ! WRITE OUT SOME TEST VARIABLES
-     !********* TEST ************'**
-!teller(:)=0
-!virr=virr+1
-!WRITE(aspec, '(i3.3)'), virr
-
-!if (readclouds) then
-!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement.txt'
-!else
-!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement_old.txt'
-!endif
-!
-!OPEN(UNIT=118, FILE=fnameH,FORM='FORMATTED',STATUS = 'UNKNOWN')
-!do kz_inv=1,nz-1
-!  kz=nz-kz_inv+1
-!  !kz=91
-!  do jy=0,nymin1
-!     do ix=0,nxmin1
-!          if (clouds(ix,jy,kz,n).eq.1) teller(1)=teller(1)+1 ! no precipitation cloud
-!          if (clouds(ix,jy,kz,n).eq.2) teller(2)=teller(2)+1 ! convp dominated rainout
-!          if (clouds(ix,jy,kz,n).eq.3) teller(3)=teller(3)+1 ! lsp dominated rainout
-!          if (clouds(ix,jy,kz,n).eq.4) teller(4)=teller(4)+1 ! convp dominated washout
-!          if (clouds(ix,jy,kz,n).eq.5) teller(5)=teller(5)+1 ! lsp dominated washout
-!          
-!        !  write(*,*) height(kz),teller
-!     end do
-!  end do
-!  write(118,*) height(kz),teller
-!  teller(:)=0
-!end do
-!teller(:)=0
-!write(*,*) teller 
-!write(*,*) aspec
-!
-!fnameA=trim(zhgpath)//trim(aspec)//'cloudV.txt'
-!fnameB=trim(zhgpath)//trim(aspec)//'cloudT.txt'
-!fnameC=trim(zhgpath)//trim(aspec)//'cloudB.txt'
-!fnameD=trim(zhgpath)//trim(aspec)//'cloudW.txt'
-!fnameE=trim(zhgpath)//trim(aspec)//'old_cloudV.txt'
-!fnameF=trim(zhgpath)//trim(aspec)//'lsp.txt'
-!fnameG=trim(zhgpath)//trim(aspec)//'convp.txt'
-!if (readclouds) then
-!OPEN(UNIT=111, FILE=fnameA,FORM='FORMATTED',STATUS = 'UNKNOWN')
-!OPEN(UNIT=112, FILE=fnameB,FORM='FORMATTED',STATUS = 'UNKNOWN')
-!OPEN(UNIT=113, FILE=fnameC,FORM='FORMATTED',STATUS = 'UNKNOWN')
-!OPEN(UNIT=114, FILE=fnameD,FORM='FORMATTED',STATUS = 'UNKNOWN')
-!else
-!OPEN(UNIT=115, FILE=fnameE,FORM='FORMATTED',STATUS = 'UNKNOWN')
-!OPEN(UNIT=116, FILE=fnameF,FORM='FORMATTED',STATUS = 'UNKNOWN')
-!OPEN(UNIT=117, FILE=fnameG,FORM='FORMATTED',STATUS = 'UNKNOWN')
-!endif
-!
-!do ix=0,nxmin1
-!if (readclouds) then
-!write(111,*) (icloud_stats(ix,jy,1,n),jy=0,nymin1)
-!write(112,*) (icloud_stats(ix,jy,2,n),jy=0,nymin1)
-!write(113,*) (icloud_stats(ix,jy,3,n),jy=0,nymin1)
-!write(114,*) (icloud_stats(ix,jy,4,n),jy=0,nymin1)
-!else
-!write(115,*) (cloudsh(ix,jy,n),jy=0,nymin1)    !integer
-!write(116,*) (lsprec(ix,jy,1,n),jy=0,nymin1)   !7.83691406E-02 
-!write(117,*) (convprec(ix,jy,1,n),jy=0,nymin1) !5.38330078E-02
-!endif
-!end do
-!
-!if (readclouds) then
-!CLOSE(111)
-!CLOSE(112)
-!CLOSE(113)
-!CLOSE(114)
-!else
-!CLOSE(115)
-!CLOSE(116)
-!CLOSE(117)
-!endif
-!
-!END ********* TEST *************** END
-! WRITE OUT SOME TEST VARIABLES
-!END ********* TEST *************** END
-
-
-! PS 2012
-!      lsp=lsprec(ix,jy,1,n)
-!      convp=convprec(ix,jy,1,n)
-!          prec=lsp+convp
-!          if (lsp.gt.convp) then !  prectype='lsp'
-!            lconvectprec = .false.
-!          else ! prectype='cp'
-!            lconvectprec = .true.
-!           endif
-!      else ! windfields does not contain cloud data 
-!          rhmin = 0.90 ! standard condition for presence of clouds
-!PS       note that original by Sabine Eckhart was 80%
-!PS       however, for T<-20 C we consider saturation over ice
-!PS       so I think 90% should be enough          
-!          icloudbot(ix,jy,n)=icmv
-!          icloudtop=icmv ! this is just a local variable
-!98        do kz=1,nz
-!            pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
-!            rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
-!ps            if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
-!            if (rh .gt. rhmin) then
-!              if (icloudbot(ix,jy,n) .eq. icmv) then
-!                icloudbot(ix,jy,n)=nint(height(kz))
-!              endif
-!              icloudtop=nint(height(kz)) ! use int to save memory
-!            endif
-    ! end do
-!PS try to get a cloud thicker than 50 m 
-!PS if there is at least .01 mm/h  - changed to 0.002 and put into
-!PS parameter precpmin        
-!          if ((icloudbot(ix,jy,n) .eq. icmv .or. &
-!              icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
-!              prec .gt. precmin) then
-!            rhmin = rhmin - 0.05
-!            if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
-!   end if
-
-!PS is based on looking at a limited set of comparison data
-!          if (lconvectprec .and. icloudtop .lt. 6000 .and. &
-!             prec .gt. precmin) then 
-!
-!            if (convp .lt. 0.1) then
-!              icloudbot(ix,jy,n) = 500
-!              icloudtop =         8000
-!            else
-!              icloudbot(ix,jy,n) = 0
-!              icloudtop =      10000
-!            endif
-!          endif
-!          if (icloudtop .ne. icmv) then
-!            icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
-!          else
-!            icloudthck(ix,jy,n) = icmv
-!          endif
-!PS  get rid of too thin clouds      
-!          if (icloudthck(ix,jy,n) .lt. 50) then
-!            icloudbot(ix,jy,n)=icmv
-!            icloudthck(ix,jy,n)=icmv
-!          endif
-!hg__________________________________
-!           rcw(ix,jy)=2E-7*prec**0.36
-!           rpc(ix,jy)=prec
-!hg end______________________________
-
-!      endif !hg read clouds
-
-
-
-
-!eso measure CPU time
-!  call mpif_mtime('verttransform',1)
-
-!eso print out the same measure as in Leo's routine
-    ! write(*,*) 'verttransform: ', &
-    !      sum(tt(:,:,:,n)*tt(:,:,:,n)), &
-    !      sum(uu(:,:,:,n)*uu(:,:,:,n)),sum(vv(:,:,:,n)*vv(:,:,:,n)),&
-    !      sum(qv(:,:,:,n)*qv(:,:,:,n)),sum(pv(:,:,:,n)*pv(:,:,:,n)),&
-    !      sum(rho(:,:,:,n)*rho(:,:,:,n)),sum(drhodz(:,:,:,n)*drhodz(:,:,:,n)),&
-    !      sum(ww(:,:,:,n)*ww(:,:,:,n)), &
-    !      sum(clouds(:,:,:,n)), sum(cloudsh(:,:,n)),sum(idx),sum(pinmconv)
-end subroutine verttransform
-
Index: src/verttransform_ecmwf.f90
===================================================================
--- src/verttransform_ecmwf.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
+++ src/verttransform_ecmwf.f90	(revision de4c5e9885c62eece05be19e1c0165ade923e201)
@@ -53,5 +53,6 @@
 ! Marian Harustak, 12.5.2017 
 !     - Renamed from verttransform to verttransform_ecmwf
-!
+!*****************************************************************************
+! Date: 2017-05-30 modification of a bug in ew. Don Morton (CTBTO project)   *
 !*****************************************************************************
 !                                                                            *
@@ -96,4 +97,7 @@
 
   logical :: init = .true.
+  logical :: init_w = .false.
+  logical :: init_r = .true.
+
 
   !ZHG SEP 2014 tests  
@@ -108,6 +112,6 @@
   ! CHARACTER(LEN=3)  :: aspec
   ! integer :: virr=0
-  real :: tot_cloud_h
-  real :: dbg_height(nzmax) 
+  !real :: tot_cloud_h
+  !real :: dbg_height(nzmax) 
 !ZHG
 
@@ -128,4 +132,17 @@
   if (init) then
 
+
+    if (init_r) then
+
+        open(333,file='heights.txt', &
+          form='formatted')
+        do kz=1,nuvz
+            read(333,*) height(kz)
+        end do
+        close(333)
+        write(*,*) 'height read'
+    else
+
+
 ! Search for a point with high surface pressure (i.e. not above significant topography)
 ! Then, use this point to construct a reference z profile, to be used at all times
@@ -166,4 +183,14 @@
     end do
 
+    if (init_w) then
+        open(333,file='heights.txt', &
+          form='formatted')
+        do kz=1,nuvz
+              write(333,*) height(kz)
+        end do
+        close(333)
+    endif
+
+    endif ! init
 
 ! Determine highest levels that can be within PBL
@@ -183,5 +210,5 @@
     init=.false.
 
-    dbg_height = height
+!    dbg_height = height
 
   endif
@@ -603,5 +630,5 @@
         convp=convprec(ix,jy,1,n)
         prec=lsp+convp
-        tot_cloud_h=0
+!        tot_cloud_h=0
 ! Find clouds in the vertical
         do kz=1, nz-1 !go from top to bottom
@@ -609,5 +636,5 @@
 ! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3 
             clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
-            tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz)) 
+!            tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz)) 
             
 !            icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]
