Subversion Repositories lagranto.ocean

c***********************************************************************
c CF-netCDF / Reading
c***********************************************************************

c     -------------------------------------------------------------------
c     cdfopn
c     -------------------------------------------------------------------

      subroutine cdfopn(filnam,cdfid,ierror)

      use netcdf
      
      character*(*) filnam
      integer       cdfid
      integer       ierror

      ierror = NF90_OPEN(TRIM(filnam),nf90_nowrite, cdfid)
      IF ( ierror /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierror)

      return
      end


c     -------------------------------------------------------------------
c     clscdf
c     -------------------------------------------------------------------

      subroutine clscdf (cdfid, ierror)

      use netcdf
      
      integer      cdfid
      integer      ierror

      ierror = NF90_CLOSE(cdfid)
      IF( ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

      end


c     -------------------------------------------------------------------
c     getdef
c     -------------------------------------------------------------------

      subroutine getdef (cdfid, varnam, ndim, misdat, 
     &                      vardim, varmin, varmax, stag, ierror)

      use netcdf

      integer        cdfid
      character *(*) varnam
      integer        ndim
      real           misdat
      integer        vardim(4)
      real           varmin(4)
      real           varmax(4)
      real           stag(4)
      integer        ierror

      real           lat(5000)
      real           lon(5000)
      real           lev(5000)
      real           tim(5000)
      integer        varid
      integer        lonid
      integer        latid
      integer        timid
      integer        dimids (nf90_max_var_dims)
      character*80   dimname(nf90_max_var_dims)
      integer        nlon
      integer        nlat
      integer        nlev
      integer        ntim
      character*(80) attname
      integer        i

c     Get <varid>
      ierror = NF90_INQ_VARID(cdfid,varnam,varid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Get <ndim>, <dimids> and <dimname>
      ierror = nf90_inquire_variable(cdfid, varid, ndims  = ndim)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_inquire_variable(cdfid, varid, 
     >                                   dimids = dimids(1:ndim))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      do i=1,ndim
         ierror = nf90_inquire_dimension(cdfid, dimids(i), 
     >                               name = dimname(i) )
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      enddo

      if ( (ndim.ne.3).and.(ndim.ne.4).and.(ndim.ne.1) ) then
          print*,' ERROR: getdef only support 3d, 4d and 1d level arrays'
         print*,'1D=lev (ORA); 3D=time,lat,lon; 4D=time,level,lat,lon'
          stop
      endif

c     Get the ID for <longitude>, <latitude>, <level> and <time>
      if (ndim.eq.4) then
        ierror = nf90_inq_dimid(cdfid,dimname(1), lonid)
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
        ierror = nf90_inq_dimid(cdfid,dimname(2), latid)
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
        ierror = nf90_inq_dimid(cdfid,dimname(3), levid)
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
        ierror = nf90_inq_dimid(cdfid,dimname(4), timid)
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      else if ( ndim.eq.3 ) then
        ierror = nf90_inq_dimid(cdfid,dimname(1), lonid)
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
        ierror = nf90_inq_dimid(cdfid,dimname(2), latid)
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
        ierror = nf90_inq_dimid(cdfid,dimname(3), timid)
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      else if ( ndim.eq.1 ) then
        ierror = nf90_inq_dimid(cdfid,dimname(1), levid)
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      endif

c     Read lon and set <vardim(1)>, <varmin(1)>, <varmax(1)>, <stag(1)>
      if ( ndim.eq.4 .or. ndim.eq.3 ) then
        ierror = nf90_inquire_dimension(cdfid, lonid, len = nlon)
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
        ierror = NF90_INQ_VARID(cdfid,dimname(1),lonid)
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
        ierror = NF90_GET_VAR(cdfid,lonid,lon(1:nlon))
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 
        vardim(1) = nlon
        varmin(1) = lon( 1    )
        varmax(1) = lon( nlon )
        stag  (1) = 0.
      else if ( ndim.eq.1 ) then
        vardim(1) = 1
        varmin(1) = 1
        varmax(1) = 1
        stag  (1) = 0.
      end if

c     Read lat and set <vardim2(2)>, <varmin(2)>, <varmax(2)>, <stag(2)>
      if ( ndim.eq.4 .or. ndim.eq.3 ) then
        ierror = nf90_inquire_dimension(cdfid, latid, len = nlat)
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
        ierror = NF90_INQ_VARID(cdfid,dimname(2),latid)
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
        ierror = NF90_GET_VAR(cdfid,latid,lat(1:nlat))
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 
        vardim(2) = nlat
        varmin(2) = lat( 1    )
        varmax(2) = lat( nlat )
        stag  (2) = 0. 
      else if ( ndim.eq.1 ) then
        vardim(2) = 1
        varmin(2) = 1
        varmax(2) = 1
        stag  (2) = 0.
      end if
      
c     Read lev and set <vardim2(3)>, <varmin(3)>, <varmax(3)>, <stag(3)>
      if ( ndim.eq.4 ) then
         ierror = nf90_inquire_dimension(cdfid, levid, len = nlev)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = NF90_INQ_VARID(cdfid,dimname(3),levid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = NF90_GET_VAR(cdfid,levid,lev(1:nlev))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 
         vardim(3) = nlev
         varmin(3) = 1050.
         varmax(3) = 0.
         stag  (3) = -0.5
      else if (ndim.eq.3) then
         vardim(3) = 1
         varmin(3) = 1
         varmax(3) = 1
         stag  (3) = 0.
      else if ( ndim.eq.1 ) then ! valid for 1d lev array
         ierror = nf90_inquire_dimension(cdfid, levid, len = nlev)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = NF90_INQ_VARID(cdfid,dimname(1),levid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = NF90_GET_VAR(cdfid,levid,lev(1:nlev))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
        vardim(3) = nlev
        varmin(3) = lev(1)
        varmax(3) = lev(nlev)
        stag  (3) = 0.
      endif 

c     Read time and set <vardim2(4)>, <varmin(4)>, <varmax(4)>, <stag(4)>
      if ( ndim.eq.4 ) then
         ierror = nf90_inquire_dimension(cdfid, timid, len = ntim)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = NF90_INQ_VARID(cdfid,dimname(4),timid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = NF90_GET_VAR(cdfid,timid,tim(1:ntim))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 
         vardim(4) = ntim
         varmin(4) = tim( 1    )
         varmax(4) = tim( ntim )
         stag  (4) = 0.0
      else if ( ndim.eq.3 ) then 
        ierror = nf90_inquire_dimension(cdfid, timid, len = ntim)
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
        ierror = NF90_INQ_VARID(cdfid,dimname(3),timid)
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
        ierror = NF90_GET_VAR(cdfid,timid,tim(1:ntim))
        IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 
        vardim(4) = ntim
        varmin(4) = tim( 1    )
        varmax(4) = tim( ntim )
        stag  (4) = 0.0
      else if ( ndim.eq.1 ) then
        vardim(4) = 1
        varmin(4) = 1
        varmax(4) = 1
        stag  (4) = 0.     
      endif
      
c     Read missing data value <misdat>
      ierror  = nf90_get_att(cdfid, varid, "_FillValue", misdat)
      IF (ierror /= nf90_NoErr) THEN
         misdat = -999.99
         ierror = nf90_NoErr
      ENDIF

c     Do some checks
      if ( varmax(1).lt.varmin(1) ) then
         print*,' Error: Wrong order of longitudes'
         ierror = 1
         return
      endif
      if ( varmax(2).lt.varmin(2) ) then
         print*,' Error: Wrong order of latitudes'
         ierror = 1
         return
      endif

         
      end


c     -------------------------------------------------------------------
c     getgrid
c     -------------------------------------------------------------------

      subroutine getgrid(cdfid,dx,dy,ierror)

      use netcdf

      integer   cdfid
      real      dx
      real      dy
      integer   ierror

      integer   lonid
      integer   latid
      real      lat(5000)
      real      lon(5000)
      integer   nlon
      integer   nlat

c     Get length of lon and lat fields
      ierror = nf90_inq_dimid(cdfid,'lon', lonid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_inq_dimid(cdfid,'lat', latid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_inquire_dimension(cdfid, lonid, len = nlon)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_inquire_dimension(cdfid, latid, len = nlat)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Read lon and calculate <dx>
      ierror = nf90_inq_varid(cdfid,'lon', lonid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = NF90_GET_VAR(cdfid,lonid,lon(1:nlon))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 
      dx = ( lon(nlon) - lon(1) ) / ( real(nlon) - 1.)

c     Read lat and calculate <dy>
      ierror = nf90_inq_varid(cdfid,'lat', latid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = NF90_GET_VAR(cdfid,latid,lat(1:nlat))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 
      dy = ( lat(nlat) - lat(1) ) / ( real(nlat) - 1.)

      end


c     -------------------------------------------------------------------
c     getpole
c     -------------------------------------------------------------------

      subroutine getpole(cdfid,pollon,pollat,ierror)

      use netcdf

      integer   cdfid
      real      pollon
      real      pollat
      integer   ierror

      integer   varid
      integer   ncvid

      ierror  = nf90_get_att(cdfid,nf90_global, "pollon", pollon)
      IF(ierror /= nf90_NoErr) pollon = 0.
      ierror  = nf90_get_att(cdfid,nf90_global, "pollat", pollat)
      IF(ierror /= nf90_NoErr) pollat = 90.

      end


c     -------------------------------------------------------------------
c     getcfn
c     -------------------------------------------------------------------

      subroutine getcfn(cdfid,cfn,ierror)
      
      use netcdf

      integer      cdfid
      character*80 cfn
      integer      ierror

      ierror  = nf90_get_att(cdfid,nf90_global,
     >                    "constants_file_name ", cfn(1:80) )
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

      end


c     -------------------------------------------------------------------
c     getdat
c     -------------------------------------------------------------------

      subroutine getdat(cdfid, varnam, time, level, dat, ierror)

      use netcdf

      integer       cdfid
      character*(*) varnam
      real          time
      integer       level
      real          dat(*)
      integer       ierror

      integer       varid
      integer       ndim
      real          misdat
      integer       vardim(4)
      real          varmin(4)
      real          varmax(4)
      real          stag(4)
      integer       timid
      integer       lonid
      integer       latid
      integer       ntim
      real          tim(5000)
      integer       dimids (nf90_max_var_dims)
      character*80  dimname(nf90_max_var_dims)
      integer       itime
      integer       corner(4)
      integer       count (4)
      real,allocatable, dimension (:,:,:,:) :: tmp
      real,allocatable, dimension (:,:,:)   :: tmp3
      real,allocatable, dimension (:,:)     :: tmp2
      real,allocatable, dimension (:)       :: tmp1
      integer       stat
      integer       i,j,k
      character*80  dir

c     Get variable definition
      call getdef (cdfid, varnam, ndim, misdat, 
     &                vardim, varmin, varmax, stag, ierror)

      if ( (ndim.ne.3).and.(ndim.ne.4).and.(ndim.ne.1) ) then
         print*,' ERROR: getdat only support 1d, 3d and 4d arraysd'
         print*,' 1D=lev (ORA); 3D=time,lat,lon; 4D=time,level,lat,lon'
         stop
      endif

c     Get <varid>
      ierror = NF90_INQ_VARID(cdfid,varnam,varid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Get <dimids> and <dimname>
      ierror = nf90_inquire_variable(cdfid, varid, 
     >                                   dimids = dimids(1:ndim))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      do i=1,ndim
         ierror = nf90_inquire_dimension(cdfid, dimids(i), 
     >                               name = dimname(i) )
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      enddo

c     Get the ID for <longitude>, <latitude>, <level> and <time>
      if ( ndim.eq.4) then 
         ierror = nf90_inq_dimid(cdfid,dimname(1), lonid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_inq_dimid(cdfid,dimname(2), latid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_inq_dimid(cdfid,dimname(3), levid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_inq_dimid(cdfid,dimname(4), timid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      else if (ndim.eq.3)then
         ierror = nf90_inq_dimid(cdfid,dimname(1), lonid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_inq_dimid(cdfid,dimname(2), latid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_inq_dimid(cdfid,dimname(3), timid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      else if (ndim.eq.1)then
         ierror = nf90_inq_dimid(cdfid,dimname(1), levid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      endif

c     Get all times for variable
      if ( ndim.eq.4 )  then
         ierror = nf90_inquire_dimension(cdfid, timid, len = ntim)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = NF90_INQ_VARID(cdfid,dimname(4),timid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = NF90_GET_VAR(cdfid,timid,tim(1:ntim))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 
      elseif ( ndim.eq.3 ) then
         ierror = nf90_inquire_dimension(cdfid, timid, len = ntim)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = NF90_INQ_VARID(cdfid,dimname(3),timid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = NF90_GET_VAR(cdfid,timid,tim(1:ntim))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 
         
      elseif ( ndim.eq.1 ) then         
         ntim=1
         tim(1:ntim)=time         
      endif
      
c     Find the appropriate time index
      itime = 0
      do i=1,ntim
         if ( tim(i).eq.time) itime = i
      enddo
      if ( itime.eq.0 ) then
         ierror = 1
         print*,'Error: Unknown time in getdat'
         return
      endif
      

c     Define the data volume to be read
      if ( ndim.eq.4 .or. ndim.eq.3) then

      corner(1) = 1
      count (1) = vardim(1)
      
      corner(2) = 1
      count (2) = vardim(2)
      
      end if
      
      if ( ndim.eq.4 ) then
         if (level.eq.0) then
            corner(3) = 1
            count (3) = vardim(3)
         elseif ( (level.ge.1).and.(level.le.vardim(3)) ) then
            corner(3) = vardim(3)-level+1
            count (3) = 1
         else
            ierror = 1
            return
         endif

         corner(4) = itime
         count (4) = 1

      elseif ( ndim.eq.3 ) then
         
         corner(3) = itime
         count (3) = 1
         
      elseif (ndim.eq.1) then
      
        corner(3)  = 1
        count(3)   = vardim(3)
      
      endif

c     Read data
       if ( ndim.eq.4 ) then
        allocate(tmp(count(1),count(2),count(3),count(4)),stat=stat)
        ierror = nf90_get_var(cdfid,varid,tmp,
     >           start = (/ corner(1),corner(2),corner(3),corner(4) /),
     >           count = (/ count (1),count (2),count (3),count (4) /) ) 

      else if ( ndim.eq.3 ) then
        allocate(tmp3(count(1),count(2),count(3)),stat=stat)
        allocate(tmp(count(1),count(2),1,count(3)),stat=stat)
        ierror = nf90_get_var(cdfid,varid,tmp3,
     >           start = (/ corner(1),corner(2),corner(3) /),
     >           count = (/ count (1),count (2),count (3) /) )
               do j = 1, count(2)
            do i = 1, count(1)
               tmp(i,j,1,1) = tmp3(i,j,1)
            enddo
           enddo
        count(4) = count(3)
        count(3) = 1 
        
       else if ( ndim.eq.1 ) then
       
        allocate(tmp1(count(3)),stat=stat)
        ierror = nf90_get_var(cdfid,varid,tmp1,
     >         start = (/ corner(3) /),
     >         count = (/ count (3)/) ) 
     


      endif

c     Copy data to output array: change orientation of vertical axis
c      if ( ndim.eq.4 .or. ndim.eq.3) then

c      do i=1,count(1)
c         do j=1,count(2)
c            if ( count(3).gt.1) then

c               do k=1,count(3)
c                     dat(i+count(1)*(j-1)+count(1)*count(2)*(k-1)) = 
c     >                                      tmp(i,j,count(3)-k+1,1)
c               enddo

c            else 
c               dat(i+count(1)*(j-1)) = tmp(i,j,1,1)
c            endif
c
c         enddo
c        enddo

c     Copy data to output array: don't change orientation of vertical axis
      if ( ndim.eq.4 .or. ndim.eq.3) then

      do i=1,count(1)
         do j=1,count(2)
            if ( count(3).gt.1) then

               do k=1,count(3)
                     dat(i+count(1)*(j-1)+count(1)*count(2)*(k-1)) = 
     >                                      tmp(i,j,k,1)
               enddo

            else 
               dat(i+count(1)*(j-1)) = tmp(i,j,1,1)
            endif

         enddo
       enddo
          
      else if  (ndim.eq.1) then ! for ORA
        do i=1,count(3)
            dat(i) = tmp1(i)
        end do     
      
      end if
               
c     Deallocate memory
      deallocate(tmp,stat=stat)

      end

c     -------------------------------------------------------------------
c     gettimes
c     -------------------------------------------------------------------

      subroutine gettimes(cdfid,times,ntimes,ierror)

      use netcdf
      
      integer cdfid
      real    times(*)
      integer ntimes
      integer ierror
      
      integer timid

c     Get all times 
      ierror = nf90_inq_dimid(cdfid,'time', timid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_inquire_dimension(cdfid, timid, len = ntimes)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = NF90_INQ_VARID(cdfid,'time',timid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = NF90_GET_VAR(cdfid,timid,times(1:ntimes))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 

      end

c     -------------------------------------------------------------------
c     getlevs
c     -------------------------------------------------------------------

      subroutine getlevs(cdfid,nlev,aklev,bklev,aklay,bklay,ierror)

      use netcdf
      
      integer     cdfid
      integer     nlev
      real        aklev(nlev)
      real        bklev(nlev)
      real        aklay(nlev)
      real        bklay(nlev)
      integer     ierror

      integer      aklayid
      integer      bklayid
      integer      aklevid
      integer      bklevid
      integer      dimids (1)
      character*80 dimname(1)
      integer      nz
      integer      nz1
      integer      dimid
      real         aklev1(nlev+1)
      real         bklev1(nlev+1)
      character*80 dir
      real         lev0,lev1
      character*80 unit
      real         tmp(nlev)

c     Get <varid>
      ierror = NF90_INQ_VARID(cdfid,'hyam',aklayid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = NF90_INQ_VARID(cdfid,'hybm',bklayid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = NF90_INQ_VARID(cdfid,'hyai',aklevid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = NF90_INQ_VARID(cdfid,'hybi',bklevid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Get dimension of the two arrays
      ierror = nf90_inquire_variable(cdfid, aklayid, dimids = dimids)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_inquire_dimension(cdfid, dimids(1),name =dimname(1))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_inq_dimid(cdfid,dimname(1), dimid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_inquire_dimension(cdfid, dimid, len = nz)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

      ierror = nf90_inquire_variable(cdfid, aklevid, dimids = dimids)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_inquire_dimension(cdfid, dimids(1),name =dimname(1))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_inq_dimid(cdfid,dimname(1), dimid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_inquire_dimension(cdfid, dimid, len = nz1)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Check dimensions
      if ( (nz.ne.nlev).or.(nz1.ne.(nlev+1)) ) then
         print*,' Error: incompatible number of levels (getlevs)'
         ierror = 1
         return
      endif
      
c     Get levels
      ierror = NF90_GET_VAR(cdfid,aklayid,aklay(1:nz))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 
      ierror = NF90_GET_VAR(cdfid,bklayid,bklay(1:nz))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 
      ierror = NF90_GET_VAR(cdfid,aklevid,aklev1(1:nz1))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 
      ierror = NF90_GET_VAR(cdfid,bklevid,bklev1(1:nz1))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 

c     Set the appropriate levels
      do i=1,nz
         aklev(i) = aklev1(i)
         bklev(i) = bklev1(i)
      enddo

c     Adjust unit (in hPa)
      ierror = NF90_INQ_VARID(cdfid,'hyam',aklayid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror  = nf90_get_att(cdfid, aklayid, "units", unit)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      if ( unit.eq.'Pa' ) then
         do i=1,nz
            aklay(i) = 0.01 * aklay(i)
         enddo
      endif

      ierror = NF90_INQ_VARID(cdfid,'hyai',aklevid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror  = nf90_get_att(cdfid, aklevid, "units", unit)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      if ( unit.eq.'Pa' ) then
         do i=1,nz
            aklev(i) = 0.01 * aklev(i)
         enddo
      endif

c     Adjust index order (up)
      do i=1,nz
         tmp(i) = aklev(i)
      enddo
      do i=1,nz
         aklev(i) = tmp(nz-i+1)
      enddo
      
      do i=1,nz
         tmp(i) = aklay(i)
      enddo
      do i=1,nz
         aklay(i) = tmp(nz-i+1)
      enddo
      
      do i=1,nz
         tmp(i) = bklev(i)
      enddo
      do i=1,nz
         bklev(i) = tmp(nz-i+1)
      enddo
         
      do i=1,nz
         tmp(i) = bklay(i)
      enddo
      do i=1,nz
         bklay(i) = tmp(nz-i+1)
      enddo
         
      end


c     -------------------------------------------------------------------
c     getvars
c     -------------------------------------------------------------------

      subroutine getvars(cdfid,nvars,vnam,ierror)

      use netcdf
      
      integer           cdfid
      integer           nvars
      character*80      vnam(*)
      integer           ierror

      integer           i

c     Get number of variables
      ierror = nf90_inquire(cdfid, nVariables = nvars)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      
c     Get variable names
      do i=1,nvars
          ierror = nf90_inquire_variable(cdfid, i, name = vnam(i) )
          IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      enddo

      return
      end


c     -------------------------------------------------------------------
c     getstart
c     -------------------------------------------------------------------

      subroutine getstart(cdfid,idate,ierror)
      
      use netcdf
      
      integer cdfid
      integer idate(5)
      integer ierror

      integer varid

c     Read data      
      ierror  = nf90_get_att(cdfid,nf90_global, "starty", idate(1))
      IF(ierror /= nf90_NoErr) idate(1) = 0
      ierror  = nf90_get_att(cdfid,nf90_global, "startm", idate(2))
      IF(ierror /= nf90_NoErr) idate(2) = 0
      ierror  = nf90_get_att(cdfid,nf90_global, "startd", idate(3))
      IF(ierror /= nf90_NoErr) idate(3) = 0
      ierror  = nf90_get_att(cdfid,nf90_global, "starth", idate(4))
      IF(ierror /= nf90_NoErr) idate(4) = 0
      ierror  = nf90_get_att(cdfid,nf90_global, "starts", idate(5))
      IF(ierror /= nf90_NoErr) idate(5) = 0

      end


c***********************************************************************
c CF-netCDF / Writing
c***********************************************************************

c     -------------------------------------------------------------------
c     crecdf
c     -------------------------------------------------------------------

      subroutine crecdf (filnam, cdfid, phymin, phymax, ndim, cfn, 
     &                      ierror) 
      
      use netcdf
      
      character*(* ) filnam
      integer        cdfid
      real           phymin(4)
      real           phymax(4)
      integer        ndim
      character*(*)  cfn
      integer        ierror

c     Create the netCDF file
      ierror = nf90_create(path  = filnam, 
     >                     cmode = nf90_clobber, 
     >                     ncid  = cdfid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Global attributes
      ierror = nf90_put_att(cdfid, nf90_global, 
     >                            "Conventions", "CF")
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, nf90_global, 
     >                             "constants_file_name", trim(cfn))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, nf90_global, 
     >                             "institution", 'IACETH')
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, nf90_global, 
     >                             "lonmin", real(phymin(1)))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, nf90_global, 
     >                             "lonmax", real(phymax(1)))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, nf90_global, 
     >                             "latmin", real(phymin(2)))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, nf90_global, 
     >                             "latmax", real(phymax(2)))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, nf90_global, 
     >                             "levmin", real(phymin(3)))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, nf90_global, 
     >                             "levmax", real(phymax(3)))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     End definition mode
      ierror = nf90_enddef(cdfid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

      end


c     -------------------------------------------------------------------
c     putdef
c     -------------------------------------------------------------------

      subroutine putdef (cdfid, newvar, ndimension, misdat, 
     &                      vardim, varmin, varmax, stag, ierror)

      use netcdf
      
      integer        cdfid
      character*(*)  newvar
      integer        ndimension
      real           misdat
      integer        vardim(4)
      real           varmin(4)
      real           varmax(4)
      real           stag(4)
      integer        ierror

      real          eps
      parameter     (eps=0.0001)
      
      integer       ndim
      integer       nvar
      character*80  dimname(100)
      character*80  varname(100)
      integer       dimid
      integer       dimids(4)
      integer       nval
      real          val(1000)
      integer       varid
      integer       varids(4)
      integer       writelon,writelat,writelev,writetim
      character*80  lonname,latname,levname,timname

c     Unit/name table
      character*80  tablefile
      parameter     (tablefile='nil')
      character*80  table_varname (1000)
      character*80  table_longname(1000)
      character*80  table_unit    (1000)
      integer       fid
      integer       ntable

c     Read table information
      if ( tablefile.ne.'nil' ) then
         ntable = 0
         fid    = 10
         open(fid,file=tablefile)
 100      ntable = ntable + 1
          read(fid,*) table_varname (ntable),
     >                        table_longname(ntable),
     >                        table_unit    (ntable)

          print*,table_varname(ntable),table_longname(ntable),
     >             table_unit(ntable)
          
          goto 100
 101     close(fid)
      endif

c     Set the default names for dimensions
      lonname = 'lon'
      latname = 'lat'
      levname = 'lev'
      timname = 'time'

c     Get number of dimensions and of variables
      ierror = nf90_inquire(cdfid, nDimensions    = ndim,
     >                             nVariables     = nvar )
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Get a list of all dimensions
      if ( ndim.gt.0 ) then
         do i=1,ndim
            ierror = nf90_inquire_dimension(cdfid,i,name = dimname(i))
            IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         enddo
      endif

c     Get a list of all variables
      if ( nvar.gt.0 ) then
         do i=1,nvar
            ierror = nf90_inquire_variable(cdfid,i,name = varname(i))
            IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         enddo
      endif

c     Check whether the new variable name is in conflict with the others
      do i=1,nvar
         if ( newvar.eq.varname(i) ) then
            print*,trim(newvar),' already exists on the file'
            stop
         endif
      enddo

c     Set IDs for dimensions 
      do i=1,4
         dimids(i) = -1
      enddo
      writelon = 0
      writelat = 0
      writelev = 0
      writetim = 0

c     Check whether there is already a time dimension (tim)
      ierror = nf90_inquire(cdfid, unlimitedDimId = dimids(4) )
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Check whether existing dimensions can be used
      do i=1,ndim

         if (i.ne.dimids(4)) then

c        Read dimension values
         ierror = nf90_inquire_dimension(cdfid,i,len = nval)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)         
         ierror = NF90_INQ_VARID(cdfid,dimname(i),varid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 
         ierror = NF90_GET_VAR(cdfid,varid,val(1:nval))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 

c        Check whether it can be used for the first dimension (lon)
         if ( (vardim(1).eq.nval).and.
     >        (abs(varmin(1)-val(1)   ).lt.eps).and.
     >        (abs(varmax(1)-val(nval)).lt.eps) )
     >   then
            dimids(1) = i
         endif

c        Check whether it can be used for the second dimension (lat)
         if ( (vardim(2).eq.nval).and.
     >        (abs(varmin(2)-val(1)   ).lt.eps).and.
     >        (abs(varmax(2)-val(nval)).lt.eps) )
     >   then
            dimids(2) = i
         endif
         
c        Check whether it can be used for the third dimension (lev)
         if ( (vardim(3).eq.nval) ) then
            dimids(3) = i
         endif
         
         endif

      enddo

c     Decide what dimension names to use
      if ( dimids(1).eq.-1 ) then
         do i=1,ndim
            if ( dimname(i).eq.lonname ) then
               lonname = trim(lonname)//'_'//trim(newvar)
            endif
         enddo
      endif
      if ( dimids(2).eq.-1 ) then
         do i=1,ndim
            if ( dimname(i).eq.latname ) then
               latname = trim(latname)//'_'//trim(newvar)
            endif
         enddo
      endif
      if ( dimids(3).eq.-1 ) then
         do i=1,ndim
            if ( dimname(i).eq.levname ) then
               levname = trim(levname)//'_'//trim(newvar)
            endif
         enddo
      endif

c     Enter definition mode
      ierror = nf90_redef(cdfid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Define dimensions and corresponding variables if necessary
      if ( dimids(1).eq.-1 ) then
         dimid = ndim + 1
         ndim  = ndim + 1
         ierror = nf90_def_dim(cdfid, lonname, vardim(1), dimids(1))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_def_var(cdfid, lonname,
     >                                 NF90_FLOAT,dimids(1),varids(1))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(1), 
     >                               "long_name", "longitude") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(1), 
     >                               "units", "degrees_east") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(1), 
     >                               "standard_name", "longitude") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(1), 
     >                               "axis", "X") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         writelon = 1
      endif
      if ( dimids(2).eq.-1 ) then
         dimid = ndim + 1
         ndim  = ndim + 1
         ierror = nf90_def_dim(cdfid, latname, vardim(2), dimids(2))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_def_var(cdfid, latname,
     >                                NF90_FLOAT,dimids(2),varids(2))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(2), 
     >                               "long_name", "latitude") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(2), 
     >                               "units", "degrees_north") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(2), 
     >                               "standard_name", "latitude") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(2), 
     >                               "axis", "Y") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         writelat = 1
      endif
      if ( dimids(3).eq.-1 ) then
         dimid = ndim + 1
         ndim  = ndim + 1
         ierror = nf90_def_dim(cdfid, levname, vardim(3), dimids(3))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_def_var(cdfid, levname,
     >                                NF90_FLOAT,dimids(3),varids(3))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(3), 
     >                 "long_name", "hybrid level at layer midpoints") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(3), 
     >                 "units", "level") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(3), 
     >                 "standard_name", "hybrid_sigma_pressure") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(3), 
     >                 "positive", "down") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(3), 
     >                 "formula", "hyam hybm (mlev=hyam+hybm*aps)") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(3), 
     >                 "formula_terms", "ap: hyam b: hybm ps: aps") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         writelev = 1
      endif
      if ( dimids(4).eq.-1 ) then
         dimid = ndim + 1
         ndim  = ndim + 1
         ierror = nf90_def_dim(cdfid,timname, nf90_unlimited, dimids(4))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_def_var(cdfid, timname,
     >                                NF90_FLOAT,dimids(4),varids(4))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(4), 
     >                 "units", "hours") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_att(cdfid, varids(4), 
     >                 "calendar", "proleptic_gregorian") 
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         
         writetim = 1
      endif

c     Define the variable
      ierror = nf90_def_var(cdfid,newvar,
     >                      NF90_FLOAT,dimids,varid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Write variable attributes
      ierror = nf90_put_att(cdfid, varid, 
     >                 "long_name", "unknown (please add with NCO)") 
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, varid, 
     >                 "units", "unknown (please add with NCO)") 
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, varid, 
     >                 "_FillValue", real(-999.99) ) 
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)


c     Leave definition mode
      ierror = nf90_enddef(cdfid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Write dimension arrays
      if ( writelon.eq.1 ) then
         do i=1,vardim(1)
            val(i) = varmin(1) + real(i-1) *
     >                 ( varmax(1) - varmin(1) ) / real(vardim(1)-1)
         enddo
         ierror = nf90_inq_varid(cdfid,lonname,varid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_var(cdfid, varid, val(1:vardim(1)))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      endif
      if ( writelat.eq.1 ) then
         do i=1,vardim(2)
            val(i) = varmin(2) + real(i-1) *
     >                 ( varmax(2) - varmin(2) ) / real(vardim(2)-1)
         enddo
         ierror = nf90_inq_varid(cdfid,latname,varid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_var(cdfid, varid, val(1:vardim(2)))
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      endif

      if ( writelev.eq.1 ) then
         do i=1,vardim(3)
            val(i) = real(i)
         enddo
         ierror = nf90_inq_varid(cdfid,levname,varid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_put_var(cdfid, varid, val(1:vardim(3)) )
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      endif

      end


c     -------------------------------------------------------------------
c     puttimes
c     -------------------------------------------------------------------

      subroutine puttimes(cdfid,times,ntimes,ierror)

      use netcdf
      
      integer   cdfid
      real      times(*)
      integer   ntimes
      integer   ierror

      integer   timid
      integer   otimes

c     Get number of times 
      ierror = nf90_inq_dimid(cdfid,'time', timid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_inquire_dimension(cdfid, timid, len = otimes)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      
      if (otimes.lt.ntimes) then
         print *,'Warning: puttimes is extending times-array'
      else if (nfiltim.gt.ntimes) then
         print *,'Warning: puttimes does not cover range of times-array'
      endif

c     Write times
      ierror = NF90_INQ_VARID(cdfid,'time',timid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = NF90_PUT_VAR(cdfid,timid,times(1:ntimes))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 

c     Sync
      ierror = nf90_sync(cdfid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 

      end


c     -------------------------------------------------------------------
c     putstart
c     -------------------------------------------------------------------

      subroutine putstart(cdfid,idate,ierror)

      use netcdf

      integer      cdfid
      integer      idate(5)
      integer      ierror

c     Enter definition mode
      ierror = nf90_redef(cdfid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Write data      
      ierror  = nf90_put_att(cdfid,nf90_global, "starty", idate(1))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror  = nf90_put_att(cdfid,nf90_global, "startm", idate(2))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror  = nf90_put_att(cdfid,nf90_global, "startd", idate(3))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror  = nf90_put_att(cdfid,nf90_global, "starth", idate(4))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror  = nf90_put_att(cdfid,nf90_global, "starts", idate(5))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Exit definition mode
      ierror = nf90_enddef(cdfid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Sync
      ierror = nf90_sync(cdfid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 

      end


c     -------------------------------------------------------------------
c     cdfwopn
c     -------------------------------------------------------------------

      subroutine cdfwopn(filnam,cdfid,ierror)
      
      use netcdf
      
      character*(*) filnam
      integer       cdfid
      integer       ierror

      ierror = NF90_OPEN(TRIM(filnam),nf90_write, cdfid)
      IF ( ierror /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierror)

      return
      end


c     -------------------------------------------------------------------
c     putdat
c     -------------------------------------------------------------------

      subroutine putdat(cdfid, varnam, time, level, dat, ierror)

      use netcdf
      
      integer        cdfid
      character*(*)  varnam
      real           time
      integer        level
      real           dat(*)
      integer        ierror

      real           eps
      parameter      (eps=0.0001)

      integer        ndims
      real           misdat
      integer        vardim(4)
      real           varmin(4)
      real           varmax(4)
      real           stag(4)
      integer        varid
      integer        ntimes
      real           times(100)
      integer        itime
      integer        corner(4)
      integer        count(4)
      integer        stat
      real,allocatable, dimension (:,:,:,:) :: tmp

c     Get definitions of data
      call getdef (cdfid, trim(varnam), ndims, misdat, 
     &                vardim, varmin, varmax, stag, ierr)
      if (ierr.ne.0)  print *,'*ERROR* in getdef in putdat'

c     Get id of the variable
      ierror = NF90_INQ_VARID(cdfid,varnam,varid)
      IF (ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Get a list of all times 
      call gettimes(cdfid,times,ntimes,ierror)

c     Get the time index
      itime = -1
      do i=1,ntimes
         if ( abs(time-times(i)).lt.eps) itime = i
      enddo

c     Start a new time step, if necessary
      if ( itime.eq.-1 ) then
         ntimes        = ntimes + 1
         times(ntimes) = time
         itime         = ntimes
         call puttimes(cdfid,times,ntimes,ierror)
      endif

C     Define data volume to write on the NetCDF file in index space
      corner(1) = 1               
      corner(2) = 1
      count(1)  = vardim(1)       
      count(2)  = vardim(2)
      if (level.eq.0) then
        corner(3) = 1
        count(3)  = vardim(3)
      else
        corner(3) = vardim(3)-level+1
        count(3)  = 1
      endif
      corner(4) = itime
      count(4)  = 1

c     Allocate temporary array
      allocate(tmp(count(1),count(2),count(3),count(4)),stat=stat)

c     Copy data to output array: change orientation of vertical axis
      if ( count(3).gt.1 ) then
         do i=1,count(1)
            do j=1,count(2)
               do k=1,count(3)
                  tmp(i,j,count(3)-k+1,1) =
     >               dat(i+count(1)*(j-1)+count(1)*count(2)*(k-1)) 
               enddo
            enddo
         enddo
      else
         do i=1,count(1)
            do j=1,count(2)
               tmp(i,j,1,1) = dat(i+count(1)*(j-1)) 
            enddo
         enddo
      endif
               
c     Write data
      ierror = nf90_put_var(cdfid, varid, tmp,
     >           start = (/ corner(1),corner(2),corner(3),corner(4) /),
     >           count = (/ count (1),count (2),count (3),count (4) /) ) 

c     Deallocate memory
      deallocate(tmp,stat=stat)

c     Sync
      ierror = nf90_sync(cdfid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 

      end


c     -------------------------------------------------------------------
c     wricst
c     -------------------------------------------------------------------

      subroutine wricst(cstnam,datar,aklev,bklev,aklay,bklay,stdate)

      use netcdf
      
      character*80  cstnam
      integer       datar(14)
      real          aklev(1000)
      real          bklev(1000)
      real          aklay(1000)
      real          bklay(1000)
      integer       stdate(5)

      integer       writeakbk
      parameter     (writeakbk=0)

      integer       ierror
      integer       cdfid
      real          phymin(4),phymax(4)
      integer       ndim
      character*80  cfn
      real          pollon,pollatq
      real          dellon,dellat
      integer       dimid
      integer       dimids(4)
      integer       varids(4)

      goto 200

c     Get the number of levels
      nz = datar(9)                     

c     Try to open the netCDF file - create it, if necesarry
      ierror = NF90_OPEN(TRIM(cstnam),nf90_write, cdfid)
      IF ( ierror /= nf90_NoErr ) then      
         phymin(1) = real(datar(4))/1000.
         phymin(2) = max(real(datar(5))/1000.,-90.)
         phymin(3) = 1
         phymax(1) = real(datar(6))/1000.
         phymax(2) = min(real(datar(3))/1000.,90.)
         phymax(3) = nz
         ndim      = 4
         cfn       = cstnam
         call crecdf (cstnam, cdfid, phymin, phymax, ndim, cfn, 
     &               ierror ) 
         call putstart (cdfid,stdate,ierror)         
         ierror = nf90_def_dim(cdfid,'nx', datar(1), dimid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_def_dim(cdfid,'ny', datar(2), dimid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
         ierror = nf90_def_dim(cdfid,'ny', datar(3), dimid)
         IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 
      endif

c     Enter definition mode
      ierror = nf90_redef(cdfid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Write pole longitude and latitude
      pollon = real(datar(13))/1000
      if (datar(14).gt.0) then
        pollat=min(real(datar(14))/1000.,90.)
      else
        pollat=max(real(datar(14))/1000.,-90.)
      endif

      ierror = nf90_put_att(cdfid, nf90_global, 
     >                               "pollon", real(pollon))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, nf90_global, 
     >                               "pollat", real(pollat))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

C     Write grid resolution
      dellon = real(datar(7))/1000.
      dellat = real(datar(8))/1000.

      ierror = nf90_put_att(cdfid, nf90_global, 
     >                               "dellon", real(dellon))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, nf90_global, 
     >                               "dellat", real(dellat))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Skip ak,bk definition if writing ot requested
      if (writeakbk.eq.0) goto 100

c     Define dimensions for nhym, nhym
      ierror = nf90_def_dim(cdfid,'nhyi', nz, dimids(1))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_def_dim(cdfid,'nhym', nz, dimids(2))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Define variables hyai,hyam,hybi,hybm
      ierror = nf90_def_var(cdfid, 'hyai',
     >                       NF90_FLOAT,dimids(1),varids(1))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, varids(1), 
     >         "long_name", "hybrid A coefficient at layer interfaces") 
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, varids(1), 
     >          "units", "hPa") 
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

      ierror = nf90_def_var(cdfid, 'hybi',
     >                       NF90_FLOAT,dimids(1),varids(2))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, varids(2), 
     >         "long_name", "hybrid B coefficient at layer interfaces") 
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, varids(2), 
     >          "units", "1") 
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

      ierror = nf90_def_var(cdfid, 'hyam',
     >                       NF90_FLOAT,dimids(2),varids(3))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, varids(3), 
     >         "long_name", "hybrid A coefficient at layer midpoints") 
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, varids(3), 
     >          "units", "hPa") 
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

      ierror = nf90_def_var(cdfid, 'hybm',
     >                       NF90_FLOAT,dimids(2),varids(4))
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, varids(4), 
     >         "long_name", "hybrid B coefficient at layer midpoints") 
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_att(cdfid, varids(4), 
     >          "units", "1") 
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Exit definition mode
      ierror = nf90_enddef(cdfid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Write ak,bk
      ierror = nf90_put_var(cdfid, varids(1), aklay(1:nz) )
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_var(cdfid, varids(2), bklay(1:nz) )
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_var(cdfid, varids(3), aklev(1:nz) )
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)
      ierror = nf90_put_var(cdfid, varids(4), bklev(1:nz) )
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Continue
 100  continue

c     Sync
      ierror = nf90_sync(cdfid)
      IF(ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror) 

c     Close netCDF
      call clscdf(cdfid,ierror)
      IF (ierror /= nf90_NoErr) PRINT *,NF90_STRERROR(ierror)

c     Exit
 200  continue

      end