Index: src/calcmatrix.f90
===================================================================
--- src/calcmatrix.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/calcmatrix.f90	(revision 9ca6e380ca9ee104eba37662c43b9e13b165195c)
@@ -67,22 +67,26 @@
 
   phconv(1) = psconv
-  ! Emanuel subroutine needs pressure in hPa, therefore convert all pressures
-  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) )
+! Emanuel subroutine needs pressure in hPa, therefore convert all pressures
+  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
+    do kuvz = 2,nuvz
+      k = kuvz-1
+      pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv)
+      phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv)
+      dpr(k) = phconv(k) - phconv(kuvz)
+      qsconv(k) = f_qvsat( pconv(k), tconv(k) )
+    end do
+  else
+! JMA / SH Bugfix phconv was set in loop with access to undefined pconv(nuvz)
+    phconv(2:nuvz-1) = 0.5*(pconv(2:nuvz-1)+pconv(1:nuvz-2))
+    phconv(nuvz) = pconv(nuvz-1)
+    dpr(1:nuvz-1) = phconv(1:nuvz-1)-phconv(2:nuvz)
 
-  ! initialize mass fractions
-    do kk=1,nconvlev
-      fmassfrac(k,kk)=0.
+    do k = 1,nuvz-1
+      qsconv(k) = f_qvsat( pconv(k), tconv(k) )
     end do
-  end do
-
+  end if
+    
+! initialize mass fractions
+  fmassfrac(1:nuvz-1,1:nconvlev)=0.0
 
   !note that Emanuel says it is important
@@ -93,47 +97,47 @@
   !******************
 
-    cbmfold = cbmf
-  ! Convert pressures to hPa, as required by Emanuel scheme
-  !********************************************************
+  cbmfold = cbmf
+! Convert pressures to hPa, as required by Emanuel scheme
+!********************************************************
 !!$    do k=1,nconvlev     !old
-    do k=1,nconvlev+1      !bugfix
-      pconv_hpa(k)=pconv(k)/100.
-      phconv_hpa(k)=phconv(k)/100.
+  do k=1,nconvlev+1      !bugfix
+    pconv_hpa(k)=pconv(k)/100.
+    phconv_hpa(k)=phconv(k)/100.
+  end do
+  phconv_hpa(nconvlev+1)=phconv(nconvlev+1)/100.
+  call convect(nconvlevmax, nconvlev, delt, iflag, &
+       precip, wd, tprime, qprime, cbmf)
+
+! do not update fmassfrac and cloudbase massflux
+! if no convection takes place or
+! if a CFL criterion is violated in convect43c.f
+  if (iflag .ne. 1 .and. iflag .ne. 4) then
+    cbmf=cbmfold
+    goto 200
+  endif
+
+! do not update fmassfrac and cloudbase massflux
+! if the old and the new cloud base mass
+! fluxes are zero
+  if (cbmf.le.0..and.cbmfold.le.0.) then
+    cbmf=cbmfold
+    goto 200
+  endif
+
+! Update fmassfrac
+! account for mass displaced from level k to level k
+
+  lconv = .true.
+  do k=1,nconvtop
+    rlevmass = dpr(k)/ga
+    summe = 0.
+    do kk=1,nconvtop
+      fmassfrac(k,kk) = delt*fmass(k,kk)
+      summe = summe + fmassfrac(k,kk)
     end do
-    phconv_hpa(nconvlev+1)=phconv(nconvlev+1)/100.
-    call convect(nconvlevmax, nconvlev, delt, iflag, &
-         precip, wd, tprime, qprime, cbmf)
+    fmassfrac(k,k)=fmassfrac(k,k) + rlevmass - summe
+  end do
 
-  ! do not update fmassfrac and cloudbase massflux
-  ! if no convection takes place or
-  ! if a CFL criterion is violated in convect43c.f
-   if (iflag .ne. 1 .and. iflag .ne. 4) then
-     cbmf=cbmfold
-     goto 200
-   endif
-
-  ! do not update fmassfrac and cloudbase massflux
-  ! if the old and the new cloud base mass
-  ! fluxes are zero
-   if (cbmf.le.0..and.cbmfold.le.0.) then
-     cbmf=cbmfold
-     goto 200
-   endif
-
-  ! Update fmassfrac
-  ! account for mass displaced from level k to level k
-
-   lconv = .true.
-    do k=1,nconvtop
-      rlevmass = dpr(k)/ga
-      summe = 0.
-      do kk=1,nconvtop
-        fmassfrac(k,kk) = delt*fmass(k,kk)
-        summe = summe + fmassfrac(k,kk)
-      end do
-      fmassfrac(k,k)=fmassfrac(k,k) + rlevmass - summe
-    end do
-
-200   continue
+200 continue
 
 end subroutine calcmatrix
