Subversion Repositories lagranto.20cr

Rev

Blame | Last modification | View Log | Download | RSS feed 

c     ************************************************************
c     * This package provides input routines to read the wind    *
c     * and other fields from IVE necdf files. The routines are  *
c     *                                                          *
c     * 1) input_open  : to open a data file                     *
c     * 2) input_grid  : to read the grid information, including *
c     *                  the vertical levels                     *
c     * 3) input_wind  : to read the wind components             *
c     * 4) input_close : to close an input file                  *
c     *                                                          *
c     * The file is characterised by an filename <filename> and  *
c     * a file identifier <fid>. The horizontal grid is given by *
c     * <xmin,xmax,ymin,ymax,dx,dy,nx,ny> where the pole of the  *
c     * rotated grid is given by <pollon,pollat>. The vertical   *
c     * grid is characterised by the surface pressure <ps> and   *
c     * the pressure at staggered <slev> and unstaggered <ulev>  *
c     * levels. The number of levels is given by <nz>. Finally,  *
c     * the retrieval of the wind <field> with name <fieldname>  *
c     * is characterised by a <time> and a missing data value    *
c     * <mdv>.                                                   *
c     *                                                          *
c     * Author: Michael Sprenger, Autumn 2008                    *
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. 

      use netcdf

      implicit none

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

c     Declaration of auxiliary variables
      integer      ierr

      ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)

      end


c     ------------------------------------------------------------
c     Read information about the grid
c     ------------------------------------------------------------
      
      subroutine input_grid 
     >                   (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >                    time,pollon,pollat,p3,ps,nz,ak,bk,stagz,
     >                    timecheck)

      use netcdf

c     Read grid information at <time> from file with identifier <fid>. 
c     The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>
c     with pole position at <pollon,pollat> and grid dimension <nx,ny>.
c     The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either
c     on the staggered or unstaggered grid (with <stagz> as the flag).
c     The surface pressure is given in <ps(nx,ny)>. If <fid> is negative, 
c     only the grid dimensions and grid parameters (xmin...pollat,nz) are 
c     determined and returned (this is needed for dynamical allocation of 
c     memory).

      implicit none

c     Declaration of subroutine parameters 
      integer      fid                 ! File identifier
      real         xmin,xmax,ymin,ymax ! Domain size
      real         dx,dy               ! Horizontal resolution
      integer      nx,ny,nz            ! Grid dimensions
      real         pollon,pollat       ! Longitude and latitude of pole
      real         p3(nx,ny,nz)        ! Staggered levels
      real         ps(nx,ny)           ! Surface pressure
      real         time                ! Time of the grid information
      real         ak(nz),bk(nz)       ! Ak and Bk for layers or levels
      real         stagz               ! Vertical staggering (0 or -0.5)
      character*80 fieldname           ! Variable from which to take grid info
      character*80 timecheck           ! Either 'yes' or 'no'

c     Numerical epsilon
      real         eps
      parameter    (eps=0.0001)

c     Auxiliary varaibles
      integer      ierr       
      integer      i,j,k
      character*80 varname
      character*80 newname
      real         lon(1000),lat(1000),lev(1000),lev2(1000)
      integer      varid,lonid,latid,levid
      integer      nz2
      real         tmp(nx-1,ny)
      integer      indx
      real         lon1
      integer      itmp

c     Inquire dimensions and grid constants if <fid> is negative
      if (fid.lt.0) then

c       Longitude
        ierr = nf90_inq_dimid(-fid,'lon', lonid)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = nf90_inquire_dimension(-fid, lonid, len = nx)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = NF90_INQ_VARID(-fid,'lon',varid)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = NF90_GET_VAR(-fid,varid,lon(1:nx))
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

c       Latitude
        ierr = nf90_inq_dimid(-fid,'lat', latid)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = nf90_inquire_dimension(-fid, latid, len = ny)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = NF90_INQ_VARID(-fid,'lat',varid)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = NF90_GET_VAR(-fid,varid,lat(1:ny))
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

