Index: src/FLEXPART.f90
===================================================================
--- src/FLEXPART.f90	(revision c8fc72499e835239b052197045c24d51657bb79f)
+++ src/FLEXPART.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -33,4 +33,12 @@
   !                                                                            *
   !*****************************************************************************
+  ! Changes:                                                                   *
+  !   Unified ECMWF and GFS builds                                             *
+  !   Marian Harustak, 12.5.2017                                               *
+  !     - Added detection of metdata format using gributils routines           *
+  !     - Distinguished calls to ecmwf/gfs gridcheck versions based on         *
+  !       detected metdata format                                              *
+  !     - Passed metdata format down to timemanager                            *
+  !*****************************************************************************
   !                                                                            *
   ! Variables:                                                                 *
@@ -46,4 +54,5 @@
   use netcdf_output_mod, only: writeheader_netcdf
   use random_mod, only: gasdev1
+  use class_gribfile
 
   implicit none
@@ -52,4 +61,7 @@
   integer :: idummy = -320
   character(len=256) :: inline_options  !pathfile, flexversion, arg2
+  integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
+  integer :: detectformat
+
 
 
@@ -172,4 +184,20 @@
   call readavailable
 
+  ! Detect metdata format
+  !**********************
+
+  metdata_format = detectformat()
+
+  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+    print *,'ECMWF metdata detected'
+  elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
+    print *,'NCEP metdata detected'
+  else
+    print *,'Unknown metdata format'
+    return
+  endif
+
+
+
   ! If nested wind fields are used, allocate arrays
   !************************************************
@@ -188,5 +216,9 @@
   endif
 
-  call gridcheck
+  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+    call gridcheck_ecmwf
+  else
+    call gridcheck_gfs
+  end if
 
   if (verbosity.gt.1) then   
@@ -411,5 +443,5 @@
   endif
 
-  call timemanager
+  call timemanager(metdata_format)
 
 ! NIK 16.02.2005 
Index: src/FLEXPART_MPI.f90
===================================================================
--- src/FLEXPART_MPI.f90	(revision 5184a7c7fc4985dfd336c4fd7fb8acacc26caef4)
+++ src/FLEXPART_MPI.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -33,4 +33,12 @@
   !                                                                            *
   !*****************************************************************************
+  ! Changes:                                                                   *
+  !   Unified ECMWF and GFS builds                                             *
+  !   Marian Harustak, 12.5.2017                                               *
+  !     - Added detection of metdata format using gributils routines           *
+  !     - Distinguished calls to ecmwf/gfs gridcheck versions based on         *
+  !       detected metdata format                                              *
+  !     - Passed metdata format down to timemanager                            *
+  !*****************************************************************************
   !                                                                            *
   ! Variables:                                                                 *
@@ -47,4 +55,5 @@
   use netcdf_output_mod, only: writeheader_netcdf
   use random_mod, only: gasdev1
+  use class_gribfile
 
   implicit none
@@ -53,4 +62,7 @@
   integer :: idummy = -320
   character(len=256) :: inline_options  !pathfile, flexversion, arg2
+  integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
+  integer :: detectformat
+
 
 
@@ -197,4 +209,20 @@
   call readavailable
 
+  ! Detect metdata format
+  !**********************
+
+  metdata_format = detectformat()
+
+  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+    print *,'ECMWF metdata detected'
+  elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
+    print *,'NCEP metdata detected'
+  else
+    print *,'Unknown metdata format'
+    return
+  endif
+
+
+
   ! If nested wind fields are used, allocate arrays
   !************************************************
@@ -213,5 +241,10 @@
   endif
 
-  call gridcheck
+  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+    call gridcheck_ecmwf
+  else
+    call gridcheck_gfs
+  end if
+
 
   if (verbosity.gt.1 .and. lroot) then   
@@ -456,5 +489,5 @@
 
 
-  call timemanager
+  call timemanager(metdata_format)
 
 
Index: src/calcmatrix.f90
===================================================================
--- src/calcmatrix.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/calcmatrix.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-subroutine calcmatrix(lconv,delt,cbmf)
+subroutine calcmatrix(lconv,delt,cbmf,metdata_format)
   !                        o    i    o
   !*****************************************************************************
@@ -30,11 +30,20 @@
   !  Petra Seibert, Bernd C. Krueger, 2000-2001                                *
   !                                                                            *
+  !*****************************************************************************
+  ! Changes:                                                                   *
   !  changed by C. Forster, November 2003 - February 2004                      *
   !  array fmassfrac(nconvlevmax,nconvlevmax) represents                       *
   !  the convective redistribution matrix for the particles                    *
   !                                                                            *
+  !   Unified ECMWF and GFS builds                                             *
+  !   Marian Harustak, 12.5.2017                                               *
+  !     - Merged calcmatrix and calcmatrix_gfs into one routine using if-then  *
+  !       for meteo-type dependent code                                        *
+  !*****************************************************************************
+  !                                                                            *
   !  lconv        indicates whether there is convection in this cell, or not   *
   !  delt         time step for convection [s]                                 *
   !  cbmf         cloud base mass flux                                         *
+  !  metdata_format format of metdata (ecmwf/gfs)                              *
   !                                                                            *
   !*****************************************************************************
@@ -43,8 +52,10 @@
   use com_mod
   use conv_mod
+  use class_gribfile
 
   implicit none
 
   real :: rlevmass,summe
+  integer :: metdata_format
 
   integer :: iflag, k, kk, kuvz
@@ -77,6 +88,10 @@
   do kuvz = 2,nuvz
     k = kuvz-1
+    if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv)
     phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv)
+    else
+      phconv(kuvz) =  0.5*(pconv(kuvz)+pconv(k))
+    endif
     dpr(k) = phconv(k) - phconv(kuvz)
     qsconv(k) = f_qvsat( pconv(k), tconv(k) )
