Subversion Repositories lagranto.icon

Rev

Rev 8 | Blame | Compare with Previous | 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 height <zb> and     *
c     * the model level height (given on the W grid). The number *
c     * of levels is given by <nz>. Finally, the retrieval of the* 
c     * wind <field> with name <fieldname>  is characterised by  *
c     a <time> and a missing data value <mdv>.                   *
c     *                                                          *
c     * Author: Michael Sprenger, Spring 2011/2012               *
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

c     Open netcdf file
      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,polgam,z3,zb,nz,stagz,timecheck)

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 <z3(nx,ny,nz)> gives the vertical coordinates, either
c     on the staggered or unstaggered grid (with <stagz> as the flag).
c     The surface height is given in <zb(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).

      use netcdf
      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,polgam ! Longitude and latitude of pole
      real         z3(nx,ny,nz)         ! Staggered levels
      real         zb(nx,ny)            ! Surface pressure
      real         time                 ! Time of the grid information
      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 and physical parameters
      real          eps                ! Numerical epsilon
      parameter    (eps=0.001)

c     Name of the constants file
      character*80  constfile
      parameter     (constfile = 'ICONCONST')

c     netCDF variables
      real         lon(5000)
      real         lat(5000)

c     Auxiliary variables
      integer      ierr   
      integer      varid
      integer      dimid
      integer      dimids(100)
      character*80 dimname(100)   
      integer      ndim
      integer      i,j,k
      real         tmp3(nx,ny,nz)
      real         stag3(nx,ny,nz+1)
      integer      cdfid 
       
c     fid can be negative, take its absolute value
      cdfid = abs(fid)
       
c     Open ICONCONST      
      ierr = NF90_OPEN(constfile,nf90_nowrite, cdfid)
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
 
