Index: branches/petra/src/readwind_nests.f90
===================================================================
--- branches/petra/src/readwind_nests.f90	(revision 36)
+++ branches/petra/src/readwind_nests.f90	(revision 37)
@@ -34,9 +34,13 @@
   !                                                                            *
   !*****************************************************************************
-  !  Changes, Bernd C. Krueger, Feb. 2001:                                     *
+  !  CHANGES:
+  !
+  !  2/2001, Bernc C. Krueger:                                              *
   !        Variables tthn and qvhn (on eta coordinates) in common block        *
-  !  CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with ECMWF grib_api    *
-  !  CHANGE: 03/12/2008, Harald Sodemann, update to f90 with ECMWF grib_api    *
-  ! 2/2015, PS: add missing paramter iret in call to grib subr
+  !  11/01/2008, Harald Sodemann, GRIB1/2 input with ECMWF grib_api            *
+  !  03/12/2008, Harald Sodemann, update to f90 with ECMWF grib_api            *
+  !  2/2015, PS: add missing paramter iret in call to grib subr
+  !  3/2015, PS: add verbosity output
+  !  3/2015, PS: add GRIB2 coding for deta/dt dp/deta used by ZAMG
   !
   !*****************************************************************************
@@ -84,4 +88,6 @@
   character(len=24) :: gribErrorMsg = 'Error reading grib file'
   character(len=20) :: gribFunction = 'readwind_nests'
+  
+  call verboutput(verbosity,1,'      entered readwind_nests')
 
   do l=1,numbnests
@@ -100,5 +106,8 @@
   !
 
-5   call grib_open_file(ifile,path(numpath+2*(l-1)+1) &
+5 continue
+  call verboutput(verbosity,1,'        call grib_open_file '//&
+    path(numpath+2*(l-1)+1)(1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,indj)))
+  call grib_open_file(ifile,path(numpath+2*(l-1)+1) &
          (1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,indj)),'r')
   if (iret.ne.GRIB_SUCCESS) then
@@ -127,101 +136,106 @@
   if (gribVer.eq.1) then ! GRIB Edition 1
 
-  !print*,'GRiB Edition 1'
-  !read the grib2 identifiers
-  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'level',isec1(8),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-
-  !change code for etadot to code for omega
-  if (isec1(6).eq.77) then
-    isec1(6)=135
-  endif
-
-  conversion_factor=1.
-
+    !print*,'GRiB Edition 1'
+    
+    ! read the grib2 identifiers
+    
+    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'level',isec1(8),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+
+    !change code for etadot to code for omega
+    if (isec1(6).eq.77) then
+      isec1(6)=135
+    endif
+
+    conversion_factor=1.
 
   else
 
