WebSVN – lagranto.ecmwf – Diff – /trunk/lib/ioinp.f


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
 
      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
      integer      nvars
      character*80 vars(100)
      integer        dimids (nf90_max_var_dims),dimid
      character*80   dimname(nf90_max_var_dims)
      character*80   stdname
      real,allocatable, dimension (:)     :: lon,lat,lev
      real,allocatable, dimension (:)     :: times
      real,allocatable, dimension (:,:)   :: tmp2
      real,allocatable, dimension (:,:,:) :: tmp3
      real,allocatable, dimension (:)     :: aktmp,bktmp
      character*80  units
      character*80  leveltype
      integer       nakbktmp
      integer       vertical_swap
 
c     Auxiliary variables
      integer      ierr       
      integer      i,j,k
      integer      isok
      real         tmp(200)
      character*80 varname
      real         rtime
      integer      varid
      integer      cdfid
      integer      stat
      real         delta
      integer      closear
      real         maxps,minps
 
c     ---- Read data from netCDF file as they are ---------------------
 
c     Set file identifier
      if (fid.lt.0) then
        cdfid = -fid
      else 
        cdfid = fid
      endif
 
c     Special handling if 3D pressure is
      if ( fieldname.eq.'P' ) then
         fieldname = 'U'
      endif
 
c     Get number of dimensions of variable -> ndim
      ierr = NF90_INQ_VARID(cdfid,fieldname,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)
      if ( ndim.ne.4 ) then
          print*,' ERROR: netCDF variables need to be 4D'
          print*,'      ',trim(fieldname)
          stop
      endif
 
c     Get dimensions -> vardim(1:ndim),dimname(1:ndim)
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
      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)
           ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      enddo
 
c     Get dimension of AK,BK
      varname = 'nhym'
      ierr = NF90_INQ_DIMID(cdfid,varname,dimid)
      ierr = nf90_inquire_dimension(cdfid, dimid,len=nakbktmp)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 
c     Check whether the list of dimensions is OK
      if ( ( dimname(1).ne.'lon'  ).or.
     >     ( dimname(2).ne.'lat'  ).or. 
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
     >     ( dimname(4).ne.'time' ) )
     >then
        print*,' ERROR: the dimensions of the variable are not correct'
        print*,'        expected -> lon / lat / lev / time'
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
        stop
      endif
 
c     Check about the type of vertical levels
      varname=dimname(3)
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_att(cdfid, varid, "standard_name", leveltype)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      if ( (leveltype.ne.'hybrid_sigma_pressure').and.
     >     (leveltype.ne.'air_pressure'         ) )
     >then
         print*,' ERROR: input netCDF data must be on hybrid-sigma'
         print*,'        or air pressure levels!',trim(leveltype)
         stop
      endif
 
c     Allocate memory for reading arrays
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
      allocate(lon(vardim(1)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lon     ***' 
      allocate(lat(vardim(2)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lat     ***' 
      allocate(lev(vardim(3)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
      allocate(times(vardim(4)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array times   ***'
      allocate(aktmp(nakbktmp),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
      allocate(bktmp(nakbktmp),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
 
c     Get domain longitudes, latitudes and levels
      varname = dimname(1)
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,lon)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname = dimname(2)
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,lat)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname = dimname(3)
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,lev)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
c     Get ak and bk
      varname='hyam'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,aktmp)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname='hybm'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,bktmp)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 
c     Check that unit of ak is in hPa - if necessary correct it
      varname='hyam'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_att(cdfid, varid, "units", units)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      if ( units.eq.'Pa' ) then
         do k=1,nakbktmp
            aktmp(k) = 0.01 * aktmp(k)
         enddo
      endif
 
c     Check that unit of lev is in hPa - if necessary correct it
      if ( leveltype.eq.'air_pressure' ) then
         varname='lev'
         ierr = NF90_INQ_VARID(cdfid,varname,varid)
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
         ierr = nf90_get_att(cdfid, varid, "units", units)
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
         if ( units.eq.'Pa' ) then
            do k=1,vardim(3)
               lev(k) = 0.01 * lev(k)
            enddo
         endif
      endif
 
c     Decide whether to swap vertical levels - highest pressure at index 1
      vertical_swap = 1
      if ( leveltype.eq.'hybrid_sigma_pressure') then
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
     >       (aktmp(2) + bktmp(2) * 1000.) )
     >  then
          vertical_swap = 0
        endif
      elseif ( leveltype.eq.'air_pressure') then
        if ( lev(1).gt.lev(2) ) then
          vertical_swap = 0
        endif
      endif
c      print*,' Vertical SWAP P -> ', vertical_swap
 
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
      do i=1,vardim(4)
        if (abs(time-times(i)).lt.eps) then
               isok = 1
               rtime = times(i)
        elseif (timecheck.eq.'no') then
               isok = 1
               rtime = times(1)
        endif
      enddo
      if ( isok.eq.0 ) then
         print*,' ERROR: time ',rtime,' not found on netCDF file' 
         stop
      endif
 
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,tmp2)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
    
c     Check that surface pressure is in hPa
      maxps = -1.e19
      minps =  1.e19
      do i=1,vardim(1)
        do j=1,vardim(2)
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
             if (tmp2(i,j).lt.minps) minps = tmp2(i,j)
        enddo
      enddo
      if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then
         print*,' ERROR: surface pressre PS must be in hPa'
         print*,'       ',maxps,minps
         stop
      endif
 
c     ---- Define output of subroutine --------------------------------
 
c     If not full list of vertical levels, reduce AK,BK arrays
      if ( (leveltype.eq.'hybrid_sigma_pressure').and.
     >     (nakbktmp.ne.vardim(3) ) )
     >then
         print*,' WARNING: only subset of vertical levels used...'
         do k=1,vardim(3)
            if ( vertical_swap.eq.1 ) then
               aktmp(k) = aktmp( k+nakbktmp-vardim(3) )
               bktmp(k) = bktmp( k+nakbktmp-vardim(3) )
            endif
         enddo
      endif
 
c     Set the grid dimensions and constants
      nx      = vardim(1)
      ny      = vardim(2)
      nz      = vardim(3)
      xmin    = lon(1)
      ymin    = lat(1)
      xmax    = lon(nx)
      ymax    = lat(ny)
      dx      = (xmax-xmin)/real(nx-1)
      dy      = (ymax-ymin)/real(ny-1)
      pollon  = 0.
      pollat  = 90.
      stagz   = 0.
      delta   = xmax-xmin-360.
      if (abs(delta+dx).lt.eps) then
          xmax    = xmax + dx
          nx      = nx + 1
          closear = 1
      else
          closear = 0
      endif
 
c     Save the output arrays (if fid>0) - close arrays on request
      if ( fid.gt.0 ) then
 
c        Calculate layer pressures
         if (leveltype.eq.'hybrid_sigma_pressure' ) then
            do i=1,vardim(1)
              do j=1,vardim(2)
                 do k=1,vardim(3)
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
                 enddo
              enddo
           enddo
         elseif (leveltype.eq.'air_pressure' ) then
           do i=1,vardim(1)
              do j=1,vardim(2)
                 do k=1,vardim(3)
                  tmp3(i,j,k)=lev(k)
                 enddo
              enddo
           enddo
         endif
 
c        Get PS - close array on demand
         do j=1,vardim(2)
           do i=1,vardim(1)
             ps(i,j) = tmp2(i,j)
           enddo
           if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)
         enddo
 
c        Get P3 - close array on demand + vertical swap
         do j=1,vardim(2)
           do k=1,vardim(3)
             do i=1,vardim(1)
               if ( vertical_swap.eq.1 ) then
                  p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)
               else
                  p3(i,j,k) = tmp3(i,j,k)
               endif
             enddo
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
           enddo
         enddo
 
c        Get AK,BK - vertical swap on demand
         if ( leveltype.eq.'hybrid_sigma_pressure' ) then
           do k=1,vardim(3)
              if ( vertical_swap.eq.1 ) then
                 ak(k) = aktmp(vardim(3)-k+1)
                 bk(k) = bktmp(vardim(3)-k+1)
              endif
           enddo
         elseif (leveltype.eq.'air_pressure' ) then
           do k=1,vardim(3)
              if ( vertical_swap.eq.1 ) then
                 ak(k) = lev(vardim(3)-k+1)
                 bk(k) = 0.
              else
                ak(k) = lev(k)
                bk(k) = 0.
              endif
           enddo
         endif
 
      endif
 
 
      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
      real,allocatable, dimension (:)     :: lon,lat,lev
      real,allocatable, dimension (:,:)   :: tmp2
      real,allocatable, dimension (:,:,:) :: tmp3
      real,allocatable, dimension (:)     :: aktmp,bktmp
      character*80  leveltype
      integer       vertical_swap
      character*80  units
      integer       nakbktmp
      integer       dimid
 
c     Auxiliary variables
      integer      isok
      integer      i,j,k
      integer      nz1
      real         rtime
      integer      closear
      integer      stat
      real         delta
 
c     Init mdv
      mdv = -999.
 
c     Get the number of dimensions -> ndim
      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)
 
