c     ************************************************************
c     * This package provides input routines to read if grid     *
c     * informations or field from WRF output			 *
c     *		                                                 *
c     * Author: Sebastian Schemm, ETH, 2013 			 *
c     * V.1.0 only support for ideal channel implemented         *
c     *      for input_grid_wrf				         *
c     *								 *
c     ************************************************************

c     ------------------------------------------------------------
c     Open input file
c     ------------------------------------------------------------

      subroutine input_open (fid,filename)

c     Open the input file with filename <filename> and return the
c     file identifier <fid> for further reference. 

      implicit none

      include 'netcdf.inc'

c     Declaration of subroutine parameters
      integer      fid              ! File identifier
      character*80 filename         ! Filename

c     Declaration of auxiliary variables
      integer      ierr

c     Open IVE netcdf file
      fid = ncopn(filename,NCWRITE,ierr)
      if (ierr.ne.0) goto 900
      
c     Exception handling
      return
      
 900  print*,'Cannot open input file ',trim(filename)
      stop

      end

c     ------------------------------------------------------------
c     Read information about the grid
c     ------------------------------------------------------------
	
	  subroutine input_grid_wrf(fid,
     >				  xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz)

  	  implicit none

      include 'netcdf.inc'

c     Declaration of subroutine parameters
      integer      fid                 ! File identifier
	  integer	   latid,lonid,levid
      real         xmin,xmax,ymin,ymax ! Domain size
      real         dx,dy               ! Horizontal resolution
      integer      nx,ny,nz            ! Grid dimensions
	  integer	   status


C	  get grid info
	  status = NF_INQ_DIMID(fid, 'south_north', LATID)
	  status = NF_INQ_DIMID(fid, 'west_east',   LONID)
	  status = NF_INQ_DIMID(fid, 'bottom_top',  LEVID)
	  if (status .ne. NF_NOERR) print*,"Error in reading grid atts"

	  status = NF_INQ_DIMLEN(fid, LONID, nx)
	  status = NF_INQ_DIMLEN(fid, LATID, ny)
	  status = NF_INQ_DIMLEN(fid, LEVID, nz)
	  if (status .ne. NF_NOERR) print*,"Error in reading grid size"

	  status = NF_GET_ATT_REAL (fid, NF_GLOBAL, 'DX', dx)
	  status = NF_GET_ATT_REAL (fid, NF_GLOBAL, 'DY', dy)
	  if (status .ne. NF_NOERR) print*,"Error in reading dx, dy"

c	  setup a pseudo grid
      xmin = 0.0
      ymin = 0.0
	  xmax = xmin+real(nx-1)*dx
      ymax = ymin+real(ny-1)*dy

	  end

c     ------------------------------------------------------------
c     Read variables
c     ------------------------------------------------------------

      subroutine input_var_wrf(fid,varname,field,n1,n2,n3)

      implicit none

      include 'netcdf.inc'

c     Declaration of subroutine parameters
      integer		  fid, varid
      integer      	  n1,n2,n3
      character*80	  varname
      real            field(n1,n2,n3)

c     Parmeters
      real            gearth         ! Earth's acceleration
      parameter       (gearth=9.81)
      real            addtheta       ! Offset for potential temperature
      parameter       (addtheta=300.)

c     Auxiliary variables
      real,allocatable,dimension (:,:,:)  :: temp,temp1,temp2
      integer         ndims
      integer         is2d
      integer		  levid,lonid,latid
      integer      	  nx,ny,nz
      integer		  status,i,j,k
      character*2	  stag

c     ------------ Get grid info ---------------------------------------
      status = NF_INQ_DIMID(fid, 'south_north', LATID)
      status = NF_INQ_DIMID(fid, 'west_east',   LONID)
      status = NF_INQ_DIMID(fid, 'bottom_top',  LEVID)
      if (status .ne. NF_NOERR) print*,"Error in reading grid atts"

      status = NF_INQ_DIMLEN(fid, LONID, nx)
      status = NF_INQ_DIMLEN(fid, LATID, ny)
      status = NF_INQ_DIMLEN(fid, LEVID, nz)
      if (status .ne. NF_NOERR) print*,"Error in reading grid size"