-  !print*,'GRiB Edition 2'
-  !read the grib2 identifiers
-  call grib_get_int(igrib,'discipline',discipl,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'parameterCategory',parCat,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'parameterNumber',parNum,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'level',valSurf,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new readwind.f90
-  call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new readwind.f90
-
-  !print*,discipl,parCat,parNum,typSurf,valSurf
-
-  !convert to grib1 identifiers
-  isec1(6)=-1
-  isec1(7)=-1
-  isec1(8)=-1
-  isec1(8)=valSurf     ! level
-   conversion_factor=1.
-  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
-    isec1(6)=130         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
-    isec1(6)=131         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
-    isec1(6)=132         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
-    isec1(6)=133         ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
-    isec1(6)=134         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot !
-    isec1(6)=135         ! indicatorOfParameter
-  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added by mc to make it consisitent with new readwind.f90
-    isec1(6)=135         ! indicatorOfParameter    !added by mc to make it consisitent with new readwind.f90
-  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
-    isec1(6)=151         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
-    isec1(6)=165         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
-    isec1(6)=166         ! indicatorOfParameter
-  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
-    isec1(6)=167         ! indicatorOfParameter
-  elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
-    isec1(6)=168         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
-    isec1(6)=141         ! indicatorOfParameter
-    conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
-  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consisitent with new readwind.f90
-    isec1(6)=164         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consisitent with new readwind.f90
-    isec1(6)=142         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
-    isec1(6)=143         ! indicatorOfParameter
-    conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
-  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
-    isec1(6)=146         ! indicatorOfParameter
-  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
-    isec1(6)=176         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consisitent with new readwind.f90
-    isec1(6)=180         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consisitent with new readwind.f90
-    isec1(6)=181         ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
-    isec1(6)=129         ! indicatorOfParameter
-   elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consisitent with new readwind.f90
-    isec1(6)=160         ! indicatorOfParameter
-  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
-       (typSurf.eq.1)) then ! LSM
-    isec1(6)=172         ! indicatorOfParameter
-  else
-    print*,'***WARNING: undefined GRiB2 message found!',discipl, &
-         parCat,parNum,typSurf
-  endif
-  if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consisitent with new readwind.f90
-    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)  !
-!    stop
-  endif
+    !print*,'GRiB Edition 2'
+
+    ! read the grib2 identifiers
+
+    call grib_get_int(igrib,'discipline',discipl,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'parameterCategory',parCat,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'parameterNumber',parNum,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'level',valSurf,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consistent with new readwind.f90
+    call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consistent with new readwind.f90
+
+    !print*,discipl,parCat,parNum,typSurf,valSurf
+
+    !convert to grib1 identifiers
+    isec1(6)=-1
+    isec1(7)=-1
+    isec1(8)=-1
+    isec1(8)=valSurf     ! level
+     conversion_factor=1.
+    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
+      isec1(6)=130         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
+      isec1(6)=131         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
+      isec1(6)=132         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
+      isec1(6)=133         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
+      isec1(6)=134         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually detadtdpdeta
+      isec1(6)=135         ! indicatorOfParameter
+    elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually detadtdpdeta
+      isec1(6)=135         ! indicatorOfParameter !added by mc to make it consistent with new readwind.f90
+    elseif ((parCat.eq.2).and.(parNum.eq.8)) then ! Omega, actually detadtdpdeta
+      isec1(6)=135         ! indicatorOfParameter !added by mc to make it consistent with new readwind.f90
+    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
+      isec1(6)=151         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
+      isec1(6)=165         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
+      isec1(6)=166         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
+      isec1(6)=167         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
+      isec1(6)=168         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
+      isec1(6)=141         ! indicatorOfParameter
+      conversion_factor=1000. !added by mc to make it consistent with new readwind.f90
+    elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consistent with new readwind.f90
+      isec1(6)=164         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consistent with new readwind.f90
+      isec1(6)=142         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
+      isec1(6)=143         ! indicatorOfParameter
+      conversion_factor=1000. !added by mc to make it consistent with new readwind.f90
+    elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
+      isec1(6)=146         ! indicatorOfParameter
+    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
+      isec1(6)=176         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consistent with new readwind.f90
+      isec1(6)=180         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consistent with new readwind.f90
+      isec1(6)=181         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
+      isec1(6)=129         ! indicatorOfParameter
+     elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consistent with new readwind.f90
+      isec1(6)=160         ! indicatorOfParameter
+    elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
+         (typSurf.eq.1)) then ! LSM
+      isec1(6)=172         ! indicatorOfParameter
+    else
+      print*,'***WARNING: undefined GRiB2 message found!',discipl, &
+           parCat,parNum,typSurf
+    endif
+    if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consistent with new readwind.f90
+      write(*,*) 'parId',parId, 'isec1(6)',isec1(6)  !
+  !    stop
+    endif
 
   endif
@@ -246,13 +260,13 @@
   ! CHECK GRID SPECIFICATIONS
   if(isec2(2).ne.nxn(l)) stop &
-  'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
+    'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
   if(isec2(3).ne.nyn(l)) stop &
-  'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
-  if(isec2(12)/2-1.ne.nlev_ec) stop 'READWIND: VERTICAL DISCRET&
-  IZATION NOT CONSISTENT FOR A NESTING LEVEL'
+    'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
+  if(isec2(12)/2-1.ne.nlev_ec) stop &
+    'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT FOR A NESTING LEVEL'
   endif ! ifield
 
   !HSO  get the second part of the grid dimensions only from GRiB1 messages
- if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consisitent with new readwind.f90
+ if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consistent with new readwind.f90
     call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
          xauxin,iret)
@@ -294,5 +308,5 @@
              zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
         if(isec1(6).eq.141) sdn(i,j,1,n,l)= &!! SNOW DEPTH
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90!
+             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consistent with new readwind.f90!
         if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS.
              zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
@@ -312,5 +326,5 @@
         endif
         if(isec1(6).eq.143) then                         !! CONVECTIVE PREC.
-          convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90
+          convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consistent with new readwind.f90
           if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0.
         endif
@@ -421,5 +435,8 @@
   end do
 
+  call verboutput(verbosity,1,'      leaving readwind_nests')
+
   return
+
 888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
   write(*,*) ' #### ',wfnamen(l,indj),' FOR NESTING LEVEL  #### '
@@ -431,4 +448,5 @@
   write(*,*) ' #### ',wfnamen(l,indj),'                    #### '
   write(*,*) ' #### CANNOT BE OPENED FOR NESTING LEVEL ',l,'####'
+  call verboutput(verbosity,1,'      leaving readwind')
 
 end subroutine readwind_nests