Index: src/ew.f90
===================================================================
--- src/ew.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/ew.f90	(revision 467460af7de11cc1002100a8e56de8b8108cd235)
@@ -14,5 +14,9 @@
 
   ew=0.
-  if(x.le.0.) stop 'sorry: t not in [k]'
+  if(x.le.0.) then
+    write(*,*) 'in ew.f90: x=', x
+    stop 'sorry: t not in [k]'
+  end if
+
   y=373.16/x
   a=-7.90298*(y-1.)
Index: src/gridcheck_gfs.f90
===================================================================
--- src/gridcheck_gfs.f90	(revision a756649aa5143665ab2654ea5dd5e713950bb1a2)
+++ src/gridcheck_gfs.f90	(revision 467460af7de11cc1002100a8e56de8b8108cd235)
@@ -75,10 +75,10 @@
   real :: sizesouth,sizenorth,xauxa,pint
   real :: akm_usort(nwzmax)
-  real,parameter :: eps=0.0001
+  real,parameter :: eps=spacing(2.0_4*360.0_4) 
 
   ! NCEP GFS
   real :: pres(nwzmax), help
 
-  integer :: i179,i180,i181
+  integer :: i180
 
   ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
@@ -224,7 +224,7 @@
   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
+  if((abs(xaux1).lt.eps).and.(xaux2.ge.359.)) then ! NCEP DATA FROM 0 TO
+    xaux1=-180.0                             ! 359 DEG EAST ->
+    xaux2=-180.0+360.-360./real(nxfield)    ! TRANSFORMED TO -179
   endif                                      ! TO 180 DEG EAST
   if (xaux1.gt.180) xaux1=xaux1-360.0
@@ -305,11 +305,5 @@
 
 
-  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
+  i180=nint(180./dx)    ! 0.5 deg data
 
 
@@ -321,10 +315,10 @@
       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