c     ------------ Allocate memory and set staggering mode -------------
      if (trim(varname).eq.'U') then
          stag = 'X'
          is2d = 0

      elseif (trim(varname).eq.'V') then
          stag = 'Y'
          is2d = 0

      elseif (trim(varname).eq.'W') then
          stag = 'Z'
          is2d = 0

      elseif ( (trim(varname).eq.'geopot').or.
     >         (trim(varname).eq.'Z'     ) )
     >then
          stag = 'Z'
          is2d = 0

      elseif ( (trim(varname).eq.'geopots').or.
     >         (trim(varname).eq.'ZB'     ) )
     >then
          stag = 'Z'
          is2d = 1

      elseif ( (trim(varname).eq.'pressure').or.
     >         (trim(varname).eq.'P'       ) )
     >then
	      stag = 'nil'
          is2d = 0

	  else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then
          stag = 'nil'
          is2d = 0

      else
         status = NF_INQ_VARID (fid, trim(varname), varid)
	     status = NF_GET_ATT_TEXT(fid,varid,'stagger',stag)
	     if (status .ne. NF_NOERR) then
             print*,'Error in inq:',trim(varname)
         endif
         status = NF_INQ_VARNDIMS (fid, varid, ndims)
         is2d = 0
         if ( ndims.eq.3 ) is2d = 1

      endif

      if ( stag.eq.'X' ) then
          allocate( temp (nx+1,ny,nz) )
          allocate( temp1(nx+1,ny,nz) )
          allocate( temp2(nx+1,ny,nz) )

      elseif ( stag.eq.'Y' ) then
          allocate( temp (nx,ny+1,nz) )
          allocate( temp1(nx,ny+1,nz) )
          allocate( temp2(nx,ny+1,nz) )

      elseif ( stag.eq.'Z' ) then
          allocate( temp (nx,ny,nz+1) )
          allocate( temp1(nx,ny,nz+1) )
          allocate( temp2(nx,ny,nz+1) )

      else
          allocate( temp (nx,ny,nz) )
          allocate( temp1(nx,ny,nz) )
          allocate( temp2(nx,ny,nz) )

      endif

c	  ------------ Read data ------------------------------------------

