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

Subversion Repositories lagranto.ecmwf

(root)/trunk/lib/ioinp_cdo.f/ – Rev 46 → 47