Index: src/calcpar.f90
===================================================================
--- src/calcpar.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/calcpar.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-subroutine calcpar(n,uuh,vvh,pvh)
+subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
   !                   i  i   i   o
   !*****************************************************************************
@@ -37,11 +37,23 @@
   !     new variables in call to richardson                                    *
   !                                                                            *
-  !*****************************************************************************
-  !  Changes, Bernd C. Krueger, Feb. 2001:
-  !   Variables tth and qvh (on eta coordinates) in common block
+  !   CHANGE 17/11/2005 Caroline Forster NCEP GFS version                      *
+  !                                                                            *
+  !   Changes, Bernd C. Krueger, Feb. 2001:                                    *
+  !    Variables tth and qvh (on eta coordinates) in common block              *
+  !                                                                            *
+  !   Unified ECMWF and GFS builds                                             *
+  !   Marian Harustak, 12.5.2017                                               *
+  !     - Merged calcpar and calcpar_gfs into one routine using if-then        *
+  !       for meteo-type dependent code                                        *
+  !*****************************************************************************
+
   !*****************************************************************************
   !                                                                            *
   ! Variables:                                                                 *
   ! n                  temporal index for meteorological fields (1 to 3)       *
+  ! uuh                                                                        *
+  ! vvh                                                                        *
+  ! pvh                                                                        *
+  ! metdata_format     format of metdata (ecmwf/gfs)                           * 
   !                                                                            *
   ! Constants:                                                                 *
@@ -56,11 +68,13 @@
   use par_mod
   use com_mod
+  use class_gribfile
 
   implicit none
 
-  integer :: n,ix,jy,i,kz,lz,kzmin
+  integer :: metdata_format
+  integer :: n,ix,jy,i,kz,lz,kzmin,llev,loop_start
   real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
   real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
-  real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
+  real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax),hmixdummy
   real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
   real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
@@ -111,6 +125,23 @@
   !***********************************************
 
+      if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
+        ! NCEP version: find first level above ground
+        llev = 0
+        do i=1,nuvz
+          if (ps(ix,jy,1,n).lt.akz(i)) llev=i
+        end do
+        llev = llev+1
+        if (llev.gt.nuvz) llev = nuvz-1
+        ! NCEP version
+
+        ! calculate inverse Obukhov length scale with tth(llev)
       ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
-           tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm)
+            tth(ix,jy,llev,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm,akz(llev),metdata_format)
+      else
+        llev=0
+        ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
+            tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm,akz(llev),metdata_format)
+      end if
+
       if (ol.ne.0.) then
         oli(ix,jy,1,n)=1./ol
@@ -130,7 +161,14 @@
       end do
 
+      if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
+        ! NCEP version hmix has been read in in readwind.f, is therefore not calculated here
       call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
            ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), &
-           td2(ix,jy,1,n),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus)
+             td2(ix,jy,1,n),hmixdummy,wstar(ix,jy,1,n),hmixplus,metdata_format)
+      else
+        call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
+             ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), &
+             td2(ix,jy,1,n),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus,metdata_format)
+      end if
 
       if(lsubgrid.eq.1) then
@@ -173,6 +211,6 @@
   !******************************************************
 
-  ! 1) Calculate altitudes of ECMWF model levels
-  !*********************************************
+  ! 1) Calculate altitudes of model levels
+  !***************************************
 
       tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
@@ -180,5 +218,10 @@
       pold=ps(ix,jy,1,n)
       zold=0.
-      do kz=2,nuvz
+      if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+        loop_start=2
+      else
+        loop_start=llev
+      end if
+      do kz=loop_start,nuvz
         pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)  ! pressure on model layers
         tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
@@ -199,5 +242,11 @@
   !************************************************************************
 
-      do kz=1,nuvz
+      if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+        loop_start=1
+      else
+        loop_start=llev
+      end if
+
+      do kz=loop_start,nuvz
         if (zlev(kz).ge.altmin) then
           kzmin=kz
Index: src/calcpar_nests.f90
===================================================================
--- src/calcpar_nests.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
+++ src/calcpar_nests.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-subroutine calcpar_nests(n,uuhn,vvhn,pvhn)
+subroutine calcpar_nests(n,uuhn,vvhn,pvhn,metdata_format)
   !                         i  i    i    o
   !*****************************************************************************
@@ -39,10 +39,15 @@
   !                                                                            *
   !*****************************************************************************
-  !  Changes, Bernd C. Krueger, Feb. 2001:
-  !   Variables tth and qvh (on eta coordinates) in common block
+  !  Changes, Bernd C. Krueger, Feb. 2001:                                     *
+  !   Variables tth and qvh (on eta coordinates) in common block               *
+  !                                                                            *
+  !   Unified ECMWF and GFS builds                                             *
+  !   Marian Harustak, 12.5.2017                                               *
+  !     - Added passing of metdata_format as it was needed by called routines  *
   !*****************************************************************************
   !                                                                            *
   ! Variables:                                                                 *
   ! n                  temporal index for meteorological fields (1 to 3)       *
+  ! metdata_format     format of metdata (ecmwf/gfs)                           *
   !                                                                            *
   ! Constants:                                                                 *
@@ -60,6 +65,7 @@
   implicit none
 
+  integer :: metdata_format
   integer :: n,ix,jy,i,l,kz,lz,kzmin
-  real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
+  real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus,dummyakzllev
   real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
   real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
@@ -110,5 +116,5 @@
       ol=obukhov(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
            td2n(ix,jy,1,n,l),tthn(ix,jy,2,n,l),ustarn(ix,jy,1,n,l), &
-           sshfn(ix,jy,1,n,l),akm,bkm)
+           sshfn(ix,jy,1,n,l),akm,bkm,dummyakzllev,metdata_format)
       if (ol.ne.0.) then
         olin(ix,jy,1,n,l)=1./ol
@@ -131,5 +137,5 @@
            qvlev,ulev,vlev,nuvz,akz,bkz,sshfn(ix,jy,1,n,l), &
            tt2n(ix,jy,1,n,l),td2n(ix,jy,1,n,l),hmixn(ix,jy,1,n,l), &