c     Get the dimensions of the arrays -> varid(1:ndim)
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
      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     Check whether the list of dimensions is OK
      if ( ( dimname(1).ne.'lon'  ).or.
     >     ( dimname(2).ne.'lat'  ).or.
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
     >     ( dimname(4).ne.'time' ) )
     >then
        print*,' ERROR: the dimensions of the variable are not correct'
        print*,'        expected -> lon / lat / lev / time'
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
        stop
      endif
 
c     Get dimension of AK,BK
      varname = 'nhym'
      ierr = NF90_INQ_DIMID(fid,varname,dimid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 
c     Check about the type of vertical levels
      varname=dimname(3)
      ierr = NF90_INQ_VARID(fid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_att(fid, varid, "standard_name", leveltype)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      if ( (leveltype.ne.'hybrid_sigma_pressure').and.
     >     (leveltype.ne.'air_pressure'         ) )
     >then
         print*,' ERROR: input netCDF data must be on hybrid-sigma'
         print*,'        or air pressure levels!',trim(leveltype)
         stop
      endif
 
c     Allocate memory for reading arrays - depending on <closear>
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
      allocate(lon(vardim(1)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
      allocate(lat(vardim(2)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
      allocate(lev(vardim(3)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
      allocate(aktmp(nakbktmp),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
      allocate(bktmp(nakbktmp),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
 
c     Get domain boundaries - longitude, latitude, levels
      varname = dimname(1)
      ierr = NF90_INQ_VARID(fid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(fid,varid,lon)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname = dimname(2)
      ierr = NF90_INQ_VARID(fid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(fid,varid,lat)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname = dimname(3)
      ierr = NF90_INQ_VARID(fid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(fid,varid,lev)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 
c     Get ak and bk
      varname='hyam'
      ierr = NF90_INQ_VARID(fid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(fid,varid,aktmp)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname='hybm'
      ierr = NF90_INQ_VARID(fid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(fid,varid,bktmp)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 
c     Check that unit of ak is in hPa - if necessary correct it
      varname='hyam'
      ierr = NF90_INQ_VARID(fid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_att(fid, varid, "units", units)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      if ( units.eq.'Pa' ) then
         do k=1,nakbktmp
            aktmp(k) = 0.01 * aktmp(k)
         enddo
      endif
 
c     Check that unit of lev is in hPa - if necessary correct it
      if ( leveltype.eq.'air_pressure' ) then
         varname='lev'
         ierr = NF90_INQ_VARID(fid,varname,varid)
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
         ierr = nf90_get_att(fid, varid, "units", units)
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
         if ( units.eq.'Pa' ) then
            do k=1,vardim(3)
               lev(k) = 0.01 * lev(k)
            enddo
         endif
      endif
 
c     Decide whether to swap vertical levels
      vertical_swap = 1
      if ( leveltype.eq.'hybrid_sigma_pressure') then
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
     >       (aktmp(2) + bktmp(2) * 1000.) )
     >  then
          vertical_swap = 0
        endif
      elseif ( leveltype.eq.'air_pressure') then
        if ( lev(1).gt.lev(2) ) then
          vertical_swap = 0
        endif
      endif
 
c     Read data 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(fid,varid,tmp3)
      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,vardim(1)
            do j=1,vardim(2)
               do k=1,vardim(3)
                  tmp3(i,j,k) = tmp3(i,j,1)
               enddo
            enddo
         enddo
      endif
 
c     Decide whether to close arrays
      delta = varmax(1)-varmin(1)-360.
      if (abs(delta+dx).lt.eps) then
          closear = 1
      else
          closear = 0
      endif
 
c     Save output array - close array and swap on demand
      do j=1,vardim(2)
        do k=1,vardim(3)
          do i=1,vardim(1)
             if ( vertical_swap.eq.1 ) then
                 field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
             else
                 field(i,j,k) = tmp3(i,j,k)
             endif
          enddo
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
        enddo
      enddo
         
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
      
c     ------------------------------------------------------------
c     Get a list of variables on netCDF file
c     ------------------------------------------------------------
 
      subroutine input_getvars(fid,vnam,nvars)
 
c     List of variables on netCDF file
 
      use netcdf
 
      implicit none
 
c     Declaration of subroutine parameters
      integer      fid
      integer      nvars
      character*80 vnam(200)
 
c     Auxiliary variables
      integer ierr
      integer i
      integer nDims,nGlobalAtts,unlimdimid
 
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
      do i=1,nVars
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
      enddo
 
      end

Rev 27	Rev 29
1	`c ************************************************************`	1	`c ************************************************************`
2	`c * This package provides input routines to read the wind *`	2	`c * This package provides input routines to read the wind *`
3	`c * and other fields from IVE necdf files. The routines are *`	3	`c * and other fields from IVE necdf files. The routines are *`
4	`c * *`	4	`c * *`
5	`c * 1) input_open : to open a data file *`	5	`c * 1) input_open : to open a data file *`
6	`c * 2) input_grid : to read the grid information, including *`	6	`c * 2) input_grid : to read the grid information, including *`
7	`c * the vertical levels *`	7	`c * the vertical levels *`
8	`c * 3) input_wind : to read the wind components *`	8	`c * 3) input_wind : to read the wind components *`
9	`c * 4) input_close : to close an input file *`	9	`c * 4) input_close : to close an input file *`
10	`c * *`	10	`c * *`
11	`c * The file is characterised by an filename <filename> and *`	11	`c * The file is characterised by an filename <filename> and *`
12	`c * a file identifier <fid>. The horizontal grid is given by *`	12	`c * a file identifier <fid>. The horizontal grid is given by *`
13	`c * <xmin,xmax,ymin,ymax,dx,dy,nx,ny> where the pole of the *`	13	`c * <xmin,xmax,ymin,ymax,dx,dy,nx,ny> where the pole of the *`
14	`c * rotated grid is given by <pollon,pollat>. The vertical *`	14	`c * rotated grid is given by <pollon,pollat>. The vertical *`
15	`c * grid is characterised by the surface pressure <ps> and *`	15	`c * grid is characterised by the surface pressure <ps> and *`
16	`c * the pressure at staggered <slev> and unstaggered <ulev> *`	16	`c * the pressure at staggered <slev> and unstaggered <ulev> *`
17	`c * levels. The number of levels is given by <nz>. Finally, *`	17	`c * levels. The number of levels is given by <nz>. Finally, *`
18	`c * the retrieval of the wind <field> with name <fieldname> *`	18	`c * the retrieval of the wind <field> with name <fieldname> *`
19	`c * is characterised by a <time> and a missing data value *`	19	`c * is characterised by a <time> and a missing data value *`
20	`c * <mdv>. *`	20	`c * <mdv>. *`
21	`c * *`	21	`c * *`
22	`c * Author: Michael Sprenger, Autumn 2008 *`	22	`c * Author: Michael Sprenger, Autumn 2008 *`
23	`c ************************************************************`	23	`c ************************************************************`
24		24
25	`c ------------------------------------------------------------`	25	`c ------------------------------------------------------------`
26	`c Open input file`	26	`c Open input file`
27	`c ------------------------------------------------------------`	27	`c ------------------------------------------------------------`
28		28
29	`subroutine input_open (fid,filename)`	29	`subroutine input_open (fid,filename)`
30		30
31	`c Open the input file with filename <filename> and return the`	31	`c Open the input file with filename <filename> and return the`
32	`c file identifier <fid> for further reference.`	32	`c file identifier <fid> for further reference.`
33		33
34	`use netcdf`	34	`use netcdf`
35		35
36	`implicit none`	36	`implicit none`
37		37
38	`c Declaration of subroutine parameters`	38	`c Declaration of subroutine parameters`
39	`integer fid ! File identifier`	39	`integer fid ! File identifier`
40	`character*80 filename ! Filename`	40	`character*80 filename ! Filename`
41		41
42	`c Declaration of auxiliary variables`	42	`c Declaration of auxiliary variables`
43	`integer ierr`	43	`integer ierr`
44		44
45	`c Open netcdf file`	45	`c Open netcdf file`
46	`ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)`	46	`ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)`
47	`IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)`	47	`IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)`
48		48
49	`c Exception handling`	49	`c Exception handling`
50	`return`	50	`return`
51		51
52	`end`	52	`end`
53		53
54	`c ------------------------------------------------------------`	54	`c ------------------------------------------------------------`
55	`c Read information about the grid`	55	`c Read information about the grid`
56	`c ------------------------------------------------------------`	56	`c ------------------------------------------------------------`
57		57
58	`subroutine input_grid`	58	`subroutine input_grid`
59	`> (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`	59	`> (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`
60	`> time,pollon,pollat,p3,ps,nz,ak,bk,stagz,`	60	`> time,pollon,pollat,p3,ps,nz,ak,bk,stagz,`
61	`> timecheck)`	61	`> timecheck)`
62		62
63	`c Read grid information at <time> from file with identifier <fid>.`	63	`c Read grid information at <time> from file with identifier <fid>.`
64	`c The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>`	64	`c The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>`
65	`c with pole position at <pollon,pollat> and grid dimension <nx,ny>.`	65	`c with pole position at <pollon,pollat> and grid dimension <nx,ny>.`
66	`c The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either`	66	`c The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either`
67	`c on the staggered or unstaggered grid (with <stagz> as the flag).`	67	`c on the staggered or unstaggered grid (with <stagz> as the flag).`
68	`c The surface pressure is given in <ps(nx,ny)>. If <fid> is negative,`	68	`c The surface pressure is given in <ps(nx,ny)>. If <fid> is negative,`
69	`c only the grid dimensions and grid parameters (xmin...pollat,nz) are`	69	`c only the grid dimensions and grid parameters (xmin...pollat,nz) are`
70	`c determined and returned (this is needed for dynamical allocation of`	70	`c determined and returned (this is needed for dynamical allocation of`
71	`c memory).`	71	`c memory).`
72		72
73	`use netcdf`	73	`use netcdf`
74		74
75	`implicit none`	75	`implicit none`
76		76
77	`c Declaration of subroutine parameters`	77	`c Declaration of subroutine parameters`
78	`integer fid ! File identifier`	78	`integer fid ! File identifier`
79	`real xmin,xmax,ymin,ymax ! Domain size`	79	`real xmin,xmax,ymin,ymax ! Domain size`
80	`real dx,dy ! Horizontal resolution`	80	`real dx,dy ! Horizontal resolution`
81	`integer nx,ny,nz ! Grid dimensions`	81	`integer nx,ny,nz ! Grid dimensions`
82	`real pollon,pollat ! Longitude and latitude of pole`	82	`real pollon,pollat ! Longitude and latitude of pole`
83	`real p3(nx,ny,nz) ! Staggered levels`	83	`real p3(nx,ny,nz) ! Staggered levels`
84	`real ps(nx,ny) ! Surface pressure`	84	`real ps(nx,ny) ! Surface pressure`
85	`real time ! Time of the grid information`	85	`real time ! Time of the grid information`
86	`real ak(nz),bk(nz) ! Ak and Bk for layers or levels`	86	`real ak(nz),bk(nz) ! Ak and Bk for layers or levels`
87	`real stagz ! Vertical staggering (0 or -0.5)`	87	`real stagz ! Vertical staggering (0 or -0.5)`
88	`character*80 fieldname ! Variable from which to take grid info`	88	`character*80 fieldname ! Variable from which to take grid info`
89	`character*80 timecheck ! Either 'yes' or 'no'`	89	`character*80 timecheck ! Either 'yes' or 'no'`
90		90
91	`c Numerical and physical parameters`	91	`c Numerical and physical parameters`
92	`real eps ! Numerical epsilon`	92	`real eps ! Numerical epsilon`
93	`parameter (eps=0.001)`	93	`parameter (eps=0.001)`
94		94
95	`c Netcdf variables`	95	`c Netcdf variables`
96	`integer vardim(4)`	96	`integer vardim(4)`
97	`real varmin(4),varmax(4)`	97	`real varmin(4),varmax(4)`
98	`real mdv`	98	`real mdv`
99	`real stag(4)`	99	`real stag(4)`
100	`integer ndim`	100	`integer ndim`
101	`character*80 cstfile`	101	`character*80 cstfile`
102	`integer cstid`	102	`integer cstid`
103	`integer nvars`	103	`integer nvars`
104	`character*80 vars(100)`	104	`character*80 vars(100)`
105	`integer dimids (nf90_max_var_dims),dimid`	105	`integer dimids (nf90_max_var_dims),dimid`
106	`character*80 dimname(nf90_max_var_dims)`	106	`character*80 dimname(nf90_max_var_dims)`
107	`character*80 stdname`	107	`character*80 stdname`
108	`real,allocatable, dimension (:) :: lon,lat,lev`	108	`real,allocatable, dimension (:) :: lon,lat,lev`
109	`real,allocatable, dimension (:) :: times`	109	`real,allocatable, dimension (:) :: times`
110	`real,allocatable, dimension (:,:) :: tmp2`	110	`real,allocatable, dimension (:,:) :: tmp2`
111	`real,allocatable, dimension (:,:,:) :: tmp3`	111	`real,allocatable, dimension (:,:,:) :: tmp3`
112	`real,allocatable, dimension (:) :: aktmp,bktmp`	112	`real,allocatable, dimension (:) :: aktmp,bktmp`
113	`character*80 units`	113	`character*80 units`
114	`character*80 leveltype`	114	`character*80 leveltype`
115	`integer nakbktmp`	115	`integer nakbktmp`
116	`integer vertical_swap`	116	`integer vertical_swap`
117		117
118	`c Auxiliary variables`	118	`c Auxiliary variables`
119	`integer ierr`	119	`integer ierr`
120	`integer i,j,k`	120	`integer i,j,k`
121	`integer isok`	121	`integer isok`
122	`real tmp(200)`	122	`real tmp(200)`
123	`character*80 varname`	123	`character*80 varname`
124	`real rtime`	124	`real rtime`
125	`integer varid`	125	`integer varid`
126	`integer cdfid`	126	`integer cdfid`
127	`integer stat`	127	`integer stat`
128	`real delta`	128	`real delta`
129	`integer closear`	129	`integer closear`
130	`real maxps,minps`	130	`real maxps,minps`
131		131
132	`c ---- Read data from netCDF file as they are ---------------------`	132	`c ---- Read data from netCDF file as they are ---------------------`
133		133
134	`c Set file identifier`	134	`c Set file identifier`
135	`if (fid.lt.0) then`	135	`if (fid.lt.0) then`
136	`cdfid = -fid`	136	`cdfid = -fid`
137	`else`	137	`else`
138	`cdfid = fid`	138	`cdfid = fid`
139	`endif`	139	`endif`
140		140
141	`c Special handling if 3D pressure is`	141	`c Special handling if 3D pressure is`
142	`if ( fieldname.eq.'P' ) then`	142	`if ( fieldname.eq.'P' ) then`
143	`fieldname = 'U'`	143	`fieldname = 'U'`
144	`endif`	144	`endif`
145		145
146	`c Get number of dimensions of variable -> ndim`	146	`c Get number of dimensions of variable -> ndim`
147	`ierr = NF90_INQ_VARID(cdfid,fieldname,varid)`	147	`ierr = NF90_INQ_VARID(cdfid,fieldname,varid)`
148	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	148	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
149	`ierr = nf90_inquire_variable(cdfid, varid, ndims = ndim)`	149	`ierr = nf90_inquire_variable(cdfid, varid, ndims = ndim)`
150	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	150	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
151	`if ( ndim.ne.4 ) then`	151	`if ( ndim.ne.4 ) then`
152	`print*,' ERROR: netCDF variables need to be 4D'`	152	`print*,' ERROR: netCDF variables need to be 4D'`
153	`print*,' ',trim(fieldname)`	153	`print*,' ',trim(fieldname)`
154	`stop`	154	`stop`
155	`endif`	155	`endif`
156		156
157	`c Get dimensions -> vardim(1:ndim),dimname(1:ndim)`	157	`c Get dimensions -> vardim(1:ndim),dimname(1:ndim)`
158	`ierr = NF90_INQ_VARID(cdfid,fieldname,varid)`	158	`ierr = NF90_INQ_VARID(cdfid,fieldname,varid)`
159	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	159	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
160	`ierr = nf90_inquire_variable(cdfid, varid,`	160	`ierr = nf90_inquire_variable(cdfid, varid,`
161	`> dimids = dimids(1:ndim))`	161	`> dimids = dimids(1:ndim))`
162	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	162	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
163	`do i=1,ndim`	163	`do i=1,ndim`
164	`ierr = nf90_inquire_dimension(cdfid, dimids(i),`	164	`ierr = nf90_inquire_dimension(cdfid, dimids(i),`
165	`> name = dimname(i) )`	165	`> name = dimname(i) )`
166	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	166	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
167	`ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))`	167	`ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))`
168	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	168	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
169	`enddo`	169	`enddo`
170		170
171	`c Get dimension of AK,BK`	171	`c Get dimension of AK,BK`
172	`varname = 'nhym'`	172	`varname = 'nhym'`
173	`ierr = NF90_INQ_DIMID(cdfid,varname,dimid)`	173	`ierr = NF90_INQ_DIMID(cdfid,varname,dimid)`
174	`ierr = nf90_inquire_dimension(cdfid, dimid,len=nakbktmp)`	174	`ierr = nf90_inquire_dimension(cdfid, dimid,len=nakbktmp)`
175	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	175	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
176		176
177	`c Check whether the list of dimensions is OK`	177	`c Check whether the list of dimensions is OK`
178	`if ( ( dimname(1).ne.'lon' ).or.`	178	`if ( ( dimname(1).ne.'lon' ).or.`
179	`> ( dimname(2).ne.'lat' ).or.`	179	`> ( dimname(2).ne.'lat' ).or.`
180	`> ( dimname(3).ne.'lev' ).and.( dimname(3).ne.'lev_2' ).or.`	180	`> ( dimname(3).ne.'lev' ).and.( dimname(3).ne.'lev_2' ).or.`
181	`> ( dimname(4).ne.'time' ) )`	181	`> ( dimname(4).ne.'time' ) )`
182	`>then`	182	`>then`
183	`print*,' ERROR: the dimensions of the variable are not correct'`	183	`print*,' ERROR: the dimensions of the variable are not correct'`
184	`print*,' expected -> lon / lat / lev / time'`	184	`print*,' expected -> lon / lat / lev / time'`
185	`print*, ( trim(dimname(i))//' / ',i=1,ndim )`	185	`print*, ( trim(dimname(i))//' / ',i=1,ndim )`
186	`stop`	186	`stop`
187	`endif`	187	`endif`
188		188
-		189	`c Check about the type of vertical levels`
-		190	`varname=dimname(3)`
-		191	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
-		192	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
-		193	`ierr = nf90_get_att(cdfid, varid, "standard_name", leveltype)`
-		194	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
-		195	`if ( (leveltype.ne.'hybrid_sigma_pressure').and.`
-		196	`> (leveltype.ne.'air_pressure' ) )`
-		197	`>then`
-		198	`print*,' ERROR: input netCDF data must be on hybrid-sigma'`
-		199	`print*,' or air pressure levels!',trim(leveltype)`
-		200	`stop`
-		201	`endif`
-		202
189	`c Allocate memory for reading arrays`	203	`c Allocate memory for reading arrays`
190	`allocate(tmp2(vardim(1),vardim(2)),stat=stat)`	204	`allocate(tmp2(vardim(1),vardim(2)),stat=stat)`
191	`if (stat.ne.0) print,' error allocating array tmp2 *'`	205	`if (stat.ne.0) print,' error allocating array tmp2 *'`
192	`allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)`	206	`allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)`
193	`if (stat.ne.0) print,' error allocating array tmp3 *'`	207	`if (stat.ne.0) print,' error allocating array tmp3 *'`
194	`allocate(lon(vardim(1)),stat=stat)`	208	`allocate(lon(vardim(1)),stat=stat)`
195	`if (stat.ne.0) print,' error allocating array lon *'`	209	`if (stat.ne.0) print,' error allocating array lon *'`
196	`allocate(lat(vardim(2)),stat=stat)`	210	`allocate(lat(vardim(2)),stat=stat)`
197	`if (stat.ne.0) print,' error allocating array lat *'`	211	`if (stat.ne.0) print,' error allocating array lat *'`
198	`allocate(lev(vardim(3)),stat=stat)`	212	`allocate(lev(vardim(3)),stat=stat)`
199	`if (stat.ne.0) print,' error allocating array lev *'`	213	`if (stat.ne.0) print,' error allocating array lev *'`
200	`allocate(times(vardim(4)),stat=stat)`	214	`allocate(times(vardim(4)),stat=stat)`
201	`if (stat.ne.0) print,' error allocating array times *'`	215	`if (stat.ne.0) print,' error allocating array times *'`
202	`allocate(aktmp(nakbktmp),stat=stat)`	216	`allocate(aktmp(nakbktmp),stat=stat)`
203	`if (stat.ne.0) print,' error allocating array aktmp *'`	217	`if (stat.ne.0) print,' error allocating array aktmp *'`
204	`allocate(bktmp(nakbktmp),stat=stat)`	218	`allocate(bktmp(nakbktmp),stat=stat)`
205	`if (stat.ne.0) print,' error allocating array bktmp *'`	219	`if (stat.ne.0) print,' error allocating array bktmp *'`
206		220
207	`c Get domain longitudes, latitudes and levels`	221	`c Get domain longitudes, latitudes and levels`
208	`varname = dimname(1)`	222	`varname = dimname(1)`
209	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	223	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
210	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	224	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
211	`ierr = nf90_get_var(cdfid,varid,lon)`	225	`ierr = nf90_get_var(cdfid,varid,lon)`
212	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	226	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
213	`varname = dimname(2)`	227	`varname = dimname(2)`
214	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	228	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
215	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	229	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
216	`ierr = nf90_get_var(cdfid,varid,lat)`	230	`ierr = nf90_get_var(cdfid,varid,lat)`
217	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	231	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
218	`varname = dimname(3)`	232	`varname = dimname(3)`
219	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	233	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
220	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	234	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
221	`ierr = nf90_get_var(cdfid,varid,lev)`	235	`ierr = nf90_get_var(cdfid,varid,lev)`
222	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	236	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
223		237
224	`c Get ak and bk`	238	`c Get ak and bk`
225	`varname='hyam'`	239	`varname='hyam'`
226	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	240	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
227	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	241	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
228	`ierr = nf90_get_var(cdfid,varid,aktmp)`	242	`ierr = nf90_get_var(cdfid,varid,aktmp)`
229	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	243	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
230	`varname='hybm'`	244	`varname='hybm'`
231	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	245	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
232	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	246	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
233	`ierr = nf90_get_var(cdfid,varid,bktmp)`	247	`ierr = nf90_get_var(cdfid,varid,bktmp)`
234	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	248	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
235		249
236	`c Check that unit of ak is in hPa - if necessary correct it`	250	`c Check that unit of ak is in hPa - if necessary correct it`
237	`varname='hyam'`	251	`varname='hyam'`
238	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	252	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
239	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	253	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
240	`ierr = nf90_get_att(cdfid, varid, "units", units)`	254	`ierr = nf90_get_att(cdfid, varid, "units", units)`
241	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	255	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
242	`if ( units.eq.'Pa' ) then`	256	`if ( units.eq.'Pa' ) then`
243	`do k=1,nakbktmp`	257	`do k=1,nakbktmp`
244	`aktmp(k) = 0.01 * aktmp(k)`	258	`aktmp(k) = 0.01 * aktmp(k)`
245	`enddo`	259	`enddo`
246	`endif`	260	`endif`
247		261
-		262	`c Check that unit of lev is in hPa - if necessary correct it`
-		263	`if ( leveltype.eq.'air_pressure' ) then`
-		264	`varname='lev'`
-		265	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
-		266	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
-		267	`ierr = nf90_get_att(cdfid, varid, "units", units)`
-		268	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
-		269	`if ( units.eq.'Pa' ) then`
-		270	`do k=1,vardim(3)`
-		271	`lev(k) = 0.01 * lev(k)`
-		272	`enddo`
-		273	`endif`
-		274	`endif`
-		275
248	`c Decide whether to swap vertical levels - highest pressure at index 1`	276	`c Decide whether to swap vertical levels - highest pressure at index 1`
249	`vertical_swap = 1`	277	`vertical_swap = 1`
-		278	`if ( leveltype.eq.'hybrid_sigma_pressure') then`
250	`if ( (aktmp(1) + bktmp(1) * 1000.).gt.`	279	`if ( (aktmp(1) + bktmp(1) * 1000.).gt.`
251	`> (aktmp(2) + bktmp(2) * 1000.) )`	280	`> (aktmp(2) + bktmp(2) * 1000.) )`
252	`>then`	281	`> then`
-		282	`vertical_swap = 0`
-		283	`endif`
-		284	`elseif ( leveltype.eq.'air_pressure') then`
-		285	`if ( lev(1).gt.lev(2) ) then`
253	`vertical_swap = 0`	286	`vertical_swap = 0`
-		287	`endif`
254	`endif`	288	`endif`
-		289	`c print*,' Vertical SWAP P -> ', vertical_swap`
255		290
256	`c Get time information (check if time is correct)`	291	`c Get time information (check if time is correct)`
257	`varname = 'time'`	292	`varname = 'time'`
258	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	293	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
259	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	294	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
260	`ierr = nf90_get_var(cdfid,varid,times)`	295	`ierr = nf90_get_var(cdfid,varid,times)`
261	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	296	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
262	`isok=0`	297	`isok=0`
263	`do i=1,vardim(4)`	298	`do i=1,vardim(4)`
264	`if (abs(time-times(i)).lt.eps) then`	299	`if (abs(time-times(i)).lt.eps) then`
265	`isok = 1`	300	`isok = 1`
266	`rtime = times(i)`	301	`rtime = times(i)`
267	`elseif (timecheck.eq.'no') then`	302	`elseif (timecheck.eq.'no') then`
268	`isok = 1`	303	`isok = 1`
269	`rtime = times(1)`	304	`rtime = times(1)`
270	`endif`	305	`endif`
271	`enddo`	306	`enddo`
272	`if ( isok.eq.0 ) then`	307	`if ( isok.eq.0 ) then`
273	`print*,' ERROR: time ',rtime,' not found on netCDF file'`	308	`print*,' ERROR: time ',rtime,' not found on netCDF file'`
274	`stop`	309	`stop`
275	`endif`	310	`endif`
276		311
277	`c Read surface pressure`	312	`c Read surface pressure`
278	`varname='PS'`	313	`varname='PS'`
279	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	314	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
280	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	315	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
281	`ierr = nf90_get_var(cdfid,varid,tmp2)`	316	`ierr = nf90_get_var(cdfid,varid,tmp2)`
282	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	317	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
283		318
284	`c Check that surface pressure is in hPa`	319	`c Check that surface pressure is in hPa`
285	`maxps = -1.e19`	320	`maxps = -1.e19`
286	`minps = 1.e19`	321	`minps = 1.e19`
287	`do i=1,vardim(1)`	322	`do i=1,vardim(1)`
288	`do j=1,vardim(2)`	323	`do j=1,vardim(2)`
289	`if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)`	324	`if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)`
290	`if (tmp2(i,j).lt.minps) minps = tmp2(i,j)`	325	`if (tmp2(i,j).lt.minps) minps = tmp2(i,j)`
291	`enddo`	326	`enddo`
292	`enddo`	327	`enddo`
293	`if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then`	328	`if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then`
294	`print*,' ERROR: surface pressre PS must be in hPa'`	329	`print*,' ERROR: surface pressre PS must be in hPa'`
295	`print*,' ',maxps,minps`	330	`print*,' ',maxps,minps`
296	`stop`	331	`stop`
297	`endif`	332	`endif`
298		333
299	`c ---- Define output of subroutine --------------------------------`	334	`c ---- Define output of subroutine --------------------------------`
300		335
301	`c If not full list of vertical levels, reduce AK,BK arrays`	336	`c If not full list of vertical levels, reduce AK,BK arrays`
302	`if ( (leveltype.eq.'hybrid_sigma_pressure').and.`	337	`if ( (leveltype.eq.'hybrid_sigma_pressure').and.`
303	`> (nakbktmp.ne.vardim(3) ) )`	338	`> (nakbktmp.ne.vardim(3) ) )`
304	`>then`	339	`>then`
305	`print*,' WARNING: only subset of vertical levels used...'`	340	`print*,' WARNING: only subset of vertical levels used...'`
306	`do k=1,vardim(3)`	341	`do k=1,vardim(3)`
307	`if ( vertical_swap.eq.1 ) then`	342	`if ( vertical_swap.eq.1 ) then`
308	`aktmp(k) = aktmp( k+nakbktmp-vardim(3) )`	343	`aktmp(k) = aktmp( k+nakbktmp-vardim(3) )`
309	`bktmp(k) = bktmp( k+nakbktmp-vardim(3) )`	344	`bktmp(k) = bktmp( k+nakbktmp-vardim(3) )`
310	`endif`	345	`endif`
311	`enddo`	346	`enddo`
312	`endif`	347	`endif`
313		348
314	`c Set the grid dimensions and constants`	349	`c Set the grid dimensions and constants`
315	`nx = vardim(1)`	350	`nx = vardim(1)`
316	`ny = vardim(2)`	351	`ny = vardim(2)`
317	`nz = vardim(3)`	352	`nz = vardim(3)`
318	`xmin = lon(1)`	353	`xmin = lon(1)`
319	`ymin = lat(1)`	354	`ymin = lat(1)`
320	`xmax = lon(nx)`	355	`xmax = lon(nx)`
321	`ymax = lat(ny)`	356	`ymax = lat(ny)`
322	`dx = (xmax-xmin)/real(nx-1)`	357	`dx = (xmax-xmin)/real(nx-1)`
323	`dy = (ymax-ymin)/real(ny-1)`	358	`dy = (ymax-ymin)/real(ny-1)`
324	`pollon = 0.`	359	`pollon = 0.`
325	`pollat = 90.`	360	`pollat = 90.`
326	`stagz = 0.`	361	`stagz = 0.`
327	`delta = xmax-xmin-360.`	362	`delta = xmax-xmin-360.`
328	`if (abs(delta+dx).lt.eps) then`	363	`if (abs(delta+dx).lt.eps) then`
329	`xmax = xmax + dx`	364	`xmax = xmax + dx`
330	`nx = nx + 1`	365	`nx = nx + 1`
331	`closear = 1`	366	`closear = 1`
332	`else`	367	`else`
333	`closear = 0`	368	`closear = 0`
334	`endif`	369	`endif`
335		370
336	`c Save the output arrays (if fid>0) - close arrays on request`	371	`c Save the output arrays (if fid>0) - close arrays on request`
337	`if ( fid.gt.0 ) then`	372	`if ( fid.gt.0 ) then`
338		373
339	`c Calculate layer pressures`	374	`c Calculate layer pressures`
-		375	`if (leveltype.eq.'hybrid_sigma_pressure' ) then`
340	`do i=1,vardim(1)`	376	`do i=1,vardim(1)`
341	`do j=1,vardim(2)`	377	`do j=1,vardim(2)`
342	`do k=1,vardim(3)`	378	`do k=1,vardim(3)`
343	`tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)`	379	`tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)`
344	`enddo`	380	`enddo`
345	`enddo`	381	`enddo`
346	`enddo`	382	`enddo`
-		383	`elseif (leveltype.eq.'air_pressure' ) then`
-		384	`do i=1,vardim(1)`
-		385	`do j=1,vardim(2)`
-		386	`do k=1,vardim(3)`
-		387	`tmp3(i,j,k)=lev(k)`
-		388	`enddo`
-		389	`enddo`
-		390	`enddo`
-		391	`endif`
347		392
348	`c Get PS - close array on demand`	393	`c Get PS - close array on demand`
349	`do j=1,vardim(2)`	394	`do j=1,vardim(2)`
350	`do i=1,vardim(1)`	395	`do i=1,vardim(1)`
351	`ps(i,j) = tmp2(i,j)`	396	`ps(i,j) = tmp2(i,j)`
352	`enddo`	397	`enddo`
353	`if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)`	398	`if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)`
354	`enddo`	399	`enddo`
355		400
356	`c Get P3 - close array on demand + vertical swap`	401	`c Get P3 - close array on demand + vertical swap`
357	`do j=1,vardim(2)`	402	`do j=1,vardim(2)`
358	`do k=1,vardim(3)`	403	`do k=1,vardim(3)`
359	`do i=1,vardim(1)`	404	`do i=1,vardim(1)`
360	`if ( vertical_swap.eq.1 ) then`	405	`if ( vertical_swap.eq.1 ) then`
361	`p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)`	406	`p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)`
362	`else`	407	`else`
363	`p3(i,j,k) = tmp3(i,j,k)`	408	`p3(i,j,k) = tmp3(i,j,k)`
364	`endif`	409	`endif`
365	`enddo`	410	`enddo`
366	`if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)`	411	`if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)`
367	`enddo`	412	`enddo`
368	`enddo`	413	`enddo`
369		414
370	`c Get AK,BK - vertical swap on demand`	415	`c Get AK,BK - vertical swap on demand`
-		416	`if ( leveltype.eq.'hybrid_sigma_pressure' ) then`
371	`do k=1,vardim(3)`	417	`do k=1,vardim(3)`
372	`if ( vertical_swap.eq.1 ) then`	418	`if ( vertical_swap.eq.1 ) then`
373	`ak(k) = aktmp(vardim(3)-k+1)`	419	`ak(k) = aktmp(vardim(3)-k+1)`
374	`bk(k) = bktmp(vardim(3)-k+1)`	420	`bk(k) = bktmp(vardim(3)-k+1)`
375	`endif`	421	`endif`
376	`enddo`	422	`enddo`
-		423	`elseif (leveltype.eq.'air_pressure' ) then`
-		424	`do k=1,vardim(3)`
-		425	`if ( vertical_swap.eq.1 ) then`
-		426	`ak(k) = lev(vardim(3)-k+1)`
-		427	`bk(k) = 0.`
-		428	`else`
-		429	`ak(k) = lev(k)`
-		430	`bk(k) = 0.`
-		431	`endif`
-		432	`enddo`
-		433	`endif`
377		434
378	`endif`	435	`endif`
379		436
380		437
381	`return`	438	`return`
382		439
383	`end`	440	`end`
384		441
385	`c ------------------------------------------------------------`	442	`c ------------------------------------------------------------`
386	`c Read wind information`	443	`c Read wind information`
387	`c ------------------------------------------------------------`	444	`c ------------------------------------------------------------`
388		445
389	`subroutine input_wind (fid,fieldname,field,time,stagz,mdv,`	446	`subroutine input_wind (fid,fieldname,field,time,stagz,mdv,`
390	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,`	447	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,`
391	`> timecheck)`	448	`> timecheck)`
392		449
393	`c Read the wind component <fieldname> from the file with identifier`	450	`c Read the wind component <fieldname> from the file with identifier`
394	`c <fid> and save it in the 3d array <field>. The vertical staggering`	451	`c <fid> and save it in the 3d array <field>. The vertical staggering`
395	`c information is provided in <stagz> and gives the reference to either`	452	`c information is provided in <stagz> and gives the reference to either`
396	`c the layer or level field from <input_grid>. A consistency check is`	453	`c the layer or level field from <input_grid>. A consistency check is`
397	`c performed to have an agreement with the grid specified by <xmin,xmax,`	454	`c performed to have an agreement with the grid specified by <xmin,xmax,`
398	`c ymin,ymax,dx,dy,nx,ny,nz>.`	455	`c ymin,ymax,dx,dy,nx,ny,nz>.`
399		456
400	`use netcdf`	457	`use netcdf`
401		458
402	`implicit none`	459	`implicit none`
403		460
404	`c Declaration of variables and parameters`	461	`c Declaration of variables and parameters`
405	`integer fid ! File identifier`	462	`integer fid ! File identifier`
406	`character*80 fieldname ! Name of the wind field`	463	`character*80 fieldname ! Name of the wind field`
407	`integer nx,ny,nz ! Dimension of fields`	464	`integer nx,ny,nz ! Dimension of fields`
408	`real field(nx,ny,nz) ! 3d wind field`	465	`real field(nx,ny,nz) ! 3d wind field`
409	`real stagz ! Staggering in the z direction`	466	`real stagz ! Staggering in the z direction`
410	`real mdv ! Missing data flag`	467	`real mdv ! Missing data flag`
411	`real xmin,xmax,ymin,ymax ! Domain size`	468	`real xmin,xmax,ymin,ymax ! Domain size`
412	`real dx,dy ! Horizontal resolution`	469	`real dx,dy ! Horizontal resolution`
413	`real time ! Time`	470	`real time ! Time`
414	`character*80 timecheck ! Either 'yes' or 'no'`	471	`character*80 timecheck ! Either 'yes' or 'no'`
415		472
416	`c Numerical and physical parameters`	473	`c Numerical and physical parameters`
417	`real eps ! Numerical epsilon`	474	`real eps ! Numerical epsilon`
418	`parameter (eps=0.001)`	475	`parameter (eps=0.001)`
419	`real notimecheck ! 'Flag' for no time check`	476	`real notimecheck ! 'Flag' for no time check`
420	`parameter (notimecheck=7.26537)`	477	`parameter (notimecheck=7.26537)`
421		478
422	`c Netcdf variables`	479	`c Netcdf variables`
423	`integer ierr`	480	`integer ierr`
424	`character*80 varname`	481	`character*80 varname`
425	`integer vardim(4)`	482	`integer vardim(4)`
426	`real varmin(4),varmax(4)`	483	`real varmin(4),varmax(4)`
427	`real stag(4)`	484	`real stag(4)`
428	`integer ndim`	485	`integer ndim`
429	`real times(10)`	486	`real times(10)`
430	`integer ntimes`	487	`integer ntimes`
431	`character*80 cstfile`	488	`character*80 cstfile`
432	`integer cstid`	489	`integer cstid`
433	`real aklay(200),bklay(200),aklev(200),bklev(200)`	490	`real aklay(200),bklay(200),aklev(200),bklev(200)`
434	`real ps(nx,ny)`	491	`real ps(nx,ny)`
435	`integer dimids (nf90_max_var_dims)`	492	`integer dimids (nf90_max_var_dims)`
436	`character*80 dimname(nf90_max_var_dims)`	493	`character*80 dimname(nf90_max_var_dims)`
437	`integer varid`	494	`integer varid`
438	`integer cdfid`	495	`integer cdfid`
439	`real,allocatable, dimension (:) :: lon,lat,lev`	496	`real,allocatable, dimension (:) :: lon,lat,lev`
440	`real,allocatable, dimension (:,:) :: tmp2`	497	`real,allocatable, dimension (:,:) :: tmp2`
441	`real,allocatable, dimension (:,:,:) :: tmp3`	498	`real,allocatable, dimension (:,:,:) :: tmp3`
442	`real,allocatable, dimension (:) :: aktmp,bktmp`	499	`real,allocatable, dimension (:) :: aktmp,bktmp`
443	`character*80 leveltype`	500	`character*80 leveltype`

Subversion Repositories lagranto.ecmwf

(root)/trunk/lib/ioinp.f – Rev 27 → 29