Subversion Repositories lagranto.wrf

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) print*,"Error in inq:",trim(varname)
         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     ------------------------------------------------------------