Subversion Repositories lagranto.um

Rev

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

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

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 
     >                   (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >                    time,pollon,pollat,p3,ps,nz,ak,bk,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 <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).

      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       ! 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 and physical parameters
      real          eps                 ! Numerical epsilon
      parameter    (eps=0.001)

c     Netcdf variables
      integer      vardim(4)
      real         varmin(4),varmax(4)
      real         mdv
      real         stag(4)
      integer      ndim
      character*80 cstfile
      integer      cstid
      real         times
      integer      ntimes
      real         aklay(nz),bklay(nz),aklev(nz),bklev(nz)
      integer      nvars
      character*80 vars(100)
      integer        dimids (nf90_max_var_dims)
      character*80   dimname(nf90_max_var_dims)

c     Auxiliary varaibles
      integer      ierr       
      integer      i,j,k
      integer      isok
      real         tmp(200)
      character*80 varname
      real         rtime
      integer      varid
      integer      cdfid
      real         rmax,rmin

c     Set file identifier
      if (fid.lt.0) then
        cdfid = -fid
      else 
        cdfid = fid
      endif
      
c     Get varid
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
         
c     Get numebr of dimensions -> ndim      
      ierr = nf90_inquire_variable(cdfid, varid, ndims  = ndim)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
         
c     Get dim IDs         
      ierr = nf90_inquire_variable(cdfid, varid, 
     >                                   dimids = dimids(1:ndim))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

c     Get dimensions -> vardim(1:ndim)         
      do i=1,ndim
           ierr = nf90_inquire_dimension(cdfid, dimids(i), 
     >                               name = dimname(i) )
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
           ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      enddo
         
c     Get domain min -> varmin(1:3)        
      ierr  = nf90_get_att(cdfid, varid, "xmin", varmin(1))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr  = nf90_get_att(cdfid, varid, "ymin", varmin(2))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr  = nf90_get_att(cdfid, varid, "zmin", varmin(3))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
c     Get domain max -> varmax(1:3)     
      ierr  = nf90_get_att(cdfid, varid, "xmax", varmax(1))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr  = nf90_get_att(cdfid, varid, "ymax", varmax(2))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr  = nf90_get_att(cdfid, varid, "zmax", varmax(3))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
c     Get vertical staggering -> stagz           
      ierr  = nf90_get_att(cdfid, varid, "zstag", stagz)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
c     Get missing data value -> mdv         
      ierr  = nf90_get_att(cdfid, varid, "missing_value", mdv)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

c     Set the grid dimensions and constants
      nx   = vardim(1)
      ny   = vardim(2)
      nz   = vardim(3)
      xmin = varmin(1)
      ymin = varmin(2)
      xmax = varmax(1)
      ymax = varmax(2)
      dx   = (xmax-xmin)/real(nx-1)
      dy   = (ymax-ymin)/real(ny-1)

c     We want the longitudes within -180 ... + 180
      if ( xmin.lt.-180.) then
            xmin = xmin + 360.
            xmax = xmax + 360.
      else if ( xmax.gt.360. ) then
            xmin = xmin - 360.
            xmax = xmax - 360.
      endif

c     Get name of constants file -> cstfile
      ierr  = nf90_get_att(cdfid,nf90_global,
     >                    "constants_file_name ", cstfile )
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

c     Get pole position from constants file
      ierr = NF90_OPEN(TRIM(cstfile),nf90_nowrite, cstid)
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
      varname='pollon'
      ierr = NF90_INQ_VARID(cstid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cstid,varid,pollon)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname='pollat'
      ierr = NF90_INQ_VARID(cstid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cstid,varid,pollat)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = NF90_CLOSE(cstid)
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
c     Skip the rest if fid < 0: 
      if ( fid.lt.0 ) goto 900      
      
c     Get ak and bk from constants file
      ierr = NF90_OPEN(TRIM(cstfile),nf90_nowrite, cstid)
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
      varname='aklev'
      ierr = NF90_INQ_VARID(cstid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cstid,varid,aklev)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname='bklev'
      ierr = NF90_INQ_VARID(cstid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cstid,varid,bklev)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname='aklay'
      ierr = NF90_INQ_VARID(cstid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cstid,varid,aklay)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname='bklay'
      ierr = NF90_INQ_VARID(cstid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cstid,varid,bklay)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = NF90_CLOSE(cstid)
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
c     print*,'AKLEV : '
c     print*, (aklev(i),i=1,20)
c     print*,'BKLEV : '
c     print*,(bklev(i),i=1,20)
c     print*,'AKLAY : '
c     print*,(aklay(i),i=1,20)
c     print*,'BKLAY : '
c     print*,(bklay(i),i=1,20)
      
c     Get time information (check if time is correct)
      varname = 'time'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,times)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      isok=0
      if (abs(time-times).lt.eps) then
               isok = 1
               rtime = times
      elseif (timecheck.eq.'no') then
               isok = 1
               rtime = times
      endif
      if ( isok.eq.0 ) then
         print*,' ERROR: time ',rtime,' not found on netCDF file' 
         stop
      endif
      
c     print*,'TIME : ',rtime

c     Read surface pressure
      varname='PS'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,ps)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

c     Check that PS is in hPa - quick and dirty
      rmax = -1.e9
      rmin = +1.e9
      do i=1,nx
       do j=1,ny
        if ( (ps(i,j).gt.rmax).and.abs(ps(i,j)-mdv).gt.eps)  then
          rmax = ps(i,j)
        endif
        if ( (ps(i,j).lt.rmin).and.abs(ps(i,j)-mdv).gt.eps)  then
          rmin = ps(i,j)
        endif
       enddo
      enddo
      if ( (rmin.lt.200.).or.(rmax.gt.1500.) ) then
        print*,' ERROR: PS must be in hPa... Stop'
        stop
      endif
    
c     Calculate layer and level pressures
      do i=1,nx
         do j=1,ny
               do k=1,nz
                  if ( abs(stagz).lt.eps ) then
                     p3(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
                  else
                     p3(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
                  endif
               enddo
         enddo
      enddo

c     Set the ak and bk for the vertical grid
      do k=1,nz
            if ( abs(stagz).lt.eps ) then
               ak(k)=aklev(k)
               bk(k)=bklev(k)
            else
               ak(k)=aklay(k)
               bk(k)=bklay(k)
            endif
      enddo
      
c     Exit point
 900  continue    
   
      return
      
      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     Numerical and physical parameters
      real        eps                 ! Numerical epsilon
      parameter  (eps=0.001)
      real        notimecheck         ! 'Flag' for no time check
      parameter  (notimecheck=7.26537)

c     Netcdf variables
      integer      ierr
      character*80 varname
      integer      vardim(4)
      real         varmin(4),varmax(4)
      real         stag(4)
      integer      ndim
      real         times(10)
      integer      ntimes
      character*80 cstfile
      integer      cstid
      real         aklay(200),bklay(200),aklev(200),bklev(200)
      real         ps(nx,ny)
      integer      dimids (nf90_max_var_dims)
      character*80 dimname(nf90_max_var_dims)
      integer      varid
      integer      cdfid

c     Auxiliary variables
      integer      isok
      integer      i,j,k
      integer      nz1
      real         rtime

c     Get varid
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
         
c     Get number of vertical levels -> vardim(3)    
      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)
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      enddo

c     Read data 
      varname=fieldname
      ierr = nf90_get_var(fid,varid,field)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
  
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
      if ( vardim(3).eq.1 ) then
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  field(i,j,k) = field(i,j,1)
               enddo
            enddo
         enddo
      endif
         
c     Exit point
      return
 
      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