+        if(ix.lt.i180) then
+          oro(i180+ix,jy)=help
+          excessoro(i180+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
+          oro(ix-i180,jy)=help
+          excessoro(ix-i180,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
         endif
       end do
@@ -339,8 +333,8 @@
       do ix=0,nxfield-1
         help=zsec4(nxfield*(ny-jy-1)+ix+1)
-        if(ix.le.i180) then
-          lsm(i179+ix,jy)=help
+        if(ix.lt.i180) then
+          lsm(i180+ix,jy)=help
         else
-          lsm(ix-i181,jy)=help
+          lsm(ix-i180,jy)=help
         endif
       end do
@@ -413,5 +407,5 @@
     write(*,*)
     write(*,*)
-    write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', &
+  write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', &
          nuvz,nwz
     write(*,*)
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision 759df5f2d94133d465ce3faa86d8b557a37b0127)
+++ src/par_mod.f90	(revision aa939a938b10003d7e7976c3f5174db7add1877e)
@@ -147,5 +147,7 @@
 ! GFS
    integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138
-   integer :: nxshift=0 ! shift not fixed for the executable 
+  ! GFS 0.25
+  ! integer,parameter :: nxmax=1441,nymax=721,nuvzmax=138,nwzmax=138,nzmax=138
+  integer :: nxshift=0 ! shift not fixed for the executable 
 
 
Index: src/readreleases.f90
===================================================================
--- src/readreleases.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/readreleases.f90	(revision 467460af7de11cc1002100a8e56de8b8108cd235)
@@ -546,7 +546,7 @@
       if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
         write(*,*) 'FLEXPART MODEL ERROR'
-        write(*,*) 'Release starts before simulation begins or ends'
+        write(*,*) 'Release starts before simulation begins or ends (1)'
         write(*,*) 'after simulation stops.'
-        write(*,*) 'Make files COMMAND and RELEASES consistent.'
+        write(*,*) 'Make files COMMAND and RELEASES consistent (1).'
         stop
       endif
@@ -561,7 +561,7 @@
       if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
         write(*,*) 'FLEXPART MODEL ERROR'
-        write(*,*) 'Release starts before simulation begins or ends'
+        write(*,*) 'Release starts before simulation begins or ends (2)'
         write(*,*) 'after simulation stops.'
-        write(*,*) 'Make files COMMAND and RELEASES consistent.'
+        write(*,*) 'Make files COMMAND and RELEASES consistent (2).'
         stop
       endif
Index: src/readwind_gfs.f90
===================================================================
--- src/readwind_gfs.f90	(revision a756649aa5143665ab2654ea5dd5e713950bb1a2)
+++ src/readwind_gfs.f90	(revision 9ca6e380ca9ee104eba37662c43b9e13b165195c)
@@ -4,46 +4,46 @@
 subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
 
-  !***********************************************************************
-  !*                                                                     *
-  !*             TRAJECTORY MODEL SUBROUTINE READWIND                    *
-  !*                                                                     *
-  !***********************************************************************
-  !*                                                                     *
-  !*             AUTHOR:      G. WOTAWA                                  *
-  !*             DATE:        1997-08-05                                 *
-  !*             LAST UPDATE: 2000-10-17, Andreas Stohl                  *
-  !*             CHANGE: 01/02/2001, Bernd C. Krueger, Variables tth and *
-  !*                     qvh (on eta coordinates) in common block        *
-  !*             CHANGE: 16/11/2005, Caroline Forster, GFS data          *
-  !*             CHANGE: 11/01/2008, Harald Sodemann, Input of GRIB1/2   *
-  !*                     data with the ECMWF grib_api library            *
-  !*             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                  *
-  !*                                                                     *
-  !***********************************************************************
-  !*                                                                     *
-  !* 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                                 *
-  !*                                                                     *
-  !***********************************************************************
+!***********************************************************************
+!*                                                                     *
+!*             TRAJECTORY MODEL SUBROUTINE READWIND                    *
+!*                                                                     *
+!***********************************************************************
+!*                                                                     *
+!*             AUTHOR:      G. WOTAWA                                  *
+!*             DATE:        1997-08-05                                 *
+!*             LAST UPDATE: 2000-10-17, Andreas Stohl                  *
+!*             CHANGE: 01/02/2001, Bernd C. Krueger, Variables tth and *
+!*                     qvh (on eta coordinates) in common block        *
+!*             CHANGE: 16/11/2005, Caroline Forster, GFS data          *
+!*             CHANGE: 11/01/2008, Harald Sodemann, Input of GRIB1/2   *
+!*                     data with the ECMWF grib_api library            *
+!*             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                  *
+!*                                                                     *
+!***********************************************************************
+!*                                                                     *
+!* 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 eccodes
@@ -53,10 +53,10 @@
   implicit none
 
-  !HSO  new parameters for grib_api
+!HSO  new parameters for grib_api
   integer :: ifile
   integer :: iret
   integer :: igrib
-  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
-  !HSO end edits
+  integer :: gribVer,parCat,parNum,typSurf,discipl,valSurf
+!HSO end edits
   real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
   real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
@@ -64,5 +64,5 @@
   integer :: ii,indj,i,j,k,n,levdiff2,ifield,iumax,iwmax
 
-  ! NCEP
+! NCEP
   integer :: numpt,numpu,numpv,numpw,numprh,numpclwch
   real :: help, temp, ew
@@ -72,14 +72,15 @@
   real :: qvh2(0:nxmax-1,0:nymax-1)
 
-  integer :: i179,i180,i181
-
-  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
-  !HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input
+  integer :: i180
+
+! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
+!HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input
 
   integer :: isec1(8),isec2(3)
+  real           :: xsec18
   real(kind=4) :: zsec4(jpunp)
   real(kind=4) :: xaux,yaux,xaux0,yaux0
+  real,parameter :: eps=spacing(2.0_4*360.0_4)
   real(kind=8) :: xauxin,yauxin
-  real,parameter :: eps=1.e-4
   real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
   real :: plev1,hlev1,ff10m,fflev1
@@ -87,5 +88,5 @@
   logical :: hflswitch,strswitch
 
-  !HSO  for grib api error messages
+!HSO  for grib api error messages
   character(len=24) :: gribErrorMsg = 'Error reading grib file'
   character(len=20) :: gribFunction = 'readwind_gfs'
@@ -100,13 +101,13 @@
 
 
-  ! OPENING OF DATA FILE (GRIB CODE)
-
-  !HSO
+! OPENING OF DATA FILE (GRIB CODE)
+
+!HSO
   call grib_open_file(ifile,path(3)(1:length(3)) &
-         //trim(wfname(indj)),'r',iret)
+       //trim(wfname(indj)),'r',iret)
   if (iret.ne.GRIB_SUCCESS) then
     goto 888   ! ERROR DETECTED
   endif
-  !turn on support for multi fields messages
+!turn on support for multi fields messages
   call grib_multi_support_on
 
@@ -118,8 +119,8 @@
   numpclwch=0
   ifield=0
-10   ifield=ifield+1
-  !
-  ! GET NEXT FIELDS
-  !
+10 ifield=ifield+1
+!
+! GET NEXT FIELDS
+!
   call grib_new_from_file(ifile,igrib,iret)
   if (iret.eq.GRIB_END_OF_FILE)  then
@@ -129,5 +130,5 @@
   endif
 
-  !first see if we read GRIB1 or GRIB2
+!first see if we read GRIB1 or GRIB2
   call grib_get_int(igrib,'editionNumber',gribVer,iret)
 !  call grib_check(iret,gribFunction,gribErrorMsg)
@@ -135,153 +136,160 @@
   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)
+!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)
+!JMA / SH: isec1(8) not evaluated any more below
+!b/c with GRIB 2 this may be a real variable
+    xsec18 = real(isec1(8))
 
   else ! GRIB Edition 2
 