c       Set grid parameters and compare to expected setting
        xmin   = lon(1)
        xmax   = lon(nx)
        ymin   = lat(ny)
        ymax   = lat(1)
        dx     = (xmax-xmin)/real(nx-1)
        dy     = (ymax-ymin)/real(ny-1)
        if (
     >       ( nx              .ne.180 ).or.
     >       ( ny              .ne.91  ).or.
     >       ( abs(xmin -   0.).gt.eps ).or.
     >       ( abs(xmax - 358.).gt.eps ).or.
     >       ( abs(ymin +  90.).gt.eps ).or.
     >       ( abs(ymax -  90.).gt.eps ).or.
     >       ( abs(dx   -   2.).gt.eps ).or.
     >       ( abs(dy   -   2.).gt.eps ) )
     >  then
          print*,' ERROR: grid does not agree with expectation.. Stop'
          stop
        endif

c       Vertical levels
        ierr = nf90_inq_dimid(-fid,'level', levid)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = nf90_inquire_dimension(-fid, levid, len = nz)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = NF90_INQ_VARID(-fid,'level',varid)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = NF90_GET_VAR(-fid,varid,lev(1:nz))
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

c       Get vertical levels for omega and check consistence
        ierr = nf90_inq_dimid(-fid,'level_2', levid)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = nf90_inquire_dimension(-fid, levid, len = nz2)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = NF90_INQ_VARID(-fid,'level_2',varid)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = NF90_GET_VAR(-fid,varid,lev2(1:nz2))
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

        if (nz.gt.nz2 ) then
            itmp = nz
            nz   = nz2
            nz2  = itmp
        endif

        if ( ( nz.ne.19 ).or.(nz2.ne.24) ) then
          print*,' ERROR: grid inconsitence... level vs level_2'
          print*,'nz,nz2 ',nz,nz2
          stop
        endif
        do i=1,nz
            if ( abs( lev2(i)-lev(i) ).gt.eps ) then
                print*,' ERROR: grid inconsitence... level vs level_2'
                print*,i,lev2(i),lev(i)
                stop
            endif
        enddo

c       Set the final (expected) parameters, including closing and shifting
        nx     = 181
        ny     = 91
        nz     = 19
        xmin   = -180.
        xmax   = 180.
        ymin   = -90.
        ymax   = 90.
        dx     = 2.
        dy     = 2.
        pollon = 0.
        pollat = 90.
        stagz  = 0.

c     Get non-constant grid parameters (surface pressure and vertical grid)
      else

c       Set 3D pressure
        ierr = NF90_INQ_VARID(fid,'level_2',varid)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = NF90_GET_VAR(fid,varid,lev(1:nz))
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

        do i=1,nx
          do j=1,ny
            do k=1,nz
              p3(i,j,k) = lev(k)
            enddo
          enddo
        enddo

c       Set ak, bk
        do k=1,nz
           ak(k) = lev(k)
           bk(k) = 0.
        enddo

c       Read surface pressure (close and shift/swap domain; unit conversion)
        varname = 'pres'
        ierr = NF90_INQ_VARID(fid,varname,varid)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = NF90_GET_VAR(fid,varid,tmp(1:(nx-1),1:ny))
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

        do i=1,nx-1
            lon1 = xmin + real(i-1)*dx
            if ( lon1.lt.-eps ) lon1 = lon1 + 360.
            indx = nint( lon1 / dx + 1. )
            do j=1,ny
                ps(i,j) = 0.01 * tmp(indx,ny-j+1)
            enddo
        enddo
        do j=1,ny
            ps(nx,j) = ps(1,j)
        enddo

      endif
      
      end

c     ------------------------------------------------------------
c     Read wind and other met field
c     ------------------------------------------------------------

      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
     >                       timecheck)

      use netcdf

c     Read the wind component <fieldname> from the file with identifier
c     <fid> and save it in the 3d array <field>. The vertical staggering 
c     information is provided in <stagz> and gives the reference to either
c     the layer or level field from <input_grid>. A consistency check is
c     performed to have an agreement with the grid specified by <xmin,xmax,
c     ymin,ymax,dx,dy,nx,ny,nz>.

      implicit none

