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

Subversion Repositories lagranto.ecmwf

(root)/trunk/lib/ioinp_cdo.f – Rev 33 → 36