-           wstarn(ix,jy,1,n,l),hmixplus)
+           wstarn(ix,jy,1,n,l),hmixplus,metdata_format)
 
       if(lsubgrid.eq.1) then
Index: src/convmix.f90
===================================================================
--- src/convmix.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/convmix.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-subroutine convmix(itime)
+subroutine convmix(itime,metdata_format)
   !                     i
   !**************************************************************
@@ -31,4 +31,14 @@
   !CHANGES by A. Stohl:
   !  various run-time optimizations - February 2005
+  ! CHANGES by C. Forster, November 2005, NCEP GFS version
+  !      in the ECMWF version convection is calculated on the
+  !      original eta-levels
+  !      in the GFS version convection is calculated on the
+  !      FLEXPART levels
+  !
+  !   Unified ECMWF and GFS builds                                             
+  !   Marian Harustak, 12.5.2017                                              
+  !     - Merged convmix and convmix_gfs into one routine using if-then           
+  !       for meteo-type dependent code                                        
   !**************************************************************
 
@@ -37,4 +47,5 @@
   use com_mod
   use conv_mod
+  use class_gribfile
 
   implicit none
@@ -44,4 +55,5 @@
   integer :: jy, kpart, ktop, ngrid,kz
   integer :: igrid(maxpart), ipoint(maxpart), igridn(maxpart,maxnests)
+  integer :: metdata_format
 
   ! itime [s]                 current time
@@ -104,4 +116,5 @@
 
     ngrid=0
+    if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     do j=numbnests,1,-1
       if ( x.gt.xln(j)+eps .and. x.lt.xrn(j)-eps .and. &
@@ -111,4 +124,13 @@
       endif
     end do
+    else
+      do j=numbnests,1,-1
+        if ( x.gt.xln(j) .and. x.lt.xrn(j) .and. &
+             y.gt.yln(j) .and. y.lt.yrn(j) ) then
+          ngrid=j
+          goto 23
+        endif
+      end do
+    endif
  23   continue
 
@@ -167,4 +189,5 @@
       td2conv=(td2(ix,jy,1,mind1)*dt2+td2(ix,jy,1,mind2)*dt1)*dtt
 !!$      do kz=1,nconvlev+1      !old
+      if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
         do kz=1,nuvz-1           !bugfix
         tconv(kz)=(tth(ix,jy,kz+1,mind1)*dt2+ &
@@ -173,7 +196,17 @@
              qvh(ix,jy,kz+1,mind2)*dt1)*dtt
       end do
+      else
+        do kz=1,nuvz-1           !bugfix
+          pconv(kz)=(pplev(ix,jy,kz,mind1)*dt2+ &
+              pplev(ix,jy,kz,mind2)*dt1)*dtt
+          tconv(kz)=(tt(ix,jy,kz,mind1)*dt2+ &
+              tt(ix,jy,kz,mind2)*dt1)*dtt
+          qconv(kz)=(qv(ix,jy,kz,mind1)*dt2+ &
+              qv(ix,jy,kz,mind2)*dt1)*dtt
+        end do
+      end if
 
   ! Calculate translocation matrix
-      call calcmatrix(lconv,delt,cbaseflux(ix,jy))
+      call calcmatrix(lconv,delt,cbaseflux(ix,jy),metdata_format)
       igrold = igr
       ktop = 0
@@ -252,5 +285,5 @@
   ! calculate translocation matrix
   !*******************************
-        call calcmatrix(lconv,delt,cbasefluxn(ix,jy,inest))
+        call calcmatrix(lconv,delt,cbasefluxn(ix,jy,inest),metdata_format)
         igrold = igr
         ktop = 0
Index: src/getfields.f90
===================================================================
--- src/getfields.f90	(revision d1a8707797c5c3fabdbd530210b84b85fb18b258)
+++ src/getfields.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-subroutine getfields(itime,nstop)
+subroutine getfields(itime,nstop,metdata_format)
 !                       i     o
 !*****************************************************************************
@@ -38,7 +38,11 @@
 !                                                                            *
 !*****************************************************************************
-!  Changes, Bernd C. Krueger, Feb. 2001:
-!   Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block.
-!   Function of nstop extended.
+!  Changes, Bernd C. Krueger, Feb. 2001:                                     *
+!   Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block.        *
+!   Function of nstop extended.                                              *
+!                                                                            *
+!   Unified ECMWF and GFS builds                                             *
+!   Marian Harustak, 12.5.2017                                               *
+!     - Added passing of metdata_format as it was needed by called routines  *
 !*****************************************************************************
 !                                                                            *
@@ -58,4 +62,5 @@
 ! tt(0:nxmax,0:nymax,nuvzmax,2)   temperature [K]                            *
 ! ps(0:nxmax,0:nymax,2)           surface pressure [Pa]                      *
+! metdata_format     format of metdata (ecmwf/gfs)                           *
 !                                                                            *
 ! Constants:                                                                 *
@@ -67,8 +72,10 @@
   use par_mod
   use com_mod
+  use class_gribfile
 
   implicit none
 
   integer :: indj,itime,nstop,memaux
+  integer :: metdata_format
 
   real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
@@ -125,9 +132,17 @@
     do indj=indmin,numbwf-1
       if (ldirect*wftime(indj+1).gt.ldirect*itime) then
-        call readwind(indj+1,memind(2),uuh,vvh,wwh)
+        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+          call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
+        else
+          call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh)
+        end if
         call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
-        call calcpar(memind(2),uuh,vvh,pvh)
-        call calcpar_nests(memind(2),uuhn,vvhn,pvhn)
-        call verttransform(memind(2),uuh,vvh,wwh,pvh)
+        call calcpar(memind(2),uuh,vvh,pvh,metdata_format)
+        call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format)
+        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+          call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh)
+        else
+          call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh)
+        end if
         call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
         memtime(2)=wftime(indj+1)