-  !read the grib2 identifiers
-  call grib_get_string(igrib,'shortName',shortname,iret)
-
-  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)
-  
-!  write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname
-  !convert to grib1 identifiers
-  isec1(6)=-1
-  isec1(7)=-1
-  isec1(8)=-1
-  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T
-    isec1(6)=11          ! indicatorOfParameter
-    isec1(7)=100         ! indicatorOfTypeOfLevel
-    isec1(8)=valSurf/100 ! level, convert to hPa
-  elseif ((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.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V
-    isec1(6)=34          ! indicatorOfParameter
-    isec1(7)=100         ! indicatorOfTypeOfLevel
-    isec1(8)=valSurf/100 ! level, convert to hPa
-  elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W
-    isec1(6)=39          ! indicatorOfParameter
-    isec1(7)=100         ! indicatorOfTypeOfLevel
-    isec1(8)=valSurf/100 ! level, convert to hPa
-  elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH
-    isec1(6)=52          ! indicatorOfParameter
-    isec1(7)=100         ! indicatorOfTypeOfLevel
-    isec1(8)=valSurf/100 ! level, convert to hPa
-  elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2
-    isec1(6)=52          ! indicatorOfParameter
-    isec1(7)=105         ! indicatorOfTypeOfLevel
-    isec1(8)=2
-  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2
-    isec1(6)=11          ! indicatorOfParameter
-    isec1(7)=105         ! indicatorOfTypeOfLevel
-    isec1(8)=2
-  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10
-    isec1(6)=33          ! indicatorOfParameter
-    isec1(7)=105         ! indicatorOfTypeOfLevel
-    isec1(8)=10
-  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10
-    isec1(6)=34          ! indicatorOfParameter
-    isec1(7)=105         ! indicatorOfTypeOfLevel
-    isec1(8)=10
-  elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]:
-    isec1(6)=153         ! indicatorOfParameter
-    isec1(7)=100         ! indicatorOfTypeOfLevel
-    isec1(8)=valSurf/100 ! level, convert to hPa
-  elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
-    isec1(6)=2           ! indicatorOfParameter
-    isec1(7)=102         ! indicatorOfTypeOfLevel
-    isec1(8)=0
-  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP
-    isec1(6)=1           ! indicatorOfParameter
-    isec1(7)=1           ! indicatorOfTypeOfLevel
-    isec1(8)=0
-  elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW
-    isec1(6)=66          ! indicatorOfParameter
-    isec1(7)=1           ! indicatorOfTypeOfLevel
-    isec1(8)=0
-  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0
-    isec1(6)=11          ! indicatorOfParameter
-    isec1(7)=107         ! indicatorOfTypeOfLevel
-    isec1(8)=0.995       ! lowest sigma level
-  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0
-    isec1(6)=33          ! indicatorOfParameter
-    isec1(7)=107         ! indicatorOfTypeOfLevel
-    isec1(8)=0.995       ! lowest sigma level
-  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0
-    isec1(6)=34          ! indicatorOfParameter
-    isec1(7)=107         ! indicatorOfTypeOfLevel
-    isec1(8)=0.995       ! lowest sigma level
-  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
-  elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH
-    isec1(6)=221         ! indicatorOfParameter
-    isec1(7)=1           ! indicatorOfTypeOfLevel
-    isec1(8)=0
-  elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP
-    isec1(6)=62          ! indicatorOfParameter
-    isec1(7)=1           ! indicatorOfTypeOfLevel
-    isec1(8)=0
-  elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP
-    isec1(6)=63          ! indicatorOfParameter
-    isec1(7)=1           ! indicatorOfTypeOfLevel
-    isec1(8)=0
-  endif
+!read the grib2 identifiers
+    
+    call grib_get_string(igrib,'shortName',shortname,iret)
+
+    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)
+
+!    write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname
+!convert to grib1 identifiers
+    isec1(6:8)=-1
+    xsec18  =-1.0
+!JMA / SH: isec1(8) not evaluated any more below
+!b/c with GRIB 2 this may be a real variable
+    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T
+      isec1(6)=11          ! indicatorOfParameter
+      isec1(7)=100         ! indicatorOfTypeOfLevel
+      xsec18=valSurf/100.0 ! level, convert to hPa
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
+      isec1(6)=33          ! indicatorOfParameter
+      isec1(7)=100         ! indicatorOfTypeOfLevel
+      xsec18=valSurf/100.0 ! level, convert to hPa
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V
+      isec1(6)=34          ! indicatorOfParameter
+      isec1(7)=100         ! indicatorOfTypeOfLevel
+      xsec18=valSurf/100.0 ! level, convert to hPa
+    elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W
+      isec1(6)=39          ! indicatorOfParameter
+      isec1(7)=100         ! indicatorOfTypeOfLevel
+      xsec18=valSurf/100.0 ! level, convert to hPa
+    elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH
+      isec1(6)=52          ! indicatorOfParameter
+      isec1(7)=100         ! indicatorOfTypeOfLevel
+      xsec18=valSurf/100.0 ! level, convert to hPa
+    elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2
+      isec1(6)=52          ! indicatorOfParameter
+      isec1(7)=105         ! indicatorOfTypeOfLevel
+      xsec18=real(2)
+    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2
+      isec1(6)=11          ! indicatorOfParameter
+      isec1(7)=105         ! indicatorOfTypeOfLevel
+      xsec18=real(2)
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10
+      isec1(6)=33          ! indicatorOfParameter
+      isec1(7)=105         ! indicatorOfTypeOfLevel
+      xsec18=real(10)
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10
+      isec1(6)=34          ! indicatorOfParameter
+      isec1(7)=105         ! indicatorOfTypeOfLevel
+      xsec18=real(10)
+    elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]:
+      isec1(6)=153         ! indicatorOfParameter
+      isec1(7)=100         ! indicatorOfTypeOfLevel
+      xsec18=valSurf/100.0 ! level, convert to hPa
+    elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
+      isec1(6)=2           ! indicatorOfParameter
+      isec1(7)=102         ! indicatorOfTypeOfLevel
+      xsec18=real(0)
+    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP
+      isec1(6)=1           ! indicatorOfParameter
+      isec1(7)=1           ! indicatorOfTypeOfLevel
+      xsec18=real(0)
+    elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW
+      isec1(6)=66          ! indicatorOfParameter
+      isec1(7)=1           ! indicatorOfTypeOfLevel
+      xsec18=real(0)
+    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0
+      isec1(6)=11          ! indicatorOfParameter
+      isec1(7)=107         ! indicatorOfTypeOfLevel
+      xsec18=0.995         ! lowest sigma level
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0
+      isec1(6)=33          ! indicatorOfParameter
+      isec1(7)=107         ! indicatorOfTypeOfLevel
+      xsec18=0.995         ! lowest sigma level
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0
+      isec1(6)=34          ! indicatorOfParameter
+      isec1(7)=107         ! indicatorOfTypeOfLevel
+      xsec18=0.995         ! lowest sigma level
+    elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
+      isec1(6)=7           ! indicatorOfParameter
+      isec1(7)=1           ! indicatorOfTypeOfLevel
+      xsec18=real(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
+      xsec18=real(0)
+    elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH
+      isec1(6)=221         ! indicatorOfParameter
+      isec1(7)=1           ! indicatorOfTypeOfLevel
+      xsec18=real(0)
+    elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP
+      isec1(6)=62          ! indicatorOfParameter
+      isec1(7)=1           ! indicatorOfTypeOfLevel
+      xsec18=real(0)
+    elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP
+      isec1(6)=63          ! indicatorOfParameter
+      isec1(7)=1           ! indicatorOfTypeOfLevel
+      xsec18=real(0)
+    endif
 
   endif ! gribVer
 
   if (isec1(6).ne.-1) then
-  !  get the size and data of the values array
+!  get the size and data of the values array
     call grib_get_real4_array(igrib,'values',zsec4,iret)
-!    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_check(iret,gribFunction,gribErrorMsg)
   endif
 
   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', &
-       xauxin,iret)
-!  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
-       yauxin,iret)
-!  call grib_check(iret,gribFunction,gribErrorMsg)
-  xaux=xauxin+real(nxshift)*dx
-  yaux=yauxin
-
-  ! CHECK GRID SPECIFICATIONS
+!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', &
+         xauxin,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
+         yauxin,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    xaux=xauxin+real(nxshift)*dx
+    yaux=yauxin
+
+! CHECK GRID SPECIFICATIONS
 
     if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
     if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