c     Declaration of variables and parameters
      integer      fid                 ! File identifier
      character*80 fieldname           ! Name of the wind field
      integer      nx,ny,nz            ! Dimension of fields
      real         field(nx,ny,nz)     ! 3d wind field
      real         stagz               ! Staggering in the z direction
      real         mdv                 ! Missing data flag
      real         xmin,xmax,ymin,ymax ! Domain size
      real         dx,dy               ! Horizontal resolution
      real         time                ! Time
      character*80 timecheck           ! Either 'yes' or 'no'

c     Netcdf variables
      integer      ierr
      character*80 varname
      character*80 newname

c     Numerical epsilon
      real         eps
      parameter    (eps=0.0001)

c     Auxiliary variables
      integer      i,j,k
      real         lev(1000)
      integer      varid
      real         tmp(nx-1,ny,nz)
      real         lon(nx-1)
      integer      indx
      real         lon1
      integer      is2d

c     Set the correct fieldname
      newname = fieldname
      if ( fieldname.eq.'PLEV'  ) newname='P'
      if ( fieldname.eq.'PLAY'  ) newname='P'
      if ( fieldname.eq.'P'     ) newname='P'
      if ( fieldname.eq.'OMEGA' ) newname='omega'
      if ( fieldname.eq.'PS'    ) newname='pres'
      if ( fieldname.eq.'U'     ) newname='uwnd'
      if ( fieldname.eq.'V'     ) newname='vwnd'

c     Get flag for 2D field
      is2d = 0
      if ( fieldname.eq.'pres' ) is2d = 1

c     Get 3D pressure
      if ( fieldname.eq.'P' ) then

        varname = 'level'
        ierr = NF90_INQ_VARID(fid,varname,varid)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = NF90_GET_VAR(fid,varid,lev(1:nz))
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

        do i=1,nx
          do j=1,ny
            do k=1,nz
              field(i,j,k) = lev(k)
            enddo
          enddo
        enddo

        mdv  = -9.96921e+36

c     Get 3D field (close and shift/swap domain)
      elseif ( is2d.eq.0 ) then

        varname = newname
        ierr = NF90_INQ_VARID(fid,varname,varid)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = NF90_GET_VAR(fid,varid,tmp(1:(nx-1),1:ny,1:nz))
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr  = nf90_get_att(fid, varid, "_FillValue", mdv)
        IF (ierr /= nf90_NoErr) THEN
            mdv  = -9.96921e+36
            ierr = nf90_NoErr
        ENDIF

        do i=1,nx-1
          lon1 = xmin + real(i-1)*dx
          if ( lon1.lt.-eps ) lon1 = lon1 + 360.
          indx = nint( lon1 / dx + 1. )
          do k=1,nz
            do j=1,ny
              field(i,j,k) = tmp(indx,ny-j+1,k)
            enddo
            field(nx,j,k) = field(1,j,k)
          enddo
        enddo

c     Get 2D field (close and shift/swap domain)
      elseif ( is2d.eq.1 ) then

        varname = newname
        ierr = NF90_INQ_VARID(fid,varname,varid)
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr = NF90_GET_VAR(fid,varid,tmp(1:(nx-1),1:ny,1:1))
        IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
        ierr  = nf90_get_att(fid, varid, "_FillValue", mdv)
        IF (ierr /= nf90_NoErr) THEN
            mdv  = -9.96921e+36
            ierr = nf90_NoErr
        ENDIF

        do i=1,nx-1
          lon1 = xmin + real(i-1)*dx
          if ( lon1.lt.-eps ) lon1 = lon1 + 360.
          indx = nint( lon1 / dx + 1. )
          do j=1,ny
            field(i,j,1) = tmp(indx,ny-j+1,1)
            do k=2,nz
              field(i,j,k) = field(i,j,k-1)
            enddo
            field(nx,j,k) = field(1,j,k)
          enddo
        enddo

      endif

      end

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

      subroutine input_close(fid)

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

      use netcdf

      implicit none

c     Declaration of subroutine parameters
      integer fid

c     Auxiliary variables
      integer ierr

      ierr = NF90_CLOSE(fid)
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)


 
      end