@@ -152,17 +167,33 @@
            (ldirect*wftime(indj+1).gt.ldirect*itime)) then
         memind(1)=1
-        call readwind(indj,memind(1),uuh,vvh,wwh)
+        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+          call readwind_ecmwf(indj,memind(1),uuh,vvh,wwh)
+        else
+          call readwind_gfs(indj,memind(1),uuh,vvh,wwh)
+        end if
         call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn)
-        call calcpar(memind(1),uuh,vvh,pvh)
-        call calcpar_nests(memind(1),uuhn,vvhn,pvhn)
-        call verttransform(memind(1),uuh,vvh,wwh,pvh)
+        call calcpar(memind(1),uuh,vvh,pvh,metdata_format)
+        call calcpar_nests(memind(1),uuhn,vvhn,pvhn,metdata_format)
+        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+          call verttransform_ecmwf(memind(1),uuh,vvh,wwh,pvh)
+        else
+          call verttransform_gfs(memind(1),uuh,vvh,wwh,pvh)
+        end if
         call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
         memtime(1)=wftime(indj)
         memind(2)=2
-        call readwind(indj+1,memind(2),uuh,vvh,wwh)
+        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+          call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
+        else
+          call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh)
+        end if
         call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
-        call calcpar(memind(2),uuh,vvh,pvh)
-        call calcpar_nests(memind(2),uuhn,vvhn,pvhn)
-        call verttransform(memind(2),uuh,vvh,wwh,pvh)
+        call calcpar(memind(2),uuh,vvh,pvh,metdata_format)
+        call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format)
+        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+          call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh)
+        else
+          call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh)
+        end if
         call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
         memtime(2)=wftime(indj+1)
Index: src/getfields_mpi.f90
===================================================================
--- src/getfields_mpi.f90	(revision 78e62dc004da36ae5da7b3ab5f03e1bd44af4586)
+++ src/getfields_mpi.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-subroutine getfields(itime,nstop,memstat)
+subroutine getfields(itime,nstop,memstat,metdata_format)
 !                       i     o       o
 !*****************************************************************************
@@ -58,4 +58,8 @@
 !    memstat=0:      no new fields to be read
 !          
+!   Unified ECMWF and GFS builds                                             
+!   Marian Harustak, 12.5.2017                                               
+!     - Added passing of metdata_format as it was needed by called routines  
+!          
 !*****************************************************************************
 !                                                                            *
@@ -77,4 +81,5 @@
 ! tt(0:nxmax,0:nymax,nuvzmax,numwfmem)   temperature [K]                     *
 ! ps(0:nxmax,0:nymax,numwfmem)           surface pressure [Pa]               *
+! metdata_format     format of metdata (ecmwf/gfs)                           *
 !                                                                            *
 ! Constants:                                                                 *
@@ -87,7 +92,9 @@
   use com_mod
   use mpi_mod, only: lmpreader,lmp_use_reader,lmp_sync
+  use class_gribfile
 
   implicit none
 
+  integer :: metdata_format
   integer :: indj,itime,nstop,memaux,mindread
 ! mindread: index of where to read 3rd field
@@ -204,9 +211,17 @@
       if (ldirect*wftime(indj+1).gt.ldirect*itime) then
         if (lmpreader.or..not. lmp_use_reader) then
-          call readwind(indj+1,mindread,uuh,vvh,wwh)
+          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+            call readwind_ecmwf(indj+1,mindread,uuh,vvh,wwh)
+          else
+            call readwind_gfs(indj+1,mindread,uuh,vvh,wwh)
+          end if
           call readwind_nests(indj+1,mindread,uuhn,vvhn,wwhn)
-          call calcpar(mindread,uuh,vvh,pvh)
-          call calcpar_nests(mindread,uuhn,vvhn,pvhn)
-          call verttransform(mindread,uuh,vvh,wwh,pvh)
+          call calcpar(mindread,uuh,vvh,pvh,metdata_format)
+          call calcpar_nests(mindread,uuhn,vvhn,pvhn,metdata_format)
+          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+            call verttransform_ecmwf(mindread,uuh,vvh,wwh,pvh)
+          else
+            call verttransform_gfs(mindread,uuh,vvh,wwh,pvh)
+          end if
           call verttransform_nests(mindread,uuhn,vvhn,wwhn,pvhn)
         end if
@@ -230,9 +245,17 @@
         memind(1)=1
         if (lmpreader.or..not.lmp_use_reader) then
-          call readwind(indj,memind(1),uuh,vvh,wwh)
+          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+            call readwind_ecmwf(indj,memind(1),uuh,vvh,wwh)
+          else
+            call readwind_gfs(indj,memind(1),uuh,vvh,wwh)
+          end if
           call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn)
-          call calcpar(memind(1),uuh,vvh,pvh)
-          call calcpar_nests(memind(1),uuhn,vvhn,pvhn)
-          call verttransform(memind(1),uuh,vvh,wwh,pvh)
+          call calcpar(memind(1),uuh,vvh,pvh,metdata_format)
+          call calcpar_nests(memind(1),uuhn,vvhn,pvhn,metdata_format)
+          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+            call verttransform_ecmwf(memind(1),uuh,vvh,wwh,pvh)
+          else
+            call verttransform_gfs(memind(1),uuh,vvh,wwh,pvh)
+          end if
           call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
         end if
@@ -240,9 +263,17 @@
         memind(2)=2
         if (lmpreader.or..not.lmp_use_reader) then
-          call readwind(indj+1,memind(2),uuh,vvh,wwh)
+          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+            call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
+          else
+            call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh)
+          end if
           call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
-          call calcpar(memind(2),uuh,vvh,pvh)
-          call calcpar_nests(memind(2),uuhn,vvhn,pvhn)
-          call verttransform(memind(2),uuh,vvh,wwh,pvh)
+          call calcpar(memind(2),uuh,vvh,pvh,metdata_format)
+          call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format)
+          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+            call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh)
+          else
+            call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh)
+          end if
           call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
         end if