-    if(xaux.eq.0.) xaux=-179.0     ! NCEP DATA
+    if(xaux.eq.0.) xaux=-180.0     ! NCEP DATA
     xaux0=xlon0
     yaux0=ylat0
@@ -290,280 +298,268 @@
     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
-  !HSO end of edits
-
-  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
+    if(abs(xaux-xaux0).gt.eps) then
+      write (*, *) xaux, xaux0
+      stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
+    endif
+    if(abs(yaux-yaux0).gt.eps) then
+      write (*, *) yaux, yaux0
+      stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
+    end if
+  endif
+!HSO end of edits
+
+  i180=nint(180./dx)
 
   if (isec1(6).ne.-1) then
 
-  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.i180) then
-          tth(i179+i,j,numpt,n)=help
-        else
-          tth(i-i181,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.i180) then
-          uuh(i179+i,j,numpu)=help
-        else
-          uuh(i-i181,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.i180) then
-          vvh(i179+i,j,numpv)=help
-        else
-          vvh(i-i181,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.i180) then
-          qvh(i179+i,j,numprh,n)=help
-        else
-          qvh(i-i181,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.i180) then
-          ps(i179+i,j,1,n)=help
-        else
-          ps(i-i181,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.i180) then
-          wwh(i179+i,j,numpw)=help
-        else
-          wwh(i-i181,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.i180) then
-          sd(i179+i,j,1,n)=help
-        else
-          sd(i-i181,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.i180) then
-          msl(i179+i,j,1,n)=help
-        else
-          msl(i-i181,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.i180) then
-          tcc(i179+i,j,1,n)=help
-        else
-          tcc(i-i181,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.i180) then
-          u10(i179+i,j,1,n)=help
-        else
-          u10(i-i181,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.i180) then
-          v10(i179+i,j,1,n)=help
-        else
-          v10(i-i181,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.i180) then
-          tt2(i179+i,j,1,n)=help
-        else
-          tt2(i-i181,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.i180) then
-          td2(i179+i,j,1,n)=help
-        else
-          td2(i-i181,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.i180) then
-          lsprec(i179+i,j,1,n)=help
-        else
-          lsprec(i-i181,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.i180) then
-          convprec(i179+i,j,1,n)=help
-        else
-          convprec(i-i181,j,1,n)=help
-        endif
-      endif
-      if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
-  ! TOPOGRAPHY
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.i180) then
-          oro(i179+i,j)=help
-          excessoro(i179+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
-        else
-          oro(i-i181,j)=help
-          excessoro(i-i181,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.i180) then
-          lsm(i179+i,j)=help
-        else
-          lsm(i-i181,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.i180) then
-          hmix(i179+i,j,1,n)=help
-        else
-          hmix(i-i181,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.i180) then
-          qvh2(i179+i,j)=help
-        else
-          qvh2(i-i181,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.i180) then
-          tlev1(i179+i,j)=help
-        else
-          tlev1(i-i181,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.i180) then
-          ulev1(i179+i,j)=help
-        else
-          ulev1(i-i181,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.i180) then
-          vlev1(i179+i,j)=help
-        else
-          vlev1(i-i181,j)=help
-        endif
-      endif
+! write (*, *) 'nxfield: ', nxfield, i180
+
+    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
+            numpt=minloc(abs(xsec18*100.0-akz),dim=1)
+          endif
+          help=zsec4(nxfield*(ny-j-1)+i+1)
+          if (help.eq.0) then
+            write (*, *) 'i, j: ', i, j
+            stop 'help == 0.0'
+          endif
+          if(i.lt.i180) then
+            tth(i180+i,j,numpt,n)=help
+          else
+            tth(i-i180,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
+            numpu=minloc(abs(xsec18*100.0-akz),dim=1)
+          endif
+          help=zsec4(nxfield*(ny-j-1)+i+1)
+          if(i.lt.i180) then
+            uuh(i180+i,j,numpu)=help
+          else
+            uuh(i-i180,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
+            numpv=minloc(abs(xsec18*100.0-akz),dim=1)
+          endif
+          help=zsec4(nxfield*(ny-j-1)+i+1)
+          if(i.lt.i180) then
+            vvh(i180+i,j,numpv)=help
+          else
+            vvh(i-i180,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
+            numprh=minloc(abs(xsec18*100.0-akz),dim=1)
+          endif
+          help=zsec4(nxfield*(ny-j-1)+i+1)
+          if(i.lt.i180) then
+            qvh(i180+i,j,numprh,n)=help
+          else
+            qvh(i-i180,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.lt.i180) then
+            ps(i180+i,j,1,n)=help
+          else
+            ps(i-i180,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)) numpw=minloc(abs(xsec18*100.0-akz),dim=1)
+          help=zsec4(nxfield*(ny-j-1)+i+1)
+          if(i.lt.i180) then
+            wwh(i180+i,j,numpw)=help
+          else
+            wwh(i-i180,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.lt.i180) then
+            sd(i180+i,j,1,n)=help
+          else
+            sd(i-i180,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.lt.i180) then
+            msl(i180+i,j,1,n)=help
+          else
+            msl(i-i180,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.lt.i180) then
+            tcc(i180+i,j,1,n)=help
+          else
+            tcc(i-i180,j,1,n)=help
+          endif
+        endif
+        if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
+             (nint(xsec18).eq.10)) then
+! 10 M U VELOCITY
+          help=zsec4(nxfield*(ny-j-1)+i+1)
+          if(i.lt.i180) then
+            u10(i180+i,j,1,n)=help
+          else
+            u10(i-i180,j,1,n)=help
+          endif
+        endif
+        if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
+             (nint(xsec18).eq.10)) then
+! 10 M V VELOCITY
+          help=zsec4(nxfield*(ny-j-1)+i+1)
+          if(i.lt.i180) then
+            v10(i180+i,j,1,n)=help
+          else
+            v10(i-i180,j,1,n)=help
+          endif
+        endif
+        if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
+             (nint(xsec18).eq.2)) then
+! 2 M TEMPERATURE
+          help=zsec4(nxfield*(ny-j-1)+i+1)
+          if(i.lt.i180) then
+            tt2(i180+i,j,1,n)=help
+          else
+            tt2(i-i180,j,1,n)=help
+          endif
+        endif
+        if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
+             (nint(xsec18).eq.2)) then
+! 2 M DEW POINT TEMPERATURE
+          help=zsec4(nxfield*(ny-j-1)+i+1)
+          if(i.lt.i180) then
+            td2(i180+i,j,1,n)=help
+          else
+            td2(i-i180,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.lt.i180) then
+            lsprec(i180+i,j,1,n)=help
+          else
+            lsprec(i-i180,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.lt.i180) then
+            convprec(i180+i,j,1,n)=help
+          else
+            convprec(i-i180,j,1,n)=help
+          endif
+        endif
+        if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
+! TOPOGRAPHY
+          help=zsec4(nxfield*(ny-j-1)+i+1)
+          if(i.lt.i180) then
+            oro(i180+i,j)=help
+            excessoro(i180+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
+          else
+            oro(i-i180,j)=help
+            excessoro(i-i180,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.lt.i180) then
+            lsm(i180+i,j)=help
+          else
+            lsm(i-i180,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.lt.i180) then
+            hmix(i180+i,j,1,n)=help
+          else
+            hmix(i-i180,j,1,n)=help
+          endif
+        endif
+        if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
+             (nint(xsec18).eq.02)) then
+! 2 M RELATIVE HUMIDITY
+          help=zsec4(nxfield*(ny-j-1)+i+1)
+          if(i.lt.i180) then
+            qvh2(i180+i,j)=help
+          else
+            qvh2(i-i180,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.lt.i180) then
+            tlev1(i180+i,j)=help
+          else
+            tlev1(i-i180,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.lt.i180) then
+            ulev1(i180+i,j)=help
+          else
+            ulev1(i-i180,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.lt.i180) then
+            vlev1(i180+i,j)=help
+          else
+            vlev1(i-i180,j)=help
+          endif
+        endif
 ! SEC & IP 12/2018 read GFS clouds
