WebSVN – lagranto.ecmwf – Diff – /trunk/lib/ioinp_cdo.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     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
      endi
 
 
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 25	Rev 27
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`	189	`c Check about the type of vertical levels`
190	`varname=dimname(3)`	190	`varname=dimname(3)`
191	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	191	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
192	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	192	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
193	`ierr = nf90_get_att(cdfid, varid, "standard_name", leveltype)`	193	`ierr = nf90_get_att(cdfid, varid, "standard_name", leveltype)`
194	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	194	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
195	`if ( (leveltype.ne.'hybrid_sigma_pressure').and.`	195	`if ( (leveltype.ne.'hybrid_sigma_pressure').and.`
196	`> (leveltype.ne.'air_pressure' ) )`	196	`> (leveltype.ne.'air_pressure' ) )`
197	`>then`	197	`>then`
198	`print*,' ERROR: input netCDF data must be on hybrid-sigma'`	198	`print*,' ERROR: input netCDF data must be on hybrid-sigma'`
199	`print*,' or air pressure levels!',trim(leveltype)`	199	`print*,' or air pressure levels!',trim(leveltype)`
200	`stop`	200	`stop`
201	`endif`	201	`endif`
202		202
203	`c Check that vardim(3)==#AK,BK for hybrid-sigmal levels`	-
204	`if ( (leveltype.eq.'hybrid_sigma_pressure').and.`	-
205	`> (dimname(3).eq.'lev') )`	-
206	`>then`	-
207	`if ( nakbktmp.ne.vardim(3) ) then`	-
208	`print*,' ERROR: for hybrid-sigma pressure levels, #AK,BK'`	-
209	`print*,' must agree with number of vertical levels'`	-
210	`print*,' ',vardim(3),nakbktmp`	-
211	`stop`	-
212	`endif`	-
213	`endif`	-
214		-
215	`c Allocate memory for reading arrays`	203	`c Allocate memory for reading arrays`
216	`allocate(tmp2(vardim(1),vardim(2)),stat=stat)`	204	`allocate(tmp2(vardim(1),vardim(2)),stat=stat)`
217	`if (stat.ne.0) print,' error allocating array tmp2 *'`	205	`if (stat.ne.0) print,' error allocating array tmp2 *'`
218	`allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)`	206	`allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)`
219	`if (stat.ne.0) print,' error allocating array tmp3 *'`	207	`if (stat.ne.0) print,' error allocating array tmp3 *'`
220	`allocate(lon(vardim(1)),stat=stat)`	208	`allocate(lon(vardim(1)),stat=stat)`
221	`if (stat.ne.0) print,' error allocating array lon *'`	209	`if (stat.ne.0) print,' error allocating array lon *'`
222	`allocate(lat(vardim(2)),stat=stat)`	210	`allocate(lat(vardim(2)),stat=stat)`
223	`if (stat.ne.0) print,' error allocating array lat *'`	211	`if (stat.ne.0) print,' error allocating array lat *'`
224	`allocate(lev(vardim(3)),stat=stat)`	212	`allocate(lev(vardim(3)),stat=stat)`
225	`if (stat.ne.0) print,' error allocating array lev *'`	213	`if (stat.ne.0) print,' error allocating array lev *'`
226	`allocate(times(vardim(4)),stat=stat)`	214	`allocate(times(vardim(4)),stat=stat)`
227	`if (stat.ne.0) print,' error allocating array times *'`	215	`if (stat.ne.0) print,' error allocating array times *'`
228	`allocate(aktmp(nakbktmp),stat=stat)`	216	`allocate(aktmp(nakbktmp),stat=stat)`
229	`if (stat.ne.0) print,' error allocating array aktmp *'`	217	`if (stat.ne.0) print,' error allocating array aktmp *'`
230	`allocate(bktmp(nakbktmp),stat=stat)`	218	`allocate(bktmp(nakbktmp),stat=stat)`
231	`if (stat.ne.0) print,' error allocating array bktmp *'`	219	`if (stat.ne.0) print,' error allocating array bktmp *'`
232		220
233	`c Get domain longitudes, latitudes and levels`	221	`c Get domain longitudes, latitudes and levels`
234	`varname = dimname(1)`	222	`varname = dimname(1)`
235	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	223	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
236	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	224	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
237	`ierr = nf90_get_var(cdfid,varid,lon)`	225	`ierr = nf90_get_var(cdfid,varid,lon)`
238	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	226	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
239	`varname = dimname(2)`	227	`varname = dimname(2)`
240	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	228	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
241	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	229	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
242	`ierr = nf90_get_var(cdfid,varid,lat)`	230	`ierr = nf90_get_var(cdfid,varid,lat)`
243	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	231	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
244	`varname = dimname(3)`	232	`varname = dimname(3)`
245	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	233	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
246	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	234	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
247	`ierr = nf90_get_var(cdfid,varid,lev)`	235	`ierr = nf90_get_var(cdfid,varid,lev)`
248	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	236	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
249		237
250	`c Get ak and bk`	238	`c Get ak and bk`
251	`varname='hyam'`	239	`varname='hyam'`
252	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	240	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
253	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	241	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
254	`ierr = nf90_get_var(cdfid,varid,aktmp)`	242	`ierr = nf90_get_var(cdfid,varid,aktmp)`
255	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	243	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
256	`varname='hybm'`	244	`varname='hybm'`
257	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	245	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
258	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	246	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
259	`ierr = nf90_get_var(cdfid,varid,bktmp)`	247	`ierr = nf90_get_var(cdfid,varid,bktmp)`
260	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	248	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
261		249
262	`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`
263	`varname='hyam'`	251	`varname='hyam'`
264	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	252	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
265	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	253	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
266	`ierr = nf90_get_att(cdfid, varid, "units", units)`	254	`ierr = nf90_get_att(cdfid, varid, "units", units)`
267	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	255	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
268	`if ( units.eq.'Pa' ) then`	256	`if ( units.eq.'Pa' ) then`
269	`do k=1,nakbktmp`	257	`do k=1,nakbktmp`
270	`aktmp(k) = 0.01 * aktmp(k)`	258	`aktmp(k) = 0.01 * aktmp(k)`
271	`enddo`	259	`enddo`
272	`endif`	260	`endif`
273		261
274	`c Check that unit of lev is in hPa - if necessary correct it`	262	`c Check that unit of lev is in hPa - if necessary correct it`
275	`if ( leveltype.eq.'air_pressure' ) then`	263	`if ( leveltype.eq.'air_pressure' ) then`
276	`varname='lev'`	264	`varname='lev'`
277	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	265	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
278	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	266	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
279	`ierr = nf90_get_att(cdfid, varid, "units", units)`	267	`ierr = nf90_get_att(cdfid, varid, "units", units)`
280	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	268	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
281	`if ( units.eq.'Pa' ) then`	269	`if ( units.eq.'Pa' ) then`
282	`do k=1,vardim(3)`	270	`do k=1,vardim(3)`
283	`lev(k) = 0.01 * lev(k)`	271	`lev(k) = 0.01 * lev(k)`
284	`enddo`	272	`enddo`
285	`endif`	273	`endif`
286	`endif`	274	`endif`
287		275
288	`c Decide whether to swap vertical levels`	276	`c Decide whether to swap vertical levels - highest pressure at index 1`
289	`vertical_swap = 1`	277	`vertical_swap = 1`
290	`if ( leveltype.eq.'hybrid_sigma_pressure') then`	278	`if ( leveltype.eq.'hybrid_sigma_pressure') then`
291	`if ( (aktmp(1) + bktmp(1) * 1000.).gt.`	279	`if ( (aktmp(1) + bktmp(1) * 1000.).gt.`
292	`> (aktmp(2) + bktmp(2) * 1000.) )`	280	`> (aktmp(2) + bktmp(2) * 1000.) )`
293	`> then`	281	`> then`
294	`vertical_swap = 0`	282	`vertical_swap = 0`
295	`endif`	283	`endif`
296	`elseif ( leveltype.eq.'air_pressure') then`	284	`elseif ( leveltype.eq.'air_pressure') then`
297	`if ( lev(1).gt.lev(2) ) then`	285	`if ( lev(1).gt.lev(2) ) then`
298	`vertical_swap = 0`	286	`vertical_swap = 0`
299	`endif`	287	`endif`
300	`endif`	288	`endif`
301	`c print*,' Vertical SWAP P -> ', vertical_swap`	289	`c print*,' Vertical SWAP P -> ', vertical_swap`
302		290
303	`c Get time information (check if time is correct)`	291	`c Get time information (check if time is correct)`
304	`varname = 'time'`	292	`varname = 'time'`
305	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	293	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
306	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	294	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
307	`ierr = nf90_get_var(cdfid,varid,times)`	295	`ierr = nf90_get_var(cdfid,varid,times)`
308	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	296	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
309	`isok=0`	297	`isok=0`
310	`do i=1,vardim(4)`	298	`do i=1,vardim(4)`
311	`if (abs(time-times(i)).lt.eps) then`	299	`if (abs(time-times(i)).lt.eps) then`
312	`isok = 1`	300	`isok = 1`
313	`rtime = times(i)`	301	`rtime = times(i)`
314	`elseif (timecheck.eq.'no') then`	302	`elseif (timecheck.eq.'no') then`
315	`isok = 1`	303	`isok = 1`
316	`rtime = times(1)`	304	`rtime = times(1)`
317	`endif`	305	`endif`
318	`enddo`	306	`enddo`
319	`if ( isok.eq.0 ) then`	307	`if ( isok.eq.0 ) then`
320	`print*,' ERROR: time ',rtime,' not found on netCDF file'`	308	`print*,' ERROR: time ',rtime,' not found on netCDF file'`
321	`stop`	309	`stop`
322	`endif`	310	`endif`
323		311
324	`c Read surface pressure`	312	`c Read surface pressure`
325	`varname='PS'`	313	`varname='PS'`
326	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	314	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
327	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	315	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
328	`ierr = nf90_get_var(cdfid,varid,tmp2)`	316	`ierr = nf90_get_var(cdfid,varid,tmp2)`
329	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	317	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
330		318
331	`c Check that surface pressure is in hPa`	319	`c Check that surface pressure is in hPa`
332	`maxps = -1.e19`	320	`maxps = -1.e19`
333	`minps = 1.e19`	321	`minps = 1.e19`
334	`do i=1,vardim(1)`	322	`do i=1,vardim(1)`
335	`do j=1,vardim(2)`	323	`do j=1,vardim(2)`
336	`if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)`	324	`if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)`
337	`if (tmp2(i,j).lt.minps) minps = tmp2(i,j)`	325	`if (tmp2(i,j).lt.minps) minps = tmp2(i,j)`
338	`enddo`	326	`enddo`
339	`enddo`	327	`enddo`
340	`if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then`	328	`if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then`
341	`print*,' ERROR: surface pressre PS must be in hPa'`	329	`print*,' ERROR: surface pressre PS must be in hPa'`
342	`print*,' ',maxps,minps`	330	`print*,' ',maxps,minps`
343	`stop`	331	`stop`
344	`endif`	332	`endif`
345		333
346	`c ---- Define output of subroutine --------------------------------`	334	`c ---- Define output of subroutine --------------------------------`
347		335
-		336	`c If not full list of vertical levels, reduce AK,BK arrays`
-		337	`if ( (leveltype.eq.'hybrid_sigma_pressure').and.`
-		338	`> (nakbktmp.ne.vardim(3) ) )`
-		339	`>then`
-		340	`print*,' WARNING: only subset of vertical levels used...'`
-		341	`do k=1,vardim(3)`
-		342	`if ( vertical_swap.eq.1 ) then`
-		343	`aktmp(k) = aktmp( k+nakbktmp-vardim(3) )`
-		344	`bktmp(k) = bktmp( k+nakbktmp-vardim(3) )`
-		345	`endif`
-		346	`enddo`
-		347	`endif`
-		348
348	`c Set the grid dimensions and constants`	349	`c Set the grid dimensions and constants`
349	`nx = vardim(1)`	350	`nx = vardim(1)`
350	`ny = vardim(2)`	351	`ny = vardim(2)`
351	`nz = vardim(3)`	352	`nz = vardim(3)`
352	`xmin = lon(1)`	353	`xmin = lon(1)`
353	`ymin = lat(1)`	354	`ymin = lat(1)`
354	`xmax = lon(nx)`	355	`xmax = lon(nx)`
355	`ymax = lat(ny)`	356	`ymax = lat(ny)`
356	`dx = (xmax-xmin)/real(nx-1)`	357	`dx = (xmax-xmin)/real(nx-1)`
357	`dy = (ymax-ymin)/real(ny-1)`	358	`dy = (ymax-ymin)/real(ny-1)`
358	`pollon = 0.`	359	`pollon = 0.`
359	`pollat = 90.`	360	`pollat = 90.`
360	`stagz = 0.`	361	`stagz = 0.`
361	`delta = xmax-xmin-360.`	362	`delta = xmax-xmin-360.`
362	`if (abs(delta+dx).lt.eps) then`	363	`if (abs(delta+dx).lt.eps) then`
363	`xmax = xmax + dx`	364	`xmax = xmax + dx`
364	`nx = nx + 1`	365	`nx = nx + 1`
365	`closear = 1`	366	`closear = 1`
366	`else`	367	`else`
367	`closear = 0`	368	`closear = 0`
368	`endif`	369	`endif`
369		370
370	`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`
371	`if ( fid.gt.0 ) then`	372	`if ( fid.gt.0 ) then`
372		373
373	`c Calculate layer pressures`	374	`c Calculate layer pressures`
374	`if (leveltype.eq.'hybrid_sigma_pressure' ) then`	375	`if (leveltype.eq.'hybrid_sigma_pressure' ) then`
375	`do i=1,vardim(1)`	376	`do i=1,vardim(1)`
376	`do j=1,vardim(2)`	377	`do j=1,vardim(2)`
377	`do k=1,vardim(3)`	378	`do k=1,vardim(3)`
378	`tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)`	379	`tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)`
379	`enddo`	380	`enddo`
380	`enddo`	381	`enddo`
381	`enddo`	382	`enddo`
382	`elseif (leveltype.eq.'air_pressure' ) then`	383	`elseif (leveltype.eq.'air_pressure' ) then`
383	`do i=1,vardim(1)`	384	`do i=1,vardim(1)`
384	`do j=1,vardim(2)`	385	`do j=1,vardim(2)`
385	`do k=1,vardim(3)`	386	`do k=1,vardim(3)`
386	`tmp3(i,j,k)=lev(k)`	387	`tmp3(i,j,k)=lev(k)`
387	`enddo`	388	`enddo`
388	`enddo`	389	`enddo`
389	`enddo`	390	`enddo`
390	`endif`	391	`endif`
391		392
392	`c Get PS - close array on demand`	393	`c Get PS - close array on demand`
393	`do j=1,vardim(2)`	394	`do j=1,vardim(2)`
394	`do i=1,vardim(1)`	395	`do i=1,vardim(1)`
395	`ps(i,j) = tmp2(i,j)`	396	`ps(i,j) = tmp2(i,j)`
396	`enddo`	397	`enddo`
397	`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)`
398	`enddo`	399	`enddo`
399		400
400	`c Get P3 - close array on demand + vertical swap`	401	`c Get P3 - close array on demand + vertical swap`
401	`do j=1,vardim(2)`	402	`do j=1,vardim(2)`
402	`do k=1,vardim(3)`	403	`do k=1,vardim(3)`
403	`do i=1,vardim(1)`	404	`do i=1,vardim(1)`
404	`if ( vertical_swap.eq.1 ) then`	405	`if ( vertical_swap.eq.1 ) then`
405	`p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)`	406	`p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)`
406	`else`	407	`else`
407	`p3(i,j,k) = tmp3(i,j,k)`	408	`p3(i,j,k) = tmp3(i,j,k)`
408	`endif`	409	`endif`
409	`enddo`	410	`enddo`
410	`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)`
411	`enddo`	412	`enddo`
412	`enddo`	413	`enddo`
413		414
414	`c Get AK,BK - vertical swap on demand`	415	`c Get AK,BK - vertical swap on demand`
415	`if ( leveltype.eq.'hybrid_sigma_pressure' ) then`	416	`if ( leveltype.eq.'hybrid_sigma_pressure' ) then`
416	`do k=1,vardim(3)`	417	`do k=1,vardim(3)`
417	`if ( vertical_swap.eq.1 ) then`	418	`if ( vertical_swap.eq.1 ) then`
418	`ak(k) = aktmp(vardim(3)-k+1)`	419	`ak(k) = aktmp(vardim(3)-k+1)`
419	`bk(k) = bktmp(vardim(3)-k+1)`	420	`bk(k) = bktmp(vardim(3)-k+1)`
420	`endif`	421	`endif`
421	`enddo`	422	`enddo`
422	`elseif (leveltype.eq.'air_pressure' ) then`	423	`elseif (leveltype.eq.'air_pressure' ) then`
423	`do k=1,vardim(3)`	424	`do k=1,vardim(3)`
424	`if ( vertical_swap.eq.1 ) then`	425	`if ( vertical_swap.eq.1 ) then`
425	`ak(k) = lev(vardim(3)-k+1)`	426	`ak(k) = lev(vardim(3)-k+1)`
426	`bk(k) = 0.`	427	`bk(k) = 0.`
427	`else`	428	`else`
428	`ak(k) = lev(k)`	429	`ak(k) = lev(k)`
429	`bk(k) = 0.`	430	`bk(k) = 0.`
430	`endif`	431	`endif`
431	`enddo`	432	`enddo`
432	`endif`	433	`endif`
433		434
434	`endif`	435	`endif`
435		436
436		437
437	`return`	438	`return`
438		439
439	`end`	440	`end`
440		441
441	`c ------------------------------------------------------------`	442	`c ------------------------------------------------------------`
442	`c Read wind information`	443	`c Read wind information`
443	`c ------------------------------------------------------------`	444	`c ------------------------------------------------------------`
444		445
445	`subroutine input_wind (fid,fieldname,field,time,stagz,mdv,`	446	`subroutine input_wind (fid,fieldname,field,time,stagz,mdv,`
446	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,`	447	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,`
447	`> timecheck)`	448	`> timecheck)`
448		449
449	`c Read the wind component <fieldname> from the file with identifier`	450	`c Read the wind component <fieldname> from the file with identifier`
450	`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`
451	`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`
452	`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`
453	`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,`
454	`c ymin,ymax,dx,dy,nx,ny,nz>.`	455	`c ymin,ymax,dx,dy,nx,ny,nz>.`
455		456
456	`use netcdf`	457	`use netcdf`
457		458
458	`implicit none`	459	`implicit none`
459		460
460	`c Declaration of variables and parameters`	461	`c Declaration of variables and parameters`
461	`integer fid ! File identifier`	462	`integer fid ! File identifier`
462	`character*80 fieldname ! Name of the wind field`	463	`character*80 fieldname ! Name of the wind field`
463	`integer nx,ny,nz ! Dimension of fields`	464	`integer nx,ny,nz ! Dimension of fields`
464	`real field(nx,ny,nz) ! 3d wind field`	465	`real field(nx,ny,nz) ! 3d wind field`
465	`real stagz ! Staggering in the z direction`	466	`real stagz ! Staggering in the z direction`
466	`real mdv ! Missing data flag`	467	`real mdv ! Missing data flag`
467	`real xmin,xmax,ymin,ymax ! Domain size`	468	`real xmin,xmax,ymin,ymax ! Domain size`
468	`real dx,dy ! Horizontal resolution`	469	`real dx,dy ! Horizontal resolution`
469	`real time ! Time`	470	`real time ! Time`
470	`character*80 timecheck ! Either 'yes' or 'no'`	471	`character*80 timecheck ! Either 'yes' or 'no'`
471		472
472	`c Numerical and physical parameters`	473	`c Numerical and physical parameters`
473	`real eps ! Numerical epsilon`	474	`real eps ! Numerical epsilon`
474	`parameter (eps=0.001)`	475	`parameter (eps=0.001)`
475	`real notimecheck ! 'Flag' for no time check`	476	`real notimecheck ! 'Flag' for no time check`
476	`parameter (notimecheck=7.26537)`	477	`parameter (notimecheck=7.26537)`
477		478
478	`c Netcdf variables`	479	`c Netcdf variables`
479	`integer ierr`	480	`integer ierr`
480	`character*80 varname`	481	`character*80 varname`
481	`integer vardim(4)`	482	`integer vardim(4)`
482	`real varmin(4),varmax(4)`	483	`real varmin(4),varmax(4)`
483	`real stag(4)`	484	`real stag(4)`
484	`integer ndim`	485	`integer ndim`
485	`real times(10)`	486	`real times(10)`
486	`integer ntimes`	487	`integer ntimes`
487	`character*80 cstfile`	488	`character*80 cstfile`

Subversion Repositories lagranto.ecmwf

(root)/trunk/lib/ioinp_cdo.f/ – Rev 25 → 27