Index: src/gridcheck_gfs.f90
===================================================================
--- src/gridcheck_gfs.f90	(revision 4c64400c622119505b289ef83be130b0b81b9590)
+++ src/gridcheck_gfs.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-subroutine gridcheck
+subroutine gridcheck_gfs
 
   !**********************************************************************
@@ -38,4 +38,8 @@
   !             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
   !                                 ECMWF grib_api                      *
+  !                                                                     *
+  !   Unified ECMWF and GFS builds                                      *
+  !   Marian Harustak, 12.5.2017                                        *
+  !     - Renamed routine from gridcheck to gridcheck_gfs               *
   !                                                                     *
   !**********************************************************************
@@ -537,3 +541,3 @@
   endif
 
-end subroutine gridcheck
+end subroutine gridcheck_gfs
Index: src/makefile
===================================================================
--- src/makefile	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/makefile	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -143,28 +143,28 @@
 ew.o			readreleases.o  \
 readdepo.o              get_vdep_prob.o   \
-get_wetscav.o	\
+get_wetscav.o		readwind_gfs.o \
 psim.o			outgrid_init.o 	\
-outgrid_init_nest.o   	\
+outgrid_init_nest.o   	calcmatrix.o \
 photo_O1D.o 		readlanduse.o \
 interpol_wind.o         readoutgrid.o \
 interpol_all.o          readpaths.o \
-getrb.o                 \
-getrc.o                 \
+getrb.o                 obukhov.o \
+getrc.o                 convmix.o \
 getvdep.o               readspecies.o \
-interpol_misslev.o      \
-scalev.o \
-pbl_profile.o           readOHfield.o\
-juldate.o               \
+interpol_misslev.o      richardson.o \
+scalev.o                verttransform_ecmwf.o \
+pbl_profile.o           readOHfield.o \
+juldate.o               verttransform_gfs.o \
 interpol_vdep.o         interpol_rain.o \
 hanna.o                 wetdepokernel.o \
-                  wetdepo.o \
+calcpar.o               wetdepo.o \
 hanna_short.o           windalign.o \
-hanna1.o                	\
-                        gridcheck_nests.o \
+hanna1.o                gridcheck_ecmwf.o \
+gridcheck_gfs.o         gridcheck_nests.o \
 readwind_nests.o        calcpar_nests.o \
 verttransform_nests.o   interpol_all_nests.o \
 interpol_wind_nests.o   interpol_misslev_nests.o \
 interpol_vdep_nests.o   interpol_rain_nests.o \
-readageclasses.o         \
+readageclasses.o        detectformat.o  \
 calcfluxes.o            fluxoutput.o \
 qvsat.o                 skplin.o \
@@ -319,6 +319,5 @@
 mpi_mod.o: com_mod.o par_mod.o unc_mod.o
 netcdf_output_mod.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
-obukhov.o: par_mod.o
-obukhov_gfs.o: par_mod.o
+obukhov.o: par_mod.o class_gribfile_mod.o
 ohreaction.o: com_mod.o oh_mod.o par_mod.o
 openouttraj.o: com_mod.o par_mod.o point_mod.o
@@ -327,5 +326,4 @@
 outgrid_init.o: com_mod.o flux_mod.o oh_mod.o outg_mod.o par_mod.o unc_mod.o
 outgrid_init_nest.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
-par_mod.o :
 part0.o: par_mod.o
 partdep.o: par_mod.o
Index: src/mpi_mod.f90
===================================================================
--- src/mpi_mod.f90	(revision b5127f9034f2a15f9a3ce226fd26edfe845311de)
+++ src/mpi_mod.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -59,4 +59,11 @@
 !                                                                            *
 !*****************************************************************************
+!                                                                            *
+! Modification by DJM, 2017-05-09 - added #ifdef USE_MPIINPLACE cpp          *
+! directive to mpif_tm_reduce_grid() to insure that MPI_IN_PLACE is          *
+! used in the MPI_Reduce() only if specifically compiled with that           *
+! directive.                                                                 *
+!                                                                            *
+!*****************************************************************************
 
   use mpi 
@@ -2421,4 +2428,25 @@
 
 ! 2) Using in-place reduction
+
+!!!--------------------------------------------------------------------
+!!! DJM - 2017-05-09 change - MPI_IN_PLACE option for MPI_Reduce() causes
+!!! severe numerical problems in some cases.  I'm guessing there are memory
+!!! overrun issues in this complex code, but have so far failed to identify
+!!! a specific problem.  And/or, when one searches the Internet for this 
+!!! issue, there is "some" hint that the implementation may be buggy.  
+!!! 
+!!! At this point, with the following CPP USE_MPIINPLACE directive, the
+!!! default behaviour will be to not use the MPI_IN_PLACE option.
+!!! Users will have to compile with -DUSE_MPIINPLACE if they want that option.
+!!! Introduction of the CPP directives also requires that the code be compiled
+!!! with the "-x f95-cpp-input" option.
+!!!
+!!! Modification of this section requires the addition of array gridunc0, which
+!!! requires an allocation in outgrid_init.f90 and initial declaration in
+!!! unc_mod.f90.
+!!!--------------------------------------------------------------------
+
+#ifdef USE_MPIINPLACE
+
     if (lroot) then
       call MPI_Reduce(MPI_IN_PLACE, gridunc, grid_size3d, mp_sp, MPI_SUM, id_root, &
@@ -2429,4 +2457,12 @@
            & mp_comm_used, mp_ierr)
     end if
+
+#else
+
+      call MPI_Reduce(gridunc, gridunc0, grid_size3d, mp_sp, MPI_SUM, id_root, &
+           & mp_comm_used, mp_ierr)
+      if (lroot) gridunc = gridunc0
+
+#endif
 
     if ((WETDEP).and.(ldirect.gt.0)) then