-      if((isec1(6).eq.153).and.(isec1(7).eq.100)) then  !! CLWCR  Cloud liquid water content [kg/kg] 
-         if((i.eq.0).and.(j.eq.0)) then
-            do ii=1,nuvz
-              if ((isec1(8)*100.0).eq.akz(ii)) numpclwch=ii
-            end do
-        endif
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.i180) then
-          clwch(i179+i,j,numpclwch,n)=help
-        else
-          clwch(i-i181,j,numpclwch,n)=help
-        endif
-        readclouds=.true.
-        sumclouds=.true.
-!        readclouds=.false.
-!       sumclouds=.false.
-      endif
-
-
+        if((isec1(6).eq.153).and.(isec1(7).eq.100)) then  !! CLWCR  Cloud liquid water content [kg/kg] 
+          if((i.eq.0).and.(j.eq.0)) then
+            numpclwch=minloc(abs(xsec18*100.0-akz),dim=1)
+          endif
+          help=zsec4(nxfield*(ny-j-1)+i+1)
+          if(i.lt.i180) then
+            clwch(i180+i,j,numpclwch,n)=help
+          else
+            clwch(i-i180,j,numpclwch,n)=help
+          endif
+          readclouds=.true.
+          sumclouds=.true.
+        endif
+
+
+      end do
     end do