c     Get <ndim>, <dimids> and <dimname>
      ierr = NF90_INQ_VARID(cdfid,'z_mc',varid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_variable(cdfid, varid,ndims  = ndim)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_variable(cdfid, varid,dimids = dimids(1:ndim))
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      do i=1,ndim
         ierr = nf90_inquire_dimension(cdfid, dimids(i), 
     >                               name = dimname(i) )
         IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      enddo
      
c     Check that dimensions are OK
      if ( ndim.ne.3 ) then
        print*,' ERROR: z_mc must be 3D ',ndim
        stop
      endif
      if ( dimname(2).ne.'lat' ) then
        print*,' ERROR: dimname(2) must be lat ',trim(dimname(2))
        stop
      endif
      if ( dimname(1).ne.'lon' ) then
        print*,' ERROR: dimname(1) must be lon ',trim(dimname(1))
        stop
      endif
             
c     Get lon coordinates        
      ierr = nf90_inq_dimid(cdfid,'lon', dimid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_dimension(cdfid, dimid, len = nx)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = NF90_INQ_VARID(cdfid,'lon',varid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = NF90_GET_VAR(cdfid,varid,lon(1:nx))
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 

c     Get lat coordinates        
      ierr = nf90_inq_dimid(cdfid,'lat', dimid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_dimension(cdfid, dimid, len = ny)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = NF90_INQ_VARID(cdfid,'lat',varid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = NF90_GET_VAR(cdfid,varid,lat(1:ny))
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 

c     Get number of vertical levels
      ierr = nf90_inquire_dimension(cdfid, dimids(3), len = nz)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

c     Handling of vertical staggering 
      nz = nz -1

c     Set parameters
      pollon = 0.
      pollat = 90.
      polgam = 0.
      xmin   = lon(1)
      ymin   = lat(1)
      xmax   = lon(nx)
      ymax   = lat(ny)
      dx     = ( xmax - xmin ) / real(nx-1) 
      dy     = ( ymax - ymin ) / real(ny-1)
      if ( xmin.lt.-180.) then
        xmin = xmin + 360.
        xmax = xmax + 360.
      else if ( xmax.gt.360. ) then
        xmin = xmin - 360.
        xmax = xmax - 360.
      endif
      stagz = 0.

c     If <fid<0>, nothing further to do - we have the grid dimensions 
      if ( fid.lt.0 ) goto 100

c     Read topography 
      ierr = NF90_INQ_VARID(cdfid,'TOPOGRAPHY',varid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = NF90_GET_VAR(cdfid,varid,zb)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 

c     Read 3D height field
      ierr = NF90_INQ_VARID(cdfid,'z_mc',varid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = NF90_GET_VAR(cdfid,varid,stag3)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 

c     Handling of vertical staggering - destagger
      do i=1,nx
         do j=1,ny
                do k=1,nz
                   tmp3(i,j,k) = 0.5 * ( stag3(i,j,k) + stag3(i,j,k+1) )
                enddo
         enddo
      enddo     

c     Check whether vertical axis is descending - vertical flip
      if ( tmp3(1,1,1).lt.tmp3(1,1,nz) ) then
        print*,' ERROR: vertical axis must be descending... Stop'
        stop
      endif
      do i=1,nx
        do j=1,ny
           do k=1,nz
              z3(i,j,k) = tmp3(i,j,nz-k+1)
           enddo
        enddo
      enddo

c     Exit point
 100  continue
   
c     Close constants file
      ierr = NF90_CLOSE(cdfid)
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)   
  
      end

c     ------------------------------------------------------------
c     Read wind information
c     ------------------------------------------------------------

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

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>.

      use netcdf
      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     Auxiliary variables
      integer      i,j,k
      integer      nlon,nlat,nlev
      real         tmp2(nx,ny)
      real         tmp3(nx,ny,nz)
      real         stag3(nx,ny,nz+1)
      integer      ierr
      integer      varid
      integer      dimid
      integer      dimids(100)
      character*80 dimname(100)   
      integer      ndim

c     ----- Get dimension and check whether OK --------------------------------- 
     
c     Get <ndim>, <dimids> and <dimname>
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_variable(fid, varid,dimids = dimids(1:ndim))
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      do i=1,ndim
         ierr = nf90_inquire_dimension(fid, dimids(i), 
     >                               name = dimname(i) )
         IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      enddo

      if ( (ndim.ne.3).and.(ndim.ne.4) ) then
        print*,' Data must be either 3D or 4D'
        stop
      endif
      
c     ----- Read 4D data ------------------------------------------------------ 

      if ( ndim.ne.4 ) goto 100

c     Check that dimensions are OK
      if ( dimname(4).ne.'time' ) then
        print*,' ERROR: dimname(4) must be time '
        print*,trim(fieldname),' / ',trim(dimname(4))
        stop
      endif
      if ( dimname(2).ne.'lat' ) then
        print*,' ERROR: dimname(2) must be lat ' 
        print*,trim(fieldname),' / ',trim(dimname(2))
        stop
      endif    
      if ( dimname(1).ne.'lon' ) then
        print*,' ERROR: dimname(1) must be lon ' 
        print*,trim(fieldname),' / ',trim(dimname(2))
        stop
      endif   
      
c     Check grid dimensions       
      ierr = nf90_inq_dimid(fid,dimname(1), dimid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_dimension(fid, dimid, len = nlon)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inq_dimid(fid,dimname(2), dimid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_dimension(fid, dimid, len = nlat)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inq_dimid(fid,dimname(3), dimid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_dimension(fid, dimid, len = nlev)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      if ( (nlon.ne.nx).or.(nlat.ne.ny).or.
     >  (nlev.ne.nz).and.(nlev.ne.1).and.(nlev.ne.(nz+1)) )
     >then
         print*,' ERROR: grid mismatch between ICONCONST and P file'
         stop
      endif

c     Get varid
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

c     Read 2D variable
      if ( nlev.eq.1 ) then
        ierr = NF90_GET_VAR(fid,varid,tmp2)
        IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 
        do i=1,nx
            do j=1,ny
               do k=1,nz
                  tmp3(i,j,k) = tmp2(i,j)
               enddo
            enddo
         enddo
      endif
      
c     Read 3D variable, destagger if necessary      
      if ( nlev.eq.nz ) then
        ierr = NF90_GET_VAR(fid,varid,tmp3)
        IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 
      elseif ( nlev.eq.(nz+1) ) then  
        ierr = NF90_GET_VAR(fid,varid,stag3)
        IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 
        do i=1,nx
          do j=1,ny
             do k=1,nz
                tmp3(i,j,k) = 0.5 * ( stag3(i,j,k) + stag3(i,j,k+1) )
             enddo
          enddo
        enddo  
      endif
     
c     Flip vertically 
      do i=1,nx
        do j=1,ny
           do k=1,nz
              field(i,j,k) = tmp3(i,j,nz-k+1)
           enddo
        enddo
      enddo

c     Exit point for reading 4D data
 100  continue

c     ----- Read 3D data ------------------------------------------------------ 

      if ( ndim.ne.3 ) goto 200

c     Check that dimensions are OK
      if ( dimname(3).ne.'time' ) then
        print*,' ERROR: dimname(3) must be time '
        print*,trim(fieldname),' / ',trim(dimname(3))
        stop
      endif
      if ( dimname(2).ne.'lat' ) then
        print*,' ERROR: dimname(2) must be lat ' 
        print*,trim(fieldname),' / ',trim(dimname(2))
        stop
      endif    
      if ( dimname(1).ne.'lon' ) then
        print*,' ERROR: dimname(1) must be lon ' 
        print*,trim(fieldname),' / ',trim(dimname(2))
        stop
      endif   
      
c     Check grid dimensions       
      ierr = nf90_inq_dimid(fid,dimname(1), dimid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_dimension(fid, dimid, len = nlon)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inq_dimid(fid,dimname(2), dimid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_dimension(fid, dimid, len = nlat)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      if ( (nlon.ne.nx).or.(nlat.ne.ny) ) then
         print*,' ERROR: grid mismatch between ICONCONST and P file'
         stop
      endif

c     Get varid
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

c     Read 2D variable
      ierr = NF90_GET_VAR(fid,varid,tmp2)
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 
      do i=1,nx
          do j=1,ny
             do k=1,nz
                field(i,j,k) = tmp2(i,j)
             enddo
          enddo
      enddo

c     Exit point for reading 4D data
 200  continue
     
      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

c     Close file
      ierr = NF90_CLOSE(fid)
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 
      end