Index: src/obukhov.f90
===================================================================
--- src/obukhov.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/obukhov.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-real function obukhov(ps,tsurf,tdsurf,tlev,ustar,hf,akm,bkm)
+real function obukhov(ps,tsurf,tdsurf,tlev,ustar,hf,akm,bkm,plev,metdata_format)
 
   !********************************************************************
@@ -27,11 +27,17 @@
   !                       Date:   1994-06-27                          *
   !                                                                   *
-  !  Update: A. Stohl, 2000-09-25, avoid division by zero by          *
-  !  setting ustar to minimum value                                   *
+  !     This program calculates Obukhov scale height from surface     *
+  !     meteorological data and sensible heat flux.                   *
   !                                                                   *
   !********************************************************************
   !                                                                   *
-  !     This program calculates Obukhov scale height from surface     *
-  !     meteorological data and sensible heat flux.                   *
+  !  Update: A. Stohl, 2000-09-25, avoid division by zero by          *
+  !  setting ustar to minimum value                                   *
+  !  CHANGE: 17/11/2005 Caroline Forster NCEP GFS version             *
+  !                                                                   *
+  !   Unified ECMWF and GFS builds                                    *
+  !   Marian Harustak, 12.5.2017                                      *
+  !     - Merged obukhov and obukhov_gfs into one routine using       *
+  !       if-then for meteo-type dependent code                       *
   !                                                                   *
   !********************************************************************
@@ -47,11 +53,15 @@
   !     akm     ECMWF vertical discretization parameter               *
   !     bkm     ECMWF vertical discretization parameter               *
+  !     plev                                                          *
+  !     metdata_format format of metdata (ecmwf/gfs)                  *
   !                                                                   *
   !********************************************************************
 
   use par_mod
+  use class_gribfile
 
   implicit none
 
+  integer :: metdata_format
   real :: akm(nwzmax),bkm(nwzmax)
   real :: ps,tsurf,tdsurf,tlev,ustar,hf,e,ew,tv,rhoa,plev
@@ -62,7 +72,9 @@
   tv=tsurf*(1.+0.378*e/ps)               ! virtual temperature
   rhoa=ps/(r_air*tv)                      ! air density
+  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
   ak1=(akm(1)+akm(2))/2.
   bk1=(bkm(1)+bkm(2))/2.
   plev=ak1+bk1*ps                        ! Pressure level 1
+  end if
   theta=tlev*(100000./plev)**(r_air/cpa) ! potential temperature
   if (ustar.le.0.) ustar=1.e-8
Index: src/outgrid_init.f90
===================================================================
--- src/outgrid_init.f90	(revision 6b22af900a57ac3fc1a984db76f4ef2cb0436af0)
+++ src/outgrid_init.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -19,4 +19,11 @@
 ! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
 !**********************************************************************
+!                                                                     *
+! DJM - 2017-05-09 - added #ifdef USE_MPIINPLACE cpp directive to     *
+! enable allocation of a gridunc0 array if required by MPI code in    *
+! mpi_mod.f90                                                         *
+!                                                                     *
+!**********************************************************************
+
 
 subroutine outgrid_init
@@ -215,7 +222,13 @@
 ! Extra field for totals at MPI root process
   if (lroot.and.mpi_mode.gt.0) then