-  end do
 
   endif
 
   if((isec1(6).eq.33).and.(isec1(7).eq.100)) then
-  ! NCEP ISOBARIC LEVELS
+! NCEP ISOBARIC LEVELS
     iumax=iumax+1
   endif
@@ -571,49 +567,31 @@
   call grib_release(igrib)
   goto 10                      !! READ NEXT LEVEL OR PARAMETER
-  !
-  ! CLOSING OF INPUT DATA FILE
-  !
-
-  !HSO close grib file
-50   continue
+!
+! CLOSING OF INPUT DATA FILE
+!
+
+!HSO close grib file
+50 continue
   call grib_close_file(ifile)
 
-  ! SENS. HEAT FLUX
+! SENS. HEAT FLUX
   sshf(:,:,1,n)=0.0     ! not available from gfs.tccz.pgrbfxx files
   hflswitch=.false.    ! Heat flux not available
-  ! SOLAR RADIATIVE FLUXES
+! SOLAR RADIATIVE FLUXES
   ssr(:,:,1,n)=0.0      ! not available from gfs.tccz.pgrbfxx files
-  ! EW SURFACE STRESS
+! EW SURFACE STRESS
   ewss=0.0         ! not available from gfs.tccz.pgrbfxx files
-  ! NS SURFACE STRESS
+! NS SURFACE STRESS
   nsss=0.0         ! not available from gfs.tccz.pgrbfxx files
   strswitch=.false.    ! stress not available
 