c	  Zonal wind : temp(nx+1,ny,nk)
	  if (trim(varname).eq.'U') then
                     
	     status = NF_INQ_VARID (fid, 'U', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq U"
	     status = NF_GET_VAR_REAL (fid, varid, temp)
	     if (status .ne. NF_NOERR) print*,"Error in reading U"

c     Meridional wind : temp(nx,ny+1,nk)
	  else if (trim(varname).eq.'V') then
                                             
	     status = NF_INQ_VARID (fid, 'V', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq V"
	     status = NF_GET_VAR_REAL (fid, varid, temp)
	     if (status .ne. NF_NOERR) print*,"Error in reading V"

c	  Vertical wind : temp(nx,ny,nz+1)
	  else if (trim(varname).eq.'W') then
                                            
	     status = NF_INQ_VARID (fid, 'W', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq W"
	     status = NF_GET_VAR_REAL (fid, varid, temp)
	     if (status .ne. NF_NOERR) print*,"Error in reading U"

c	  Geopotential height (base + perturbation) : temp(nx,ny,nz+1)
	  else if ( (trim(varname).eq.'geopot').or.
     >          (trim(varname).eq.'Z'     ) )
     >then
                                              
	     status = NF_INQ_VARID (fid, 'PHB', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq geopot"
	     status = NF_GET_VAR_REAL (fid, varid, temp1)
	     if (status .ne. NF_NOERR) print*,"Error in reading geopot"

	     status = NF_INQ_VARID (fid, 'PH', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq pgeopot"
	     status = NF_GET_VAR_REAL (fid, varid, temp2)
	     if (status .ne. NF_NOERR) print*,"Error in reading pgeopot"

	     do k = 1, nz+1
	       do j = 1, ny
	         do i = 1, nx
	           temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)
	         enddo
	       enddo
	     enddo

c	  surface geopotential: temp(nx,ny,nz+1)
	  else if ( (trim(varname).eq.'geopots').or.
     >          (trim(varname).eq.'ZB'     ) )
     >then
                                            
	     status = NF_INQ_VARID (fid, 'PHB', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq sgeopot"
	     status = NF_GET_VAR_REAL (fid, varid, temp1)
	     if (status .ne. NF_NOERR) print*,"Error in reading sgeopot"

	     status = NF_INQ_VARID (fid, 'PH', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq pgeopot"
	     status = NF_GET_VAR_REAL (fid, varid, temp2)
	     if (status .ne. NF_NOERR) print*,"Error in reading pgeopot"

	     do k = 1, nz+1
	       do j = 1, ny
	         do i = 1, nx
	           temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)
	         enddo
	       enddo
	     enddo

c	  Pressure (base + perturbation) : temp(nx,ny,nz)
	  elseif ( (trim(varname).eq.'pressure').or.
     >         (trim(varname).eq.'P'       ) )
     >then
                               
	      status = NF_INQ_VARID (fid, 'PB', varid)
	      if (status .ne. NF_NOERR) print*,"Error in inq pb"
	      status = NF_GET_VAR_REAL (fid, varid, temp1)
	      if (status .ne. NF_NOERR) print*,"Error in reading pb"

       	  status = NF_INQ_VARID (fid, 'P', varid)
	      if (status .ne. NF_NOERR) print*,"Error in inq p"
	      status = NF_GET_VAR_REAL (fid, varid, temp2)
	      if (status .ne. NF_NOERR) print*,"Error in reading p"

	      do k = 1, nz
	        do j = 1, ny
	          do i = 1, nx
	            temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)
	          enddo
	        enddo
	      enddo
	
c     Potential temperature: temp(nx,ny,nz)
	  else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then
                    
	      status = NF_INQ_VARID (fid, 'T', varid)
	      if (status .ne. NF_NOERR) print*,"Error in inq T"
	      status = NF_GET_VAR_REAL (fid, varid, temp)
	      if (status .ne. NF_NOERR) print*,"Error in reading T"

c	  Any other field (2D + 3D): temp(:,:,:), depending on staggering
	  else

	      status = NF_INQ_VARID (fid, trim(varname), varid)
	      if (status .ne. NF_NOERR) then
	         print*,"Error in inq:",trim(varname)
	 	  endif
	 	  status = NF_GET_VAR_REAL (fid, varid, temp)
	 	  if (status .ne. NF_NOERR) then
	 	     print*,"Error in reading:",trim(varname)
          endif

      endif

c     ------------ Destaggering in X, Y and Z direction ---------------
      if (trim(stag).eq.'X') then
         do k = 1, n3
	  	 do j = 1, n2
	   	 do i = 1, n1
	    	 field(i,j,k) = 0.5*(temp(i,j,k) + temp(i+1,j,k))
	   	 enddo
	  	 enddo
	 	 enddo

	  elseif (trim(stag).eq.'Y') then
	 	 do k = 1, n3
	  	 do j = 1, n2
	   	 do i = 1, n1
	    	 field(i,j,k) = 0.5*(temp(i,j,k) + temp(i,j+1,k))
	   	 enddo
	  	 enddo
	 	 enddo

	  elseif (trim(stag).eq.'Z') then
	     do k = 1, n3
	  	 do j = 1, n2
	   	 do i = 1, n1
	    	 field(i,j,k) = 0.5*(temp(i,j,k) + temp(i,j,k+1))
	   	 enddo
	  	 enddo
	 	 enddo

	  else
	  	 do k = 1, n3
	  	 do j = 1, n2
	   	 do i = 1, n1
	    	 field(i,j,k) = temp(i,j,k)
	   	 enddo
	  	 enddo
	 	 enddo

	  endif

c     ---------- Change units -----------------------------------------

c	  Pressure Pa -> hPa
	  if (trim(varname).eq.'pressure' .or. trim(varname).eq.'P') then
	    do k = 1, n3
	      do j = 1, n2
	        do i = 1, n1
	          field(i,j,k)=0.01 * field(i,j,k)
	        enddo
	      enddo
	    enddo

c     Potential temperature + 300
      else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then
	    do k = 1, n3
	      do j = 1, n2
	        do i = 1, n1
	          field(i,j,k)=field(i,j,k) + addtheta
	        enddo
	      enddo
	    enddo

c     Geopotential -> Geopotential Height (3D + Surface)
	  else if ( (trim(varname).eq.'geopot' ).or.
     >          (trim(varname).eq.'Z'      ).or.
     >          (trim(varname).eq.'geopots').or.
     >          (trim(varname).eq.'ZB'     ) )
     >then
     	do k = 1, n3
	      do j = 1, n2
	        do i = 1, n1
	          field(i,j,k)=field(i,j,k)/gearth
	        enddo
	      enddo
	    enddo

	  endif

c     ---------- Copy lowest level for 2d tracing ---------------------
      if ( is2d.eq.1 ) then
         do i=1,n1
           do j=1,n2
             do k=1,n3
                field(i,j,k) = field(i,j,1)
             enddo
           enddo
         enddo
      endif


	  end

c     ------------------------------------------------------------
c     Close input file
c     ------------------------------------------------------------

      subroutine input_close(fid)

c     Close the input file with file identifier <fid>.

      implicit none

      include 'netcdf.inc'

c     Declaration of subroutine parameters
      integer fid

c     Auxiliary variables
      integer ierr

c     Close file
      call ncclos(fid,ierr)
 
      end

c     ------------------------------------------------------------
c     ------------------------------------------------------------	