-    ! allocate(gridunc0(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, &
-    !      maxpointspec_act,nclassunc,maxageclass),stat=stat)
-    ! if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc0'
+
+#ifdef USE_MPIINPLACE
+#else
+    ! If MPI_IN_PLACE option is not used in mpi_mod.f90::mpif_tm_reduce_grid(),
+    ! then an aux array is needed for parallel grid reduction
+    allocate(gridunc0(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, &
+         maxpointspec_act,nclassunc,maxageclass),stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc0'
+#endif
     if (ldirect.gt.0) then
       allocate(wetgridunc0(0:numxgrid-1,0:numygrid-1,maxspec, &
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/par_mod.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -142,5 +142,7 @@
   
 !  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF new 
-  integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 
+!  integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 
+
+  integer,parameter :: nxmax=181,nymax=91,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 
 
 !  INTEGER,PARAMETER :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !NCEP data
@@ -149,6 +151,6 @@
 !  integer,parameter :: nxshift=0     ! for GFS
 !  integer,parameter :: nxmax=15,nymax=15,nuvzmax=140,nwzmax=140,nzmax=140
-  integer,parameter :: nxshift=359  ! for ECMWF
-!  integer,parameter :: nxshift=0     ! for GFS
+!  integer,parameter :: nxshift=359  ! for ECMWF
+  integer,parameter :: nxshift=0     ! for GFS
 
   !*********************************************
Index: src/readwind_gfs.f90
===================================================================
--- src/readwind_gfs.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/readwind_gfs.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-subroutine readwind(indj,n,uuh,vvh,wwh)
+subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
 
   !***********************************************************************
@@ -38,4 +38,8 @@
   !*             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
   !*                                 ECMWF grib_api                      *
+  !                                                                      *
+  !   Unified ECMWF and GFS builds                                       *
+  !   Marian Harustak, 12.5.2017                                         *
+  !     - Renamed routine from readwind to readwind_gfs                  *
   !*                                                                     *
   !***********************************************************************
@@ -716,3 +720,3 @@
   stop 'Execution terminated'
 
-end subroutine readwind
+end subroutine readwind_gfs
Index: src/richardson.f90
===================================================================
--- src/richardson.f90	(revision c2162ce5e4e231a653816e16d96eacf839a90235)
+++ src/richardson.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -21,5 +21,5 @@
 
 subroutine richardson(psurf,ust,ttlev,qvlev,ulev,vlev,nuvz, &
-       akz,bkz,hf,tt2,td2,h,wst,hmixplus)
+       akz,bkz,hf,tt2,td2,h,wst,hmixplus,metdata_format)
   !                        i    i    i     i    i    i    i
   ! i   i  i   i   i  o  o     o
@@ -41,4 +41,6 @@
   !     Meteor. 81, 245-269.                                                  *
   !                                                                           *
+  !****************************************************************************
+  !                                                                           *
   !     Update: 1999-02-01 by G. Wotawa                                       *
   !                                                                           *
@@ -46,4 +48,11 @@
   !     instead of first model level.                                         *
   !     New input variables tt2, td2 introduced.                              *
+  !                                                                           *
+  !     CHANGE: 17/11/2005 Caroline Forster NCEP GFS version                  * 
+  !                                                                           *
+  !     Unified ECMWF and GFS builds                                          *
+  !     Marian Harustak, 12.5.2017                                            *
+  !       - Merged richardson and richardson_gfs into one routine using       *
+  !         if-then for meteo-type dependent code                             *
   !                                                                           *
   !****************************************************************************
@@ -55,4 +64,5 @@
   ! tv                         virtual temperature                            *
   ! wst                        convective velocity scale                      *
+  ! metdata_format             format of metdata (ecmwf/gfs)                  *
   !                                                                           *
   ! Constants:                                                                *
@@ -62,8 +72,10 @@
 
   use par_mod
+  use class_gribfile
 
   implicit none
 
-  integer :: i,k,nuvz,iter
+  integer :: metdata_format
+  integer :: i,k,nuvz,iter,llev,loop_start
   real :: tv,tvold,zref,z,zold,pint,pold,theta,thetaref,ri
   real :: akz(nuvz),bkz(nuvz),ulev(nuvz),vlev(nuvz),hf,wst,tt2,td2,ew
@@ -77,4 +89,20 @@
   iter=0
 
+  if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
+    ! NCEP version: find first model level above ground
+    !**************************************************
+
+     llev = 0
+     do i=1,nuvz
+       if (psurf.lt.akz(i)) llev=i
+     end do
+     llev = llev+1
+    ! sec llev should not be 1!
+     if (llev.eq.1) llev = 2
+     if (llev.gt.nuvz) llev = nuvz-1
+    ! NCEP version
+  end if
+
+
   ! Compute virtual temperature and virtual potential temperature at
   ! reference level (2 m)
@@ -95,6 +123,10 @@
   ! Integrate z up to one level above zt
   !*************************************
-
-  do k=2,nuvz
+  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+    loop_start=2
+  else
+    loop_start=llev
+  end if
+  do k=loop_start,nuvz
     pint=akz(k)+bkz(k)*psurf  ! pressure on model layers
     tv=ttlev(k)*(1.+0.608*qvlev(k))
Index: src/timemanager.f90
===================================================================
--- src/timemanager.f90	(revision a76d954a64a27c97f4b2204adec2e05b311283f5)
+++ src/timemanager.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-subroutine timemanager
+subroutine timemanager(metdata_format)
 
   !*****************************************************************************
@@ -50,4 +50,7 @@
   !   For compatibility with MPI version,                                      *
   !   variables uap,ucp,uzp,us,vs,ws,cbt now in module com_mod                 *
+  !  Unified ECMWF and GFS builds                                              *
+  !   Marian Harustak, 12.5.2017                                               *
+  !   - Added passing of metdata_format as it was needed by called routines    *
   !*****************************************************************************
   !                                                                            *
@@ -83,4 +86,5 @@
   ! xtra1(maxpart), ytra1(maxpart), ztra1(maxpart) =                           *
   !                    spatial positions of trajectories                       *
+  ! metdata_format     format of metdata (ecmwf/gfs)                           *
   !                                                                            *
   ! Constants:                                                                 *
@@ -102,4 +106,5 @@
   implicit none
 
+  integer :: metdata_format
   integer :: j,ks,kp,l,n,itime=0,nstop,nstop1
 ! integer :: ksp
@@ -194,5 +199,5 @@
            write (*,*) 'timemanager> call convmix -- backward'
         endif         
-      call convmix(itime)
+      call convmix(itime,metdata_format)
         if (verbosity.gt.1) then
           !CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
@@ -207,5 +212,5 @@
            write (*,*) 'timemanager> call getfields'
     endif 
-    call getfields(itime,nstop1)
+    call getfields(itime,nstop1,metdata_format)
         if (verbosity.gt.1) then
           CALL SYSTEM_CLOCK(count_clock)
@@ -266,5 +271,5 @@
        write (*,*) 'timemanager> call convmix -- forward'
      endif    
-     call convmix(itime)
+     call convmix(itime,metdata_format)
    endif
 
Index: src/timemanager_mpi.f90
===================================================================
--- src/timemanager_mpi.f90	(revision b5127f9034f2a15f9a3ce226fd26edfe845311de)
+++ src/timemanager_mpi.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-subroutine timemanager
+subroutine timemanager(metdata_format)
 
 !*****************************************************************************
@@ -50,4 +50,7 @@
 !   MPI version                                                              *
 !   Variables uap,ucp,uzp,us,vs,ws,cbt now in module com_mod                 *
+!  Unified ECMWF and GFS builds                                              *
+!   Marian Harustak, 12.5.2017                                               *
+!   - Added passing of metdata_format as it was needed by called routines    *
 !*****************************************************************************
 !                                                                            *
@@ -205,5 +208,5 @@
       endif
 ! readwind process skips this step
-      if (.not.(lmpreader.and.lmp_use_reader)) call convmix(itime)
+      if (.not.(lmpreader.and.lmp_use_reader)) call convmix(itime,metdata_format)
     endif
 
@@ -218,5 +221,5 @@
     if (mp_measure_time) call mpif_mtime('getfields',0)
 
-    call getfields(itime,nstop1,memstat)
+    call getfields(itime,nstop1,memstat,metdata_format)
 
     if (mp_measure_time) call mpif_mtime('getfields',1)
Index: src/unc_mod.f90
===================================================================
--- src/unc_mod.f90	(revision 4c64400c622119505b289ef83be130b0b81b9590)
+++ src/unc_mod.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -19,4 +19,10 @@
 ! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
 !**********************************************************************
+!                                                                     *
+! DJM - 2017-05-09 - added #ifdef USE_MPIINPLACE cpp directive to     *
+! enable declaration of a gridunc0 array if required by MPI code in   *
+! mpi_mod.f90                                                         *
+!                                                                     *
+!**********************************************************************
 
 module unc_mod
@@ -27,4 +33,10 @@
 
   real,allocatable, dimension (:,:,:,:,:,:,:) :: gridunc
+#ifdef USE_MPIINPLACE
+#else
+  ! If MPI_IN_PLACE option is not used in mpi_mod.f90::mpif_tm_reduce_grid(),
+  ! then an aux array is needed for parallel grid reduction
+  real,allocatable, dimension (:,:,:,:,:,:,:) :: gridunc0
+#endif
   real,allocatable, dimension (:,:,:,:,:,:,:) :: griduncn
   real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: drygridunc
Index: src/verttransform_ecmwf.f90
===================================================================
--- src/verttransform_ecmwf.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/verttransform_ecmwf.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -144,5 +144,5 @@
 
 
-    tvold(ixm,jym)=tt2(ixm,jym,1,n)*(1.+0.378*ew*(td2(ixm,jym,1,n))/ &
+    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)
@@ -194,5 +194,5 @@
   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))/ &