-  ! CONVERT TP TO LSP (GRIB2 only)
+! CONVERT TP TO LSP (GRIB2 only)
   if (gribVer.eq.2) then
-    do j=0,nymin1
-    do i=0,nxfield-1
-     if(i.le.i180) then
-     if (convprec(i179+i,j,1,n).lt.lsprec(i179+i,j,1,n)) then ! neg precip would occur
-         lsprec(i179+i,j,1,n)= &
-              lsprec(i179+i,j,1,n)-convprec(i179+i,j,1,n)
-     else
-         lsprec(i179+i,j,1,n)=0
-     endif
-     else
-     if (convprec(i-i181,j,1,n).lt.lsprec(i-i181,j,1,n)) then
-          lsprec(i-i181,j,1,n)= &
-               lsprec(i-i181,j,1,n)-convprec(i-i181,j,1,n)
-     else
-          lsprec(i-i181,j,1,n)=0
-     endif
-     endif
-    enddo
-    enddo
-  endif
-  !HSO end edits
-
-
-  ! TRANSFORM RH TO SPECIFIC HUMIDITY
+    lsprec(0:nxfield-1, 0:nymin1, 1, n) = max( 0.0 , lsprec(0:nxfield-1,0:nymin1,1,n)-convprec(0:nxfield-1,0:nymin1,1,n) )
+  endif
+!HSO end edits
+
+
+! TRANSFORM RH TO SPECIFIC HUMIDITY
 
   do j=0,ny-1
@@ -622,4 +600,8 @@
         help=qvh(i,j,k,n)
         temp=tth(i,j,k,n)
+        if (temp .eq. 0.0) then 
+          write (*, *) i, j, k, n
+          temp = 273.0
+        endif
         plev1=akm(k)+bkm(k)*ps(i,j,1,n)
         elev=ew(temp)*help/100.0
@@ -629,15 +611,19 @@
   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
+! 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)
+      help=qvh2(i,j)
+      temp=tt2(i,j,1,n)
+      if (temp .eq. 0.0) then 
+        write (*, *) i, j, n
+        temp = 273.0
+      endif
+      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
@@ -653,7 +639,7 @@
 
 
-  ! For global fields, assign the leftmost data column also to the rightmost
-  ! data column; if required, shift whole grid by nxshift grid points
-  !*************************************************************************
+! 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
@@ -691,5 +677,5 @@
   do i=0,nxmin1
     do j=0,nymin1
-  ! Convert precip. from mm/s -> mm/hour
+! Convert precip. from mm/s -> mm/hour
       convprec(i,j,1,n)=convprec(i,j,1,n)*3600.
       lsprec(i,j,1,n)=lsprec(i,j,1,n)*3600.
@@ -699,9 +685,9 @@
 
   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
-  !***************************************************************************
+!  write(*,*) 'WARNING: No flux data contained in GRIB file ',
+!    +  wfname(indj)
+
+! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
+!***************************************************************************
 
     do i=0,nxmin1
@@ -723,9 +709,9 @@
 
   return
-888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
+888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
   write(*,*) ' #### ',wfname(indj),'                    #### '
   write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
   stop 'Execution terminated'
-999   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
+999 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
   write(*,*) ' #### ',wfname(indj),'                    #### '
   write(*,*) ' #### CANNOT BE OPENED !!!                    #### '
Index: src/richardson.f90
===================================================================
--- src/richardson.f90	(revision 92fab6545ef2a8e57cadc977df94b6b510318a1d)
+++ src/richardson.f90	(revision 467460af7de11cc1002100a8e56de8b8108cd235)
@@ -141,9 +141,13 @@
     zold=z
   end do
-  k=k-1 ! ESO: make sure k <= nuvz (ticket #139)
+  ! k=k-1 ! ESO: make sure k <= nuvz (ticket #139)
 20   continue
 
   ! Determine Richardson number between the critical levels
   !********************************************************
+  ! JMA: It may happen that k >= nuvz:
+  ! JMA: In that case: k is set to nuvz 
+
+  if (k .gt. nuvz) k = nuvz ! JMA
 
   zl1=zold