+      tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
            ps(ix,jy,1,n))
     enddo
Index: src/verttransform_gfs.f90
===================================================================
--- src/verttransform_gfs.f90	(revision 4fbe7a5609fe43bdc3940768d8f4d949237cc44a)
+++ src/verttransform_gfs.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -20,5 +20,5 @@
 !**********************************************************************
 
-subroutine verttransform(n,uuh,vvh,wwh,pvh)
+subroutine verttransform_gfs(n,uuh,vvh,wwh,pvh)
   !                      i  i   i   i   i
   !*****************************************************************************
@@ -47,4 +47,9 @@
   !  Changes, Bernd C. Krueger, Feb. 2001:
   !   Variables tth and qvh (on eta coordinates) from common block
+  !
+  !   Unified ECMWF and GFS builds                                      
+  !   Marian Harustak, 12.5.2017                                        
+  !     - Renamed routine from verttransform to verttransform_gfs
+  !
   !*****************************************************************************
   !                                                                            *
@@ -531,3 +536,3 @@
 
 
-end subroutine verttransform
+end subroutine verttransform_gfs
Index: src/writeheader.f90
===================================================================
--- src/writeheader.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
+++ src/writeheader.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -35,4 +35,12 @@
   !*****************************************************************************
   !                                                                            *
+  !  Modified to remove TRIM around the output of flexversion so that          *
+  !  it will be a constant length (defined in com_mod.f90) in output header    *
+  !                                                                            *
+  !     Don Morton, Boreal Scientific Computing                                *
+  !     07 May 2017                                                            *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
   ! Variables:                                                                 *
   !                                                                            *
@@ -67,7 +75,7 @@
 
   if (ldirect.eq.1) then
-    write(unitheader) ibdate,ibtime, trim(flexversion)
+    write(unitheader) ibdate,ibtime, flexversion
   else
-    write(unitheader) iedate,ietime, trim(flexversion)
+    write(unitheader) iedate,ietime, flexversion
   endif
 
Index: src/writeheader_nest.f90
===================================================================
--- src/writeheader_nest.f90	(revision 4fbe7a5609fe43bdc3940768d8f4d949237cc44a)
+++ src/writeheader_nest.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -35,4 +35,12 @@
   !*****************************************************************************
   !                                                                            *
+  !  Modified to remove TRIM around the output of flexversion so that          *
+  !  it will be a constant length (defined in com_mod.f90) in output header    *
+  !                                                                            *
+  !     Don Morton, Boreal Scientific Computing                                *
+  !     07 May 2017                                                            *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
   ! Variables:                                                                 *
   !                                                                            *
@@ -67,7 +75,7 @@
 
   if (ldirect.eq.1) then
-    write(unitheader) ibdate,ibtime, trim(flexversion)
+    write(unitheader) ibdate,ibtime, flexversion
   else
-    write(unitheader) iedate,ietime, trim(flexversion)
+    write(unitheader) iedate,ietime, flexversion
   endif
 
Index: src/writeheader_nest_surf.f90
===================================================================
--- src/writeheader_nest_surf.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
+++ src/writeheader_nest_surf.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -35,4 +35,12 @@
   !*****************************************************************************
   !                                                                            *
+  !  Modified to remove TRIM around the output of flexversion so that          *
+  !  it will be a constant length (defined in com_mod.f90) in output header    *
+  !                                                                            *
+  !     Don Morton, Boreal Scientific Computing                                *
+  !     07 May 2017                                                            *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
   ! Variables:                                                                 *
   !                                                                            *
@@ -67,7 +75,7 @@
 
   if (ldirect.eq.1) then
-    write(unitheader) ibdate,ibtime,trim(flexversion)
+    write(unitheader) ibdate,ibtime,flexversion
   else
-    write(unitheader) iedate,ietime,trim(flexversion)
+    write(unitheader) iedate,ietime,flexversion
   endif
 
Index: src/writeheader_surf.f90
===================================================================
--- src/writeheader_surf.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
+++ src/writeheader_surf.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
@@ -35,4 +35,12 @@
   !*****************************************************************************
   !                                                                            *
+  !  Modified to remove TRIM around the output of flexversion so that          *
+  !  it will be a constant length (defined in com_mod.f90) in output header    *
+  !                                                                            *
+  !     Don Morton, Boreal Scientific Computing                                *
+  !     07 May 2017                                                            *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
   ! Variables:                                                                 *
   !                                                                            *
@@ -67,7 +75,7 @@
 
   if (ldirect.eq.1) then
-    write(unitheader) ibdate,ibtime, trim(flexversion)
+    write(unitheader) ibdate,ibtime, flexversion
   else
-    write(unitheader) iedate,ietime, trim(flexversion)
+    write(unitheader) iedate,ietime, flexversion
   endif
 
