WebSVN – pvinversion.ecmwf – Path Comparison – / – /trunk/lib/ Rev 2 and /trunk/lib/ Rev 3

Ignore whitespace Rev 2 → Rev 3

 /trunk/lib/libcdfio.f
 ,0 → 1,1648
+      subroutine clscdf (cdfid, error)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This routine closes an open netCDF file.
+c     Aguments :
+c        cdfid  int  input   the id of the file to be closed.
+c        error  int  output  indicates possible errors found in this
+c                            routine.
+c                            error = 0   no errors detected.
+c                            error = 1   error detected.
+c     History:
+c        Nov. 91  PPM  UW  Created.
+c-----------------------------------------------------------------------
+      include "netcdf.inc"
+c     Argument declarations.
+      integer      cdfid, error
+c     Local variable declarations.
+      integer      ncopts
+c     Get current value of error options.
+      call ncgopt (ncopts)
+c     Make sure netCDF errors do not abort execution.
+      call ncpopt (NCVERBOS)
+c     Close requested file.
+      call ncclos (cdfid, error)
+c     Reset error options.
+      call ncpopt (ncopts)
+      end
+      subroutine cdfwopn(filnam,cdfid,ierr)
+C------------------------------------------------------------------------
+C     Opens the NetCDF file 'filnam' and returns its identifier cdfid.
+C     filnam    char    input   name of NetCDF file to open
+C     cdfid     int     output  identifier of NetCDF file
+C     ierr      int     output  error flag
+C------------------------------------------------------------------------
+      include 'netcdf.inc'
+      integer   cdfid,ierr
+      character*(*) filnam
+      integer   strend
+      call ncpopt(NCVERBOS)
+      cdfid=ncopn(filnam(1:strend(filnam)),NCWRITE,ierr)
+      return
+      end
+      subroutine cdfopn(filnam,cdfid,ierr)
+C------------------------------------------------------------------------
+C     Opens the NetCDF file 'filnam' and returns its identifier cdfid.
+C     filnam    char    input   name of NetCDF file to open
+C     cdfid     int     output  identifier of NetCDF file
+C     ierr      int     output  error flag
+C------------------------------------------------------------------------
+      include 'netcdf.inc'
+      integer   cdfid,ierr
+      character*(*) filnam
+      integer   strend
+      call ncpopt(NCVERBOS)
+      cdfid=ncopn(filnam(1:strend(filnam)),NCNOWRIT,ierr)
+      return
+      end
+      subroutine crecdf (filnam, cdfid, phymin, phymax, ndim, cfn,
+     &                   error)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This routine is called to create a netCDF file for use with
+c        the UWGAP plotting package.
+c           Any netCDF file written to must be closed with the call
+c        'call clscdf(cdfid,error)', where cdfid and error are
+c        as in the argumentlist below.
+c     Arguments:
+c        filnam  char  input   the user-supplied netCDF file name.
+c        cdfid   int   output  the file-identifier
+c        phymin  real  input   the minimum physical dimension of the
+c                              entire physical domain along each axis.
+c                              phymin is dimensioned (ndim)
+c        phymax  real  input   the maximum physical dimension of the
+c                              entire physical domain along each axis.
+c                              phymax is dimensioned (ndim)
+c        ndim    int   input   the number of dimensions in the file
+c                              (i.e. number of elements in phymin,
+c                              phymax)
+c        cfn     char  input   constants file name
+c                              ('0' = no constants file).
+c        error   int   output  indicates possible errors found in this
+c                              routine.
+c                              error = 0   no errors detected.
+c                              error = 1   error detected.
+c     History:
+c        Nov. 91  PPM  UW  Created cr3df.
+c        Jan. 92  CS   UW  Created crecdf.
+c-----------------------------------------------------------------------
+      include "netcdf.inc"
+c     Argument declarations.
+      integer        MAXDIM
+      parameter      (MAXDIM=4)
+      integer        ndim, error
+      character *(*) filnam,cfn
+      real           phymin(*), phymax(*)
+c     Local variable declarations.
+      character *(20) attnam
+      character *(1)  chrid(MAXDIM)
+      integer         cdfid, k, ibeg, iend, lenfil, ncopts
+      data            chrid/'x','y','z','a'/
+c     Get current value of error options, and make sure netCDF-errors do
+c     not abort execution
+      call ncgopt (ncopts)
+      call ncpopt(NCVERBOS)
+c     Initially set error to indicate no errors.
+      error = 0
+c     create the netCDF file
+      cdfid = nccre (trim(filnam), NCCLOB, error)
+      if (error.ne.0) go to 920
+c     define global attributes
+      do k=1,ndim
+        attnam(1:3)='dom'
+        attnam(4:4)=chrid(k)
+        attnam(5:7)='min'
+        attnam=attnam(1:7)
+        call ncapt(cdfid,NCGLOBAL,attnam,NCFLOAT,1,phymin(k),error)
+        if (error.gt.0) goto 920
+        attnam(1:3)='dom'
+        attnam(4:4)=chrid(k)
+        attnam(5:7)='max'
+        attnam=attnam(1:7)
+        call ncapt(cdfid,NCGLOBAL,attnam,NCFLOAT,1,phymax(k),error)
+        if (error.gt.0) goto 920
+      enddo
+c     define constants file name
+      if (cfn.ne.'0') then
+        call ncaptc (cdfid, NCGLOBAL, 'constants_file_name',
+c    &             NCCHAR, len_trim(cfn)+1, cfn // char(0) , error)
+     &             NCCHAR, len_trim(cfn), cfn , error)
+        if (error.gt.0) goto 920
+      endif
+c     End variable definitions.
+      call ncendf (cdfid, error)
+      if (error.gt.0) goto 920
+c     normal exit
+      call ncpopt (ncopts)
+      return
+c     error exit
+write (6, *) 'ERROR: An error occurred while attempting to ',
+     &             'create the data file in subroutine crecdf.'
+      call ncpopt (ncopts)
+      call ncclos (cdfid, error)
+      end
+      subroutine opncdf(filnam, cdfid,
+     &                 phymin, phymax, ndim, varnam, nvar, cfn, error)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This routine is called to open a netCDF file for read and write
+c        with the UWGAP plotting package.
+c     Arguments:
+c        filnam  char  input   the user-supplied netCDF file name.
+c        cdfid   int   output  the file-identifier
+c        phymin  real  output  the minimum physical dimension of the
+c                              entire physical domain along each axis.
+c                              phymin is dimensioned (ndim)
+c        phymax  real  output  the maximum physical dimension of the
+c                              entire physical domain along each axis.
+c                              phymax is dimensioned (ndim)
+c        ndim    int   output  the number of dimensions in the file
+c                              (i.e. number of elements in phymin,
+c                              phymax)
+c        varnam  char  output  an array containing the variable names.
+c                              varnam is dimensioned (nvar).
+c        nvar    int   output  the number of variables in the file
+c        cfn     char  output  constants file name
+c                              ('0'=no constants file).
+c        error   int   output  indicates possible errors found in this
+c                              routine.
+c                              error = 0   no errors detected.
+c                              error = 1   error detected.
+c     History:
+c        Nov. 91  PPM  UW  Created cr3df.
+c        Jan. 92  CS   UW  Created opncdf.
+c-----------------------------------------------------------------------
+      include "netcdf.inc"
+c     Argument declarations.
+      integer        MAXDIM
+      parameter      (MAXDIM=4)
+      integer        ndim, nvar, error
+      character *(*) filnam, varnam(*),cfn
+      real           phymin(*), phymax(*)
+c     Local variable declarations.
+      character *(20) attnam,vnam
+      character *(1)  chrid(MAXDIM)
+      integer         cdfid, i,k
+      integer         ncopts, ndims,ngatts,recdim
+      integer         nvdims,vartyp,nvatts,vardim(MAXDIM)
+      real            attval
+      integer         lenstr
+      data            chrid/'x','y','z','a'/
+      data            lenstr/80/
+c     Get current value of error options and make sure netCDF-errors do
+c     not abort execution
+      call ncgopt (ncopts)
+      call ncpopt(NCVERBOS)
+c     Initially set error to indicate no errors.
+      error = 0
+c     open the netCDF file for write
+      cdfid = ncopn (trim(filnam), NCWRITE, error)
+      if (error.ne.0) then
+c       try to open the netCDF file for read
+        cdfid = ncopn (trim(filnam), NCNOWRIT, error)
+        if (error.ne.0) go to 920
+      endif
+c     inquire for number of variables
+      call ncinq(cdfid,ndims,nvar,ngatts,recdim,error)
+      if (error.eq.1) goto 920
+c     read the variables
+      do i=1,nvar
+        call ncvinq(cdfid,i,varnam(i),vartyp,nvdims,vardim,
+     &                         nvatts,error)
+        if (vartyp.ne.NCFLOAT) error=1
+        if (error.gt.0) goto 920
+      enddo
+c     get global attributes
+      k=0
+continue
+        k=k+1
+        attnam(1:3)='dom'
+        attnam(4:4)=chrid(k)
+        attnam(5:7)='min'
+        attnam=attnam(1:7)
+c       switch off error message
+        call ncpopt(0)
+c       check whether dimension k is present
+        call ncagt(cdfid,NCGLOBAL,attnam,attval,error)
+        if (error.gt.0) goto 110
+        phymin(k)=attval
+        attnam(1:3)='dom'
+        attnam(4:4)=chrid(k)
+        attnam(5:7)='max'
+        attnam=attnam(1:7)
+        call ncagt(cdfid,NCGLOBAL,attnam,attval,error)
+        if (error.gt.0) goto 920
+        phymax(k)=attval
+      if (k.lt.3) goto 100
+      k=k+1
+c     define ndim-parameter
+continue
+      ndim=k-1
+      error=0
+c     switch on error messages
+      call ncpopt(NCVERBOS)
+c     get constants file name
+c      call ncagt(cdfid,NCGLOBAL,'constants_file_name',cfn,error)
+c     ! chrigel
+      call ncagtc(cdfid,NCGLOBAL,'constants_file_name',cfn,lenstr,error)
+      if (error.gt.0) cfn='0'
+c     normal exit
+      call ncpopt (ncopts)
+      return
+c     error exit
+write (6, *) 'ERROR: An error occurred while attempting to ',
+     &             'read the data file in subroutine opncdf.'
+      call ncclos (cdfid, error)
+      call ncpopt (ncopts)
+      end
+      subroutine readcdf(filnam, cdfid,
+     &                 phymin, phymax, ndim, varnam, nvar, cfn, error)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This routine is called to open a netCDF file for read
+c        with the UWGAP plotting package.
+c     Arguments:
+c        filnam  char  input   the user-supplied netCDF file name.
+c        cdfid   int   output  the file-identifier
+c        phymin  real  output  the minimum physical dimension of the
+c                              entire physical domain along each axis.
+c                              phymin is dimensioned (ndim)
+c        phymax  real  output  the maximum physical dimension of the
+c                              entire physical domain along each axis.
+c                              phymax is dimensioned (ndim)
+c        ndim    int   output  the number of dimensions in the file
+c                              (i.e. number of elements in phymin,
+c                              phymax)
+c        varnam  char  output  an array containing the variable names.
+c                              varnam is dimensioned (nvar).
+c        nvar    int   output  the number of variables in the file
+c        cfn     char  output  constants file name
+c                              ('0'=no constants file).
+c        error   int   output  indicates possible errors found in this
+c                              routine.
+c                              error = 0   no errors detected.
+c                              error = 1   error detected.
+c     History:
+c        Nov. 91  PPM  UW  Created cr3df.
+c        Jan. 92  CS   UW  Created opncdf.
+c-----------------------------------------------------------------------
+      include "netcdf.inc"
+c     Argument declarations.
+      integer        MAXDIM
+      parameter      (MAXDIM=4)
+      integer        ndim, nvar, error
+      character *(*) filnam, varnam(*),cfn
+      real           phymin(*), phymax(*)
+c     Local variable declarations.
+      character *(20) attnam
+      character *(1)  chrid(MAXDIM)
+      integer         cdfid, i,k
+      integer         ncopts, ndims,ngatts,recdim
+      integer         nvdims,vartyp,nvatts,vardim(MAXDIM)
+      real            attval
+      integer         lenstr
+      data            chrid/'x','y','z','a'/
+      data            lenstr/80/
+c     Get current value of error options.
+      call ncgopt (ncopts)
+c     make sure netCDF-errors do not abort execution
+      call ncpopt(NCVERBOS)
+c     Initially set error to indicate no errors.
+      error = 0
+c     open the netCDF file for read
+      cdfid = ncopn (trim(filnam), NCNOWRIT, error)
+      if (error.ne.0) go to 920
+c     inquire for number of variables
+      call ncinq(cdfid,ndims,nvar,ngatts,recdim,error)
+      if (error.eq.1) goto 920
+c     read the variables
+      do i=1,nvar
+        call ncvinq(cdfid,i,varnam(i),vartyp,nvdims,vardim,
+     &                         nvatts,error)
+        if (vartyp.ne.NCFLOAT) error=1
+c       print *,varnam(i),nvdims,nvatts
+        if (error.gt.0) goto 920
+      enddo
+c     get global attributes
+      k=0
+continue
+        k=k+1
+        attnam(1:3)='dom'
+        attnam(4:4)=chrid(k)
+        attnam(5:7)='min'
+        attnam=attnam(1:7)
+c       switch off error message
+        call ncpopt(0)
+c       check whether dimension k is present
+        call ncagt(cdfid,NCGLOBAL,attnam,attval,error)
+        if (error.gt.0) goto 110
+        phymin(k)=attval
+        attnam(1:3)='dom'
+        attnam(4:4)=chrid(k)
+        attnam(5:7)='max'
+        attnam=attnam(1:7)
+        call ncagt(cdfid,NCGLOBAL,attnam,attval,error)
+        if (error.gt.0) goto 920
+        phymax(k)=attval
+      if (k.lt.4) goto 100
+      k=k+1
+c     define ndim-parameter
+continue
+      ndim=k-1
+      error=0
+c     switch on error messages
+      call ncpopt(NCVERBOS)
+c     get constants file name
+c      call ncagt(cdfid,NCGLOBAL,'constants_file_name',cfn,error)
+c     ! chrigel
+      call ncagtc(cdfid,NCGLOBAL,'constants_file_name',cfn,lenstr,error)
+      if (error.gt.0) cfn='0'
+c     print *,cfn
+c     normal exit
+      call ncpopt (ncopts)
+      return
+c     error exit
+write (6, *) 'ERROR: An error occurred while attempting to ',
+     &             'read the data file in subroutine opncdf.'
+      call ncclos (cdfid, error)
+      call ncpopt (ncopts)
+      end
+      subroutine getcdf (cdfid, varnam, ndim, misdat,
+     &                       vardim, varmin, varmax, stag, dat, error)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This routine is called to get a variable and its attributes
+c        from a netCDF file for use with the UWGAP plotting package.
+c        It is assumed that the data is floating-point data. Prior to
+c        calling this routine, the file must be opened with a call to
+c        opncdf.
+c     Arguments:
+c        cdfid   int   input   file-identifier
+c                              (can be obtained by calling routine
+c                              opncdf)
+c        varnam  char  input   the user-supplied variable name.
+c                              (can be obtained by calling routine
+c                              opncdf)
+c        ndim    int   output  the number of dimensions (ndim<=4)
+c        misdat  real  output  missing data value for the variable.
+c        vardim  int   output  the dimensions of the variable.
+c                              is dimensioned at least (ndim).
+c        varmin  real  output  the location in physical space of the
+c                              origin of each variable.
+c                              is dimensioned at least Min(ndim,3).
+c        varmax  real  output  the extent of each variable in physical
+c                              space.
+c                              is dimensioned at least Min(ndim,3).
+c        stag    real  output  the grid staggering for each variable.
+c                              is dimensioned at least Min(ndim,3).
+c        dat     real  output  data-array dimensioned suffiecently
+c                              large, at least
+c                              vardim(1)* ... vardim(ndim)
+c        error   int   output  indicates possible errors found in this
+c                              routine.
+c                              error = 0   no errors detected.
+c                              error = 1   error detected.
+c     History:
+c        Nov. 91  PPM  UW  Created cr3df.
+c        Jan. 92  CS   UW  Created getcdf.
+c-----------------------------------------------------------------------
+      include "netcdf.inc"
+c     Argument declarations.
+      integer        MAXDIM
+      parameter      (MAXDIM=4)
+      character *(*) varnam
+      integer        vardim(*), ndim, error, cdfid
+      real           misdat,  stag(*), varmin(*), varmax(*), dat(*)
+c     Local variable declarations.
+      character *(20) dimnam(100),attnam
+      character *(1)  chrid(MAXDIM)
+      integer         id,i,k,corner(MAXDIM)
+      integer         ndims,nvars,ngatts,recdim,dimsiz(100)
+      integer         vartyp,nvatts, ncopts
+      data            chrid/'x','y','z','a'/
+      data            corner/1,1,1,1/
+c     Get current value of error options, and make sure netCDF-errors do
+c     not abort execution
+      call ncgopt (ncopts)
+      call ncpopt(NCVERBOS)
+c     Initially set error to indicate no errors.
+      error = 0
+c     inquire for number of dimensions
+      call ncinq(cdfid,ndims,nvars,ngatts,recdim,error)
+      if (error.eq.1) goto 920
+c     read dimension-table
+      do i=1,ndims
+        call ncdinq(cdfid,i,dimnam(i),dimsiz(i),error)
+        if (error.gt.0) goto 920
+      enddo
+c     get id of the variable
+      id=ncvid(cdfid,varnam,error)
+      if (error.eq.1) goto 910
+c     inquire about variable
+      call ncvinq(cdfid,id,varnam,vartyp,ndim,vardim,nvatts,error)
+      if (vartyp.ne.NCFLOAT) error=1
+      if (error.gt.0) goto 920
+c     Make sure ndim <= MAXDIM.
+      if (ndim.gt.MAXDIM) then
+         error = 1
+         go to 900
+      endif
+c     get dimensions from dimension-table
+      do k=1,ndim
+        vardim(k)=dimsiz(vardim(k))
+      enddo
+c     get attributes
+      do k=1,min0(ndim,3)
+c       get staggering
+        attnam(1:1)=chrid(k)
+        attnam(2:5)='stag'
+        attnam=attnam(1:5)
+        call ncagt(cdfid,id,attnam,stag(k),error)
+        if (error.gt.0) goto 920
+c       get min postion
+        attnam(1:1)=chrid(k)
+        attnam(2:4)='min'
+        attnam=attnam(1:4)
+        call ncagt(cdfid,id,attnam,varmin(k),error)
+        if (error.gt.0) goto 920
+c       get max position
+        attnam(1:1)=chrid(k)
+        attnam(2:4)='max'
+        attnam=attnam(1:4)
+        call ncagt(cdfid,id,attnam,varmax(k),error)
+        if (error.gt.0) goto 920
+      enddo
+c     get missing data value
+      call ncagt(cdfid,id,'missing_data',misdat,error)
+      if (error.gt.0) goto 920
+c     get data
+      call ncvgt(cdfid,id,corner,vardim,dat,error)
+      if (error.gt.0) goto 920
+c     normal exit
+      call ncpopt (ncopts)
+      return
+c     Error exits.
+write (6, *) 'ERROR: When calling getcdf, the number of ',
+     &             'variable dimensions must be less or equal 4.'
+      call ncpopt (ncopts)
+      call ncclos (cdfid, error)
+      return
+write (6, *) 'ERROR: The selected variable could not be found ',
+     &             'in the file by getcdf.'
+      call ncpopt (ncopts)
+      call ncclos (cdfid, error)
+      return
+write (6, *) 'ERROR: An error occurred while attempting to ',
+     &             'read the data file in subroutine getcdf.'
+      call ncpopt (ncopts)
+      call ncclos (cdfid, error)
+      return
+      end
+      subroutine putcdf (cdfid, varnam, ndim, misdat,
+     &                       vardim, varmin, varmax, stag, dat, error)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This routine is called to put a variable and its attributes
+c        onto a netCDF file for use with the UWGAP plotting package.
+c        It is assumed that the data is floating-point data. Prior to
+c        calling this routine, the file must be created (crecdf) or
+c        opened (opncdf).
+c           Any netCDF file written to must be closed with the call
+c        call ncclos(cdfid,error), where cdfid and error are
+c        as in the argumentlist below.
+c     Arguments:
+c        cdfid   int   input   file-identifier
+c                              (can be obtained by calling routine
+c                              opncdf)
+c        varnam  char  input   the user-supplied variable name.
+c                              (can be obtained by calling routine
+c                              opncdf)
+c        ndim    int   input   the number of dimensions (ndim<=4)
+c        misdat  real  input   missing data value for the variable.
+c        vardim  int   input   the dimensions of the variable.
+c                              is dimensioned at least (ndim).
+c        varmin  real  input   the location in physical space of the
+c                              origin of each variable.
+c                              is dimensioned at least Min(ndim,3).
+c        varmax  real  input   the extent of each variable in physical
+c                              space.
+c                              is dimensioned at least Min(ndim,3).
+c        stag    real  input   the grid staggering for each variable.
+c                              is dimensioned at least Min(ndim,3).
+c        dat     real  input   data-array dimensioned suffiecently
+c                              large, at least
+c                              vardim(1)* ... vardim(ndim)
+c        error   int   output  indicates possible errors found in this
+c                              routine.
+c                              error = 0   no errors detected.
+c                              error = 1   error detected.
+c     History:
+c        Nov. 91  PPM  UW  Created cr3df, wr3df.
+c        Jan. 92  CS   UW  Created putcdf.
+c-----------------------------------------------------------------------
+      include "netcdf.inc"
+c     Argument declarations.
+      integer        MAXDIM
+      parameter      (MAXDIM=4)
+      character *(*) varnam
+      integer        vardim(*), ndim, error, cdfid
+      real           misdat,  stag(*), varmin(*), varmax(*), dat(*)
+c     Local variable declarations.
+      character *(20) dimnam,attnam,dimchk
+      character *(1)  chrid(MAXDIM)
+      character *(20) dimnams(MAXNCDIM)
+      integer         dimvals(MAXNCDIM)
+      integer         numdims,numvars,numgats,dimulim
+      integer         id,did(MAXDIM),i,k,corner(MAXDIM)
+      integer         ncopts
+      integer         ibeg,iend
+      data            chrid/'x','y','z','t'/
+      data            corner/1,1,1,1/
+c     Get current value of error options.
+      call ncgopt (ncopts)
+c     make sure netCDF-errors do not abort execution
+      call ncpopt(NCVERBOS)
+c     Initially set error to indicate no errors.
+      error = 0
+c     Make sure ndim <= MAXDIM.
+      if (ndim.gt.MAXDIM) then
+         error = 1
+         go to 900
+      endif
+c     Read existing dimensions-declarations from the file
+      call ncinq(cdfid,numdims,numvars,numgats,dimulim,error)
+      if (error.ne.0) numdims=0
+      if (numdims.gt.0) then
+        do i=1,numdims
+          call ncdinq(cdfid,i,dimnams(i),dimvals(i),error)
+c          print *,dimnams(i),dimvals(i)
+        enddo
+      endif
+c     put file into define mode
+      call ncredf(cdfid,error)
+      if (error.ne.0) goto 920
+c     define the dimension
+      do k=1,ndim
+c       define the dimension-name
+        dimnam(1:3)='dim'
+        dimnam(4:4)=chrid(k)
+        dimnam(5:5)='_'
+        dimnam(6:5+len_trim(varnam))=trim(varnam)
+        dimnam=dimnam(1:5+len_trim(varnam))
+        did(k)=-1
+        if (numdims.gt.0) then
+c         check if an existing dimension-declaration can be used
+c         instead of defining a nuw dimension
+          do i=1,numdims
+            dimchk=dimnams(i)
+            if ((vardim(k).eq.dimvals(i)).and.
+     &        (dimnam(1:4).eq.dimchk(1:4))) then
+              did(k)=i
+              goto 100
+            endif
+          enddo
+continue
+        endif
+        if (did(k).lt.0) then
+c         define the dimension
+          did(k)=ncddef(cdfid,dimnam,vardim(k),error)
+          if (error.ne.0) goto 920
+        endif
+      enddo
+c     define variable
+      id=ncvdef(cdfid,varnam,NCFLOAT,ndim,did,error)
+      if (error.ne.0) goto 920
+c     define attributes
+      do k=1,min0(ndim,3)
+c       staggering
+        attnam(1:1)=chrid(k)
+        attnam(2:5)='stag'
+        attnam=attnam(1:5)
+        call ncapt(cdfid,id,attnam,NCFLOAT,1,stag(k),error)
+        if (error.gt.0) goto 920
+c       min postion
+        attnam(1:1)=chrid(k)
+        attnam(2:4)='min'
+        attnam=attnam(1:4)
+        call ncapt(cdfid,id,attnam,NCFLOAT,1,varmin(k),error)
+        if (error.gt.0) goto 920
+c       max position
+        attnam(1:1)=chrid(k)
+        attnam(2:4)='max'
+        attnam=attnam(1:4)
+        call ncapt(cdfid,id,attnam,NCFLOAT,1,varmax(k),error)
+        if (error.gt.0) goto 920
+      enddo
+c     define missing data value
+      call ncapt(cdfid,id,'missing_data',NCFLOAT,1,misdat,error)
+      if (error.gt.0) goto 920
+c     leave define mode
+      call ncendf(cdfid,error)
+      if (error.gt.0) goto 920
+c     define data
+      call ncvpt(cdfid,id,corner,vardim,dat,error)
+      if (error.gt.0) goto 920
+c     synchronyse output to disk and exit
+      call ncsnc (cdfid,error)
+      call ncpopt (ncopts)
+      return
+c     Error exits.
+write (6, *) 'ERROR: When calling putcdf, the number of ',
+     &             'variable dimensions must be less or equal 4.'
+      call ncpopt (ncopts)
+      call ncclos (cdfid, error)
+      return
+write (6, *) 'ERROR: An error occurred while attempting to ',
+     &             'write the data file in subroutine putcdf.'
+      call ncpopt (ncopts)
+      call ncclos (cdfid, error)
+      return
+      end
+c
+c
+      subroutine getdef (cdfid, varnam, ndim, misdat,
+     &                              vardim, varmin, varmax, stag, error)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This routine is called to get the dimensions and attributes of
+c        a variable from an IVE-NetCDF file for use with the IVE plotting
+c        package. Prior to calling this routine, the file must be opened
+c        with a call to opncdf.
+c     Arguments:
+c        cdfid   int   input   file-identifier
+c                              (can be obtained by calling routine
+c                              opncdf)
+c        varnam  char  input   the user-supplied variable name.
+c                              (can be obtained by calling routine
+c                              opncdf)
+c        ndim    int   output  the number of dimensions (ndim<=4)
+c        misdat  real  output  missing data value for the variable.
+c        vardim  int   output  the dimensions of the variable.
+c                              Is dimensioned at least (ndim).
+c        varmin  real  output  the location in physical space of the
+c                              origin of each variable.
+c                              Is dimensioned at least Min(3,ndim).
+c        varmax  real  output  the extend of each variable in physical
+c                              space.
+c                              Is dimensioned at least Min(3,ndim).
+c        stag    real  output  the grid staggering for each variable.
+c                              Is dimensioned at least Min(3,ndim).
+c        error   int   output  indicates possible errors found in this
+c                              routine.
+c                              error = 0   no errors detected.
+c                              error = 1   the variable is not on the file.
+c                              error =10   other errors.
+c     History:
+c       Apr. 93    Christoph Schaer (ETHZ)     Created.
+c-----------------------------------------------------------------------
+      include "netcdf.inc"
+c     Argument declarations.
+      integer        MAXDIM
+      parameter      (MAXDIM=4)
+      character *(*) varnam
+      integer        vardim(*), ndim, error, cdfid
+      real           misdat,  stag(*), varmin(*), varmax(*)
+c     Local variable declarations.
+      character *(20) dimnam(MAXNCDIM),attnam,vnam
+      character *(1)  chrid(MAXDIM)
+      integer         id,i,k
+      integer         ndims,nvars,ngatts,recdim,dimsiz(MAXNCDIM)
+      integer         vartyp,nvatts, ncopts
+      data            chrid/'x','y','z','t'/
+c     Get current value of error options.
+      call ncgopt (ncopts)
+c     make sure NetCDF-errors do not abort execution
+      call ncpopt(NCVERBOS)
+c     Initially set error to indicate no errors.
+      error = 0
+c     inquire for number of dimensions
+      call ncinq(cdfid,ndims,nvars,ngatts,recdim,error)
+      if (error.eq.1) goto 920
+c     read dimension-table
+      do i=1,ndims
+        call ncdinq(cdfid,i,dimnam(i),dimsiz(i),error)
+        if (error.gt.0) goto 920
+      enddo
+c     get id of the variable
+      id=ncvid(cdfid,varnam,error)
+      if (error.eq.1) goto 910
+c     inquire about variable
+      call ncvinq(cdfid,id,vnam,vartyp,ndim,vardim,nvatts,error)
+      if (vartyp.ne.NCFLOAT) error=1
+      if (error.gt.0) goto 920
+c     Make sure ndim <= MAXDIM.
+      if (ndim.gt.MAXDIM) then
+         error = 1
+         go to 900
+      endif
+c     get dimensions from dimension-table
+      do k=1,ndim
+        vardim(k)=dimsiz(vardim(k))
+      enddo
+c     get attributes
+      do k=1,min0(3,ndim)
+c       get min postion
+        attnam(1:1)=chrid(k)
+        attnam(2:4)='min'
+        attnam=attnam(1:4)
+        call ncagt(cdfid,id,attnam,varmin(k),error)
+        if (error.gt.0) goto 920
+c       get max position
+        attnam(1:1)=chrid(k)
+        attnam(2:4)='max'
+        attnam=attnam(1:4)
+        call ncagt(cdfid,id,attnam,varmax(k),error)
+        if (error.gt.0) goto 920
+c       get staggering
+        attnam(1:1)=chrid(k)
+        attnam(2:5)='stag'
+        attnam=attnam(1:5)
+        call ncagt(cdfid,id,attnam,stag(k),error)
+        if (error.gt.0) goto 920
+      enddo
+c     get missing data value
+      call ncagt(cdfid,id,'missing_data',misdat,error)
+      if (error.gt.0) goto 920
+c     normal exit
+      call ncpopt (ncopts)
+      return
+c     Error exits.
+write (6, *) '*ERROR*: When calling getcdf, the number of ',
+     &             'variable dimensions must be less or equal 4.'
+      call ncpopt (ncopts)
+      call ncclos (cdfid, error)
+      return
+write (6, *) '*ERROR*: The selected variable could not be found ',
+     &             'in the file by getcdf.'
+      call ncpopt (ncopts)
+      call ncclos (cdfid, error)
+      return
+write (6, *) '*ERROR*: An error occurred while attempting to ',
+     &             'read the data file in subroutine getcdf.'
+      call ncpopt (ncopts)
+      call ncclos (cdfid, error)
+      end
+      subroutine getdat(cdfid, varnam, time, level, dat, error)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This routine is called to read the data of a variable
+c        from an IVE-NetCDF file for use with the IVE plotting package.
+c        Prior to calling this routine, the file must be opened with
+c        a call to opncdf (for extension) or crecdf (for creation) or
+c        readcdf (for readonly).
+c     Arguments:
+c        cdfid   int   input   file-identifier
+c                              (must be obtained by calling routine
+c                              opncdf,readcdf  or crecdf)
+c        varnam  char  input   the user-supplied variable name (must
+c                              previously be defined with a call to
+c                              putdef)
+c        time    real  input   the user-supplied time-level of the
+c                              data to be read from the file (the time-
+c                              levels stored in the file can be obtained
+c                              with a call to gettimes).
+c        level   int   input   the horizontal level(s) to be read
+c                              to the NetCDF file. Suppose that the
+c                              variable is defined as (nx,ny,nz,nt).
+c                              level>0: the call reads the subdomain
+c                                       (1:nx,1:ny,level,itimes)
+c                              level=0: the call reads the subdomain
+c                                       (1:nx,1:ny,1:nz,itimes)
+c                              Here itimes is the time-index corresponding
+c                              to the value of 'time'.
+c        dat     real  output  data-array dimensioned sufficiently
+c                              large. The dimensions (nx,ny,nz)
+c                              of the variable must previously be defined
+c                              with a call to putdef. No previous
+c                              definition of the time-dimension is
+c                              required.
+c        error   int   output  indicates possible errors found in this
+c                              routine.
+c                              error = 0   no errors detected.
+c                              error = 1   the variable is not present on
+c                                          the file.
+c                              error = 2   the value of 'time' is not
+c                                          known.to the file.
+c                              error = 3   inconsistent value of level
+c                              error =10   another error.
+c     History:
+c       March 93    Heini Wernli (ETHZ)      Created wr2cdf.
+c       April 93    Bettina Messmer (ETHZ)   Created putdat.
+c       June  93    Christoph Schaer (ETHZ)  Created getdat
+c-----------------------------------------------------------------------
+      include "netcdf.inc"
+C     Declaration of local variables
+      character*(*) varnam
+      character*(20) chars
+      integer cdfid
+      real        dat(*)
+      real        misdat,varmin(4),varmax(4),stag(4)
+      real        time, timeval
+      integer     corner(4),edgeln(4),didtim,vardim(4),ndims
+      integer     error, ierr
+      integer     level,ntime
+      integer     idtime,idvar,iflag
+      integer     i
+      call ncpopt(NCVERBOS)
+c     access the variable
+      call getdef (cdfid, trim(varnam), ndims, misdat,
+     &                           vardim, varmin, varmax, stag, ierr)
+      if (ierr.ne.0) then
+        print *,'*ERROR* in getdef in getdat'
+        error=1
+        return
+      endif
+      idvar=ncvid(cdfid,trim(varnam),ierr)
+      if (ierr.ne.0) then
+        print *,'*ERROR* in ncvid in getdat'
+        error=1
+        return
+      endif
+C     Get times-array
+      didtim=ncdid(cdfid,'time',ierr)
+      if (ierr.ne.0) then
+        print *,'*ERROR* didtim in getdat'
+        error=10
+        return
+      endif
+      call ncdinq(cdfid,didtim,chars,ntime,ierr)
+      if (ierr.ne.0) then
+        print *,'*ERROR* in ncdinq in getdat'
+        error=10
+        return
+      endif
+      idtime=ncvid(cdfid,'time',ierr)
+      if (ierr.ne.0) then
+        print *,'*ERROR* in ncvid for time in getdat'
+        error=10
+        return
+      endif
+c     find appropriate time-index
+      iflag=0
+      do i=1,ntime
+        call ncvgt1(cdfid,idtime,i,timeval,ierr)
+        if (ierr.ne.0) print *,'*ERROR* in ncvgt1 in getdat'
+        if (time.eq.timeval) iflag=i
+      enddo
+      if (iflag.eq.0) then
+        error=2
+        print *,'Error: Unknown time in getdat'
+        return
+      endif
+C     Define data volume to be written (index space)
+      corner(1)=1
+      corner(2)=1
+      edgeln(1)=vardim(1)
+      edgeln(2)=vardim(2)
+      if (level.eq.0) then
+        corner(3)=1
+        edgeln(3)=vardim(3)
+      else if ((level.le.vardim(3)).and.(level.ge.1)) then
+        corner(3)=level
+        edgeln(3)=1
+      else
+        error=3
+        return
+      endif
+      corner(4)=iflag
+      edgeln(4)=1
+C     Read data from NetCDF file
+c      print *,'getdat vor Aufruf ncvgt'
+c      print *,'cdfid ',cdfid
+c      print *,'idvar ',idvar
+c      print *,'corner ',corner
+c      print *,'edgeln ',edgeln
+      call ncvgt(cdfid,idvar,corner,edgeln,dat,error)
+      if (error.ne.0) then
+        print *, '*ERROR* in ncvgt in getdat'
+        error=10
+      endif
+      end
+      subroutine putdat(cdfid, varnam, time, level, dat, error)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This routine is called to write the data of a variable
+c        to an IVE-NetCDF file for use with the IVE plotting package.
+c        Prior to calling this routine, the file must be opened with
+c        a call to opncdf (for extension) or crecdf (for creation), the
+c        variable must be defined with a call to putdef.
+c     Arguments:
+c        cdfid   int   input   file-identifier
+c                              (must be obtained by calling routine
+c                              opncdf or crecdf)
+c        varnam  char  input   the user-supplied variable name (must
+c                              previously be defined with a call to
+c                              putdef)
+c        time    real  input   the user-supplied time-level of the
+c                              data to be written to the file (the time-
+c                              levels stored in the file can be obtained
+c                              with a call to gettimes). If 'time' is not
+c                              yet known to the file, a knew time-level is
+c                              allocated and appended to the times-array.
+c        level   int input     the horizontal level(s) to be written
+c                              to the NetCDF file. Suppose that the
+c                              variable is defined as (nx,ny,nz,nt).
+c                              level>0: the call writes the subdomain
+c                                       (1:nx,1:ny,level,itimes)
+c                              level=0: the call writes the subdomain
+c                                       (1:nx,1:ny,1:nz,itimes)
+c                              Here itimes is the time-index corresponding
+c                              to the value of 'time'.
+c        dat     real  output  data-array dimensioned sufficiently
+c                              large. The dimensions (nx,ny,nz)
+c                              of the variable must previously be defined
+c                              with a call to putdef. No previous
+c                              definition of the time-dimension is
+c                              required.
+c        error   int output    indicates possible errors found in this
+c                              routine.
+c                              error = 0   no errors detected.
+c                              error = 1   the variable is not present on
+c                                          the file.
+c                              error = 2   the value of 'time' is new, but
+c                                          appending it would yield a non
+c                                          ascending times-array.
+c                              error = 3   inconsistent value of level
+c                              error =10   another error.
+c     History:
+c       March 93    Heini Wernli (ETHZ)      Created wr2cdf.
+c       April 93    Bettina Messmer (ETHZ)   Created putdat.
+c-----------------------------------------------------------------------
+      include "netcdf.inc"
+C     Declaration of local variables
+      character*(*) varnam
+      character*(20) chars
+      integer cdfid
+      real      dat(*)
+      real      misdat,varmin(4),varmax(4),stag(4)
+      real      time, timeval
+      data      stag/0.,0.,0.,0./
+      integer   corner(4),edgeln(4),did(4),vardim(4),ndims
+      integer   error, ierr
+      integer   level,ntime
+      integer   idtime,idvar,iflag
+      integer   i
+      call ncpopt(NCVERBOS)
+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 variable
+      idvar=ncvid(cdfid,trim(varnam),ierr)
+      if (ierr.ne.0) print *,'*ERROR* in ncvid in putdat'
+c     get times-array
+      did(4)=ncdid(cdfid,'time',ierr)
+      if (ierr.ne.0) print *,'*ERROR* did(4) in putdat'
+      call ncdinq(cdfid,did(4),chars,ntime,ierr)
+      if (ierr.ne.0) print *,'*ERROR* in ncdinq in putdat'
+      idtime=ncvid(cdfid,'time',ierr)
+      if (ierr.ne.0) print *,'*ERROR* in ncvid in putdat'
+C     Check if a new time step is starting
+      iflag=0
+      do i=1,ntime
+        call ncvgt1(cdfid,idtime,i,timeval,ierr)
+        if (ierr.ne.0) print *,'*ERROR* in ncvgt1 in putdat'
+        if (time.eq.timeval) iflag=i
+      enddo
+      if (iflag.eq.0) then              ! new time step
+        ntime=ntime+1
+        iflag=ntime
+        idtime=ncvid(cdfid,'time',ierr)
+        if (ierr.ne.0) print *, '*ERROR* in ncvid in putdat'
+        call ncvpt1(cdfid,idtime,ntime,time,ierr)
+        if (ierr.ne.0) print *, '*ERROR* in ncvpt1 in putdat'
+      endif
+C     Define data volume to write on the NetCDF file in index space
+      corner(1)=1               ! starting corner of data volume
+      corner(2)=1
+      edgeln(1)=vardim(1)       ! edge lengthes of data volume
+      edgeln(2)=vardim(2)
+      if (level.eq.0) then
+        corner(3)=1
+        edgeln(3)=vardim(3)
+      else
+        corner(3)=level
+        edgeln(3)=1
+      endif
+      corner(4)=iflag
+      edgeln(4)=1
+C     Put data on NetCDF file
+c      print *,'vor Aufruf ncvpt d.h. Daten schreiben in putdat '
+c      print *,'cdfid ',cdfid
+c      print *,'idvar ',idvar
+c      print *,'corner ',corner
+c      print *,'edgeln ',edgeln
+      call ncvpt(cdfid,idvar,corner,edgeln,dat,error)
+      if (error.ne.0) then
+        print *, '*ERROR* in ncvpt in putdat - Put data on NetCDF file'
+      endif
+C     Synchronize output to disk and close the files
+      call ncsnc(cdfid,ierr)
+      if (ierr.ne.0) print *, '*ERROR* in ncsnc in putdat'
+      end
+      subroutine putdef (cdfid, varnam, ndim, misdat,
+     &                            vardim, varmin, varmax, stag, error)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This routine is called to define the dimensions and the
+c        attributes of a variable on an IVE-NetCDF file for use with the
+c        IVE plotting package. Prior to calling this routine, the file must
+c        be opened with a call to opncdf (extend an existing file) or
+c        crecdf (create a new file).
+c     Arguments:
+c        cdfid   int   input   file-identifier
+c                              (can be obtained by calling routine
+c                              opncdf)
+c        varnam  char  input   the user-supplied variable name.
+c        ndim    int   input   the number of dimensions (ndim<=4).
+c                              Upon ndim=4, the fourth dimension of the
+c                              variable is specified as 'unlimited'
+c                              on the file (time-dimension). It can
+c                              later be extended to arbitrary length.
+c        misdat  real  input   missing data value for the variable.
+c        vardim  int   input   the dimensions of the variable.
+c                              Is dimensioned at least Min(3,ndim).
+c        varmin  real  input   the location in physical space of the
+c                              origin of each variable.
+c                              Is dimensioned at least Min(3,ndim).
+c        varmax  real  input   the extent of each variable in physical
+c                              space.
+c                              Is dimensioned at least Min(ndim).
+c        stag    real  input   the grid staggering for each variable.
+c                              Is dimensioned at least Min(3,ndim).
+c        error   int   output  indicates possible errors found in this
+c                              routine.
+c                              error = 0   no errors detected.
+c                              error =10   other errors detected.
+c     History:
+c       Apr. 93    Christoph Schaer (ETHZ)     Created.
+c-----------------------------------------------------------------------
+      include "netcdf.inc"
+c     Argument declarations.
+      integer        MAXDIM
+      parameter      (MAXDIM=4)
+      character *(*) varnam
+      integer        vardim(*), ndim, error, cdfid
+      real           misdat,  stag(*), varmin(*), varmax(*)
+c     Local variable declarations.
+      character *(20) dimnam,attnam,dimchk
+      character *(1)  chrid(MAXDIM)
+      character *(20) dimnams(MAXNCDIM)
+      integer         dimvals(MAXNCDIM)
+      integer         numdims,numvars,numgats,dimulim
+      integer         id,did(MAXDIM),idtime,i,k,ierr
+      integer         ncopts
+      integer         ibeg,iend
+      data            chrid/'x','y','z','t'/
+c     Get current value of error options.
+      call ncgopt (ncopts)
+c     make sure NetCDF-errors do not abort execution
+      call ncpopt(NCVERBOS)
+c     Initially set error to indicate no errors.
+      error = 0
+c     Make sure ndim <= MAXDIM.
+      if (ndim.gt.MAXDIM) then
+         error = 10
+         go to 900
+      endif
+c     Read existing dimensions-declarations from the file
+      call ncinq(cdfid,numdims,numvars,numgats,dimulim,error)
+      if (numdims.gt.0) then
+        do i=1,numdims
+          call ncdinq(cdfid,i,dimnams(i),dimvals(i),error)
+c         print *,dimnams(i),dimvals(i)
+        enddo
+      endif
+c     put file into define mode
+      call ncredf(cdfid,error)
+      if (error.ne.0) goto 920
+c     define spatial dimensions
+      do k=1,min0(3,ndim)
+c       define the default dimension-name
+        dimnam(1:3)='dim'
+        dimnam(4:4)=chrid(k)
+        dimnam(5:5)='_'
+        dimnam(6:5+len_trim(varnam))=trim(varnam)
+        dimnam=dimnam(1:5+len_trim(varnam))
+        did(k)=-1
+        if (numdims.gt.0) then
+c         check if an existing dimension-declaration can be used
+c         instead of defining a new dimension
+          do i=1,numdims
+            dimchk=dimnams(i)
+            if ((vardim(k).eq.dimvals(i)).and.
+     &        (dimnam(1:4).eq.dimchk(1:4))) then
+              did(k)=i
+              goto 100
+            endif
+          enddo
+continue
+        endif
+        if (did(k).lt.0) then
+c         define the dimension
+          did(k)=ncddef(cdfid,dimnam,vardim(k),error)
+          if (error.ne.0) goto 920
+        endif
+      enddo
+c     define the times-array
+      if (ndim.eq.4) then
+c       define dimension and variable 'time'
+        if (numdims.ge.4) then
+          did(4)=ncdid(cdfid,'time',ierr)
+          idtime=ncvid(cdfid,'time',ierr)
+        else
+c         this dimension must first be defined
+          did(4) = ncddef (cdfid,'time',NCUNLIM,ierr)
+          idtime = ncvdef (cdfid,'time',NCFLOAT,1,did(4),ierr)
+        endif
+      endif
+c     define variable
+      id=ncvdef(cdfid,varnam,NCFLOAT,ndim,did,error)
+      if (error.ne.0) goto 920
+c     define attributes
+      do k=1,min0(ndim,3)
+c       min postion
+        attnam(1:1)=chrid(k)
+        attnam(2:4)='min'
+        attnam=attnam(1:4)
+        call ncapt(cdfid,id,attnam,NCFLOAT,1,varmin(k),error)
+        if (error.gt.0) goto 920
+c       max position
+        attnam(1:1)=chrid(k)
+        attnam(2:4)='max'
+        attnam=attnam(1:4)
+        call ncapt(cdfid,id,attnam,NCFLOAT,1,varmax(k),error)
+        if (error.gt.0) goto 920
+c       staggering
+        attnam(1:1)=chrid(k)
+        attnam(2:5)='stag'
+        attnam=attnam(1:5)
+        call ncapt(cdfid,id,attnam,NCFLOAT,1,stag(k),error)
+        if (error.gt.0) goto 920
+      enddo
+c     define missing data value
+      call ncapt(cdfid,id,'missing_data',NCFLOAT,1,misdat,error)
+      if (error.gt.0) goto 920
+c     leave define mode
+      call ncendf(cdfid,error)
+      if (error.gt.0) goto 920
+c     synchronyse output to disk and exit
+      call ncsnc (cdfid,error)
+      call ncpopt (ncopts)
+      return
+c     Error exits.
+write (6, *) '*ERROR*: When calling putcdf, the number of ',
+     &             'variable dimensions must be less or equal 4.'
+      call ncpopt (ncopts)
+      call ncclos (cdfid, error)
+      return
+write (6, *) '*ERROR*: An error occurred while attempting to ',
+     &             'write the data file in subroutine putcdf.'
+      call ncpopt (ncopts)
+      call ncclos (cdfid, error)
+      return
+      end
+      subroutine gettimes(cdfid,times,ntimes,ierr)
+C------------------------------------------------------------------------
+C     Purpose:
+C        Get all times on the specified NetCDF file
+C     Arguments:
+C        cdfid  int  input   identifier for NetCDF file
+C        times  real output  array contains all time values on the file,
+C                            dimensioned at least times(ntimes)
+C        ntimes int  output  number of times on the file
+C        error  int  output  errorflag
+C     History:
+C        Heini Wernli, ETHZ
+C------------------------------------------------------------------------
+      include "netcdf.inc"
+      integer   ierr,i
+      real times(*)
+      integer didtim,ntimes
+      integer   cdfid,idtime
+      integer   ncopts
+      character*(20) dimnam
+c     Get current value of error options, and make sure netCDF-errors do
+c     not abort execution
+      call ncgopt (ncopts)
+      call ncpopt(NCVERBOS)
+      didtim=ncdid(cdfid,'time',ierr)   ! inquire id for time dimension
+      if (ierr.ne.0) goto 900
+      idtime=ncvid(cdfid,'time',ierr)   ! inquire id for time array
+      if (ierr.ne.0) goto 900
+      call ncdinq(cdfid,didtim,dimnam,ntimes,ierr)      ! inquire # of times
+      if (ierr.ne.0) goto 900
+      do 10 i=1,ntimes
+        call ncvgt1(cdfid,idtime,i,times(i),ierr) ! get times
+        if (ierr.ne.0) goto 900
+continue
+c     normal exit
+      call ncpopt (ncopts)
+      return
+c     error exit
+ntimes=1
+      times(1)=0.
+      call ncpopt (ncopts)
+      end
+      subroutine puttimes(cdfid,times,ntimes,ierr)
+C------------------------------------------------------------------------
+C     Purpose:
+C        Get all times on the specified NetCDF file
+C     Arguments:
+C        cdfid  int  input   identifier for NetCDF file
+C        times  real input   array contains all time values on the file,
+C                            dimensioned at least times(ntimes)
+C        ntimes int  input   number of times on the file
+C        error  int  output  errorflag
+C     History:
+C        Heini Wernli, ETHZ
+C        Christoph Schaer, ETHZ
+C     Note:
+C        This preliminary version does not define the times-array, but only
+C        overwrites or extends an existing times-array.
+C------------------------------------------------------------------------
+      integer   ierr,i
+      real times(*)
+      integer didtim,ntimes
+      integer   cdfid,idtime,nfiltim
+      integer   ncdid,ncvid
+      idtime=ncvid(cdfid,'time',ierr)   ! inquire id for time array
+      if (ierr.ne.0) return
+      didtim=ncdid(cdfid,'time',ierr)   ! inquire id for time dimension
+      if (ierr.ne.0) return
+      call ncdinq(cdfid,didtim,'time',nfiltim,ierr)   ! inquire # of times
+      if (ierr.ne.0) return
+      if (nfiltim.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
+      do 10 i=1,ntimes
+        call ncvpt1(cdfid,idtime,i,times(i),ierr)
+        if (ierr.ne.0) return
+continue
+      end
+      subroutine cpp_crecdf(filnam,filnam_len,cdfid,phymin,phymax,ndim,
+     &     cfn,cfn_len,error)
+C------------------------------------------------------------------------
+C     Purpose:
+C        allows to call crecdf from c++
+C     Arguments:
+C        see crecdf
+C        additionally: fname_len and cfn_len, the length of the
+C           strings
+C
+C
+C     History:
+C        981221  Mark A. Liniger ETHZ
+C
+C     Note:
+C
+C
+C------------------------------------------------------------------------
+      integer        filnam_len,ndim,cfn_len,error,cdfid
+      character *(*) filnam,cfn
+      real           phymin(*),phymax(*)
+      call crecdf (filnam(1:filnam_len),cdfid,phymin,phymax,ndim,
+     &     cfn(1:cfn_len),error)
+      end
+      subroutine cpp_putdef(cdfid,varnam,varnam_len,ndim,misdat,
+     &     vardim,varmin,varmax,stag,error)
+C------------------------------------------------------------------------
+C     Purpose:
+C        allows to call putdef from c++
+C     Arguments:
+C        see crecdf
+C        additionally: varnam_len, the length of the
+C           strings
+C
+C
+C     History:
+C        981221  Mark A. Liniger ETHZ
+C
+C     Note:
+C
+C
+C------------------------------------------------------------------------
+      integer        varnam_len,ndim,error,vardim(*),cdfid
+      character *(*) varnam
+      real           misdat,stag(*),varmin(*),varmax(*)
+      call putdef (cdfid, varnam(1:varnam_len), ndim, misdat,
+     &     vardim, varmin, varmax, stag, error)
+      end
+      subroutine cpp_putdat(cdfid, varnam,varnam_len,
+     &     time, level, dat, error)
+C------------------------------------------------------------------------
+C     Purpose:
+C        allows to call putdef from c++
+C     Arguments:
+C        see crecdf
+C        additionally: varnam_len, the length of the
+C           strings
+C
+C
+C     History:
+C        981221  Mark A. Liniger ETHZ
+C
+C     Note:
+C
+C
+C------------------------------------------------------------------------
+      integer        varnam_len,cdfid,error,level
+      character *(*) varnam
+      real           dat(*)
+      real           time
+      call putdat(cdfid, varnam(1:varnam_len), time, level, dat, error)
+      end
+      FUNCTION strbeg (string)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This function returns the position of the first nonblank
+c           character in the input string.  Returns 0 if the entire
+c           string is blank.
+c     Arguments:
+c        string  char  input  string to be examined.
+c     History:
+c-----------------------------------------------------------------------
+      IMPLICIT NONE
+c     Function declaration
+      INTEGER strbeg
+c     Argument declarations
+      CHARACTER*(*) string
+c     Local variable declarations.
+      INTEGER i
+c     External function declarations.
+      INTEGER len
+      DO i = 1, len(string)
+        strbeg = i
+        IF (string(i:i) .NE. ' ' .AND. string(i:i) .NE. char(0)) THEN
+          RETURN
+        ENDIF
+      ENDDO
+      strbeg = 0
+      END
+      FUNCTION strend (string)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This function returns the position of the last nonblank
+c           character in the input string.  Returns 0 if the entire
+c           string is blank.
+c     Arguments:
+c        string  char  input  string to be examined.
+c     History:
+c-----------------------------------------------------------------------
+      IMPLICIT NONE
+c     Function declaration
+      INTEGER strend
+c     Argument declarations
+      CHARACTER*(*) string
+c     Local variable declarations.
+      INTEGER i
+c     External function declarations.
+      INTEGER len
+      DO i = len(string), 1, -1
+        strend = i
+        IF (string(i:i) .NE. ' ' .AND. string(i:i) .NE. char(0)) THEN
+          RETURN
+        ENDIF
+      ENDDO
+      strend = 0
+      END

 /trunk/lib/libcdfplus.f
 ,0 → 1,1128
+      subroutine wricst(cstnam,datar,aklev,bklev,aklay,bklay,stdate)
+C------------------------------------------------------------------------
+C     Creates the constants file for NetCDF files containing ECMWF
+C     data. The constants file is compatible with the one created
+C     for EM data (with subroutine writecst).
+C
+C     Input parameters:
+C
+C     cstnam    name of constants file
+C     datar     array contains all required parameters to write file
+C               datar(1):       number of points along x
+C               datar(2):       number of points along y
+C               datar(3):       maximum latitude of data region (ymax)
+C               datar(4):       minimum longitude of data region (xmin)
+C               datar(5):       minimum latitude of data region (ymin)
+C               datar(6):       maximum longitude of data region (xmax)
+C               datar(7):       grid increment along x
+C               datar(8):       grid increment along y
+C               datar(9):       number of levels
+C               datar(10):      data type (forecast or analysis)
+C               datar(11):      data version
+C               datar(12):      constants file version
+C               datar(13):      longitude of pole of coordinate system
+C               datar(14):      latitude of pole of coordinate system
+C     aklev     array contains the aklev values
+C     bklev     array contains the bklev values
+C     aklay     array contains the aklay values
+C     bklay     array contains the bklay values
+C     stdate    array contains date (year,month,day,time,step) of first
+C               field on file (start-date), dimensionised as stdate(5)
+C------------------------------------------------------------------------
+      include "netcdf.inc"
+      integer   nchar,maxlev
+      parameter (nchar=20,maxlev=32)
+      real      aklev(maxlev),bklev(maxlev)
+      real      aklay(maxlev),bklay(maxlev)
+      real      pollat,latmin,latmax
+      integer   datar(14)
+      integer   stdate(5)
+      character*80 cstnam
+C     declarations for constants-variables
+      integer   nz
+      integer   dattyp, datver, cstver
+C     further declarations
+      integer   ierr                    ! error flag
+      integer   cdfid                   ! NetCDF id
+      integer   xid,yid,zid             ! dimension ids
+      integer   pollonid, pollatid,     ! variable ids
+     >          aklevid, bklevid, aklayid, bklayid,
+     >          lonminid, lonmaxid, latminid, latmaxid,
+     >          dellonid, dellatid,
+     >          startyid, startmid, startdid, starthid, startsid,
+     >          dattypid, datverid, cstverid
+      nz=datar(9)                       ! number of levels
+C     Set data-type and -version, version of cst-file-format
+      dattyp=datar(10)
+      datver=datar(11)
+      cstver=datar(12)
+C     Initially set error to false
+      ierr=0
+C     Create constants file
+      cdfid=nccre(trim(cstnam),NCCLOB,ierr)
+C     Define the dimensions
+      xid = ncddef (cdfid,'nx',datar(1),ierr)
+      yid = ncddef (cdfid,'ny',datar(2),ierr)
+      zid = ncddef (cdfid,'nz',datar(9),ierr)
+C     Define integer constants
+      pollonid = ncvdef(cdfid,'pollon', NCFLOAT,0,0,ierr)
+      pollatid = ncvdef(cdfid,'pollat', NCFLOAT,0,0,ierr)
+      aklevid = ncvdef (cdfid, 'aklev', NCFLOAT, 1, zid, ierr)
+      bklevid = ncvdef (cdfid, 'bklev', NCFLOAT, 1, zid, ierr)
+      aklayid = ncvdef (cdfid, 'aklay', NCFLOAT, 1, zid, ierr)
+      bklayid = ncvdef (cdfid, 'bklay', NCFLOAT, 1, zid, ierr)
+      lonminid = ncvdef (cdfid, 'lonmin', NCFLOAT, 0, 0, ierr)
+      lonmaxid = ncvdef (cdfid, 'lonmax', NCFLOAT, 0, 0, ierr)
+      latminid = ncvdef (cdfid, 'latmin', NCFLOAT, 0, 0, ierr)
+      latmaxid = ncvdef (cdfid, 'latmax', NCFLOAT, 0, 0, ierr)
+      dellonid = ncvdef (cdfid, 'dellon', NCFLOAT, 0, 0, ierr)
+      dellatid = ncvdef (cdfid, 'dellat', NCFLOAT, 0, 0, ierr)
+      startyid = ncvdef (cdfid, 'starty', NCLONG, 0, 0, ierr)
+      startmid = ncvdef (cdfid, 'startm', NCLONG, 0, 0, ierr)
+      startdid = ncvdef (cdfid, 'startd', NCLONG, 0, 0, ierr)
+      starthid = ncvdef (cdfid, 'starth', NCLONG, 0, 0, ierr)
+      startsid = ncvdef (cdfid, 'starts', NCLONG, 0, 0, ierr)
+      dattypid = ncvdef (cdfid, 'dattyp', NCLONG, 0, 0, ierr)
+      datverid = ncvdef (cdfid, 'datver', NCLONG, 0, 0, ierr)
+      cstverid = ncvdef (cdfid, 'cstver', NCLONG, 0, 0, ierr)
+C     Leave define mode
+      call ncendf(cdfid,ierr)
+C     Store levels
+      call ncvpt(cdfid, aklevid, 1, nz, aklev, ierr)
+      call ncvpt(cdfid, bklevid, 1, nz, bklev, ierr)
+      call ncvpt(cdfid, aklayid, 1, nz, aklay, ierr)
+      call ncvpt(cdfid, bklayid, 1, nz, bklay, ierr)
+C     Store position of pole (trivial for ECMWF data)
+      call ncvpt1(cdfid, pollonid, 1, real(datar(13))/1000., ierr)
+      if (datar(14).gt.0) then
+        pollat=min(real(datar(14))/1000.,90.)
+      else
+        pollat=max(real(datar(14))/1000.,-90.)
+      endif
+      call ncvpt1(cdfid, pollatid, 1, pollat, ierr)
+C     Store horizontal data borders and grid increments
+      call ncvpt1(cdfid, lonminid, 1, real(datar(4))/1000., ierr)
+      call ncvpt1(cdfid, lonmaxid, 1, real(datar(6))/1000., ierr)
+      latmin=max(real(datar(5))/1000.,-90.)
+      latmax=min(real(datar(3))/1000.,90.)
+      call ncvpt1(cdfid, latminid, 1, latmin, ierr)
+      call ncvpt1(cdfid, latmaxid, 1, latmax, ierr)
+      call ncvpt1(cdfid, dellonid, 1, real(datar(7))/1000., ierr)
+      call ncvpt1(cdfid, dellatid, 1, real(datar(8))/1000., ierr)
+C     Store date of first field on file (start-date)
+      call ncvpt1(cdfid, startyid, 1, stdate(1), ierr)
+      call ncvpt1(cdfid, startmid, 1, stdate(2), ierr)
+      call ncvpt1(cdfid, startdid, 1, stdate(3), ierr)
+      call ncvpt1(cdfid, starthid, 1, stdate(4), ierr)
+      call ncvpt1(cdfid, startsid, 1, stdate(5), ierr)
+C     Store datatype and version
+      call ncvpt1(cdfid, dattypid, 1, dattyp, ierr)
+      call ncvpt1(cdfid, datverid, 1, datver, ierr)
+C     Store version of the constants file format
+      call ncvpt1(cdfid, cstverid, 1, cstver, ierr)
+C     Store strings
+      call ncclos(cdfid,ierr)
+      return
+      end
+      subroutine writelmcst(cdfid,nx,ny,nz,pollon,pollat,lonmin,
+     &lonmax,latmin,latmax,dellon,dellat,dattyp,datver,cstver,
+     &psref,tstar,tbeta,pintf,p0top,idate)
+c     ------------------------------------------------------------------
+      implicit none
+      integer   cdfid
+c     deklarationen der constants-variablen
+      real       pollon,pollat
+      real       lonmin,lonmax,latmin,latmax,dellon,dellat
+      integer    idate(5)
+      integer    nx,ny,nz
+      integer    dattyp, datver, cstver
+      real       psref, tstar, tbeta, pintf, p0top
+      include 'netcdf.inc'
+* netcdf declaration
+      integer   iret, k
+* dimension ids
+      integer  nxdim, nydim, nzdim
+* variable ids
+      integer  startyid, startmid, startdid, starthid
+* variable shapes, corners and edge lengths
+      integer dims(1), corner(1), edges(1)
+* enter define mode
+      call ncredf(cdfid, iret)
+      startyid = ncvdef (cdfid, 'starty', NCLONG, 0, 0, iret)
+      startmid = ncvdef (cdfid, 'startm', NCLONG, 0, 0, iret)
+      startdid = ncvdef (cdfid, 'startd', NCLONG, 0, 0, iret)
+      starthid = ncvdef (cdfid, 'starth', NCLONG, 0, 0, iret)
+* store the rest as global attributes
+* store nx,ny,nz
+      call ncapt(cdfid,NCGLOBAL,'nx',NCLONG,1,nx,iret)
+      call ncapt(cdfid,NCGLOBAL,'ny',NCLONG,1,ny,iret)
+      call ncapt(cdfid,NCGLOBAL,'nz',NCLONG,1,nz,iret)
+* store pollon, pollat
+      call ncapt(cdfid,NCGLOBAL,'pollon',NCFLOAT,1,pollon,iret)
+      call ncapt(cdfid,NCGLOBAL,'pollat',NCFLOAT,1,pollat,iret)
+* store lonmin, etc
+      call ncapt(cdfid,NCGLOBAL,'lonmin',NCFLOAT,1,lonmin,iret)
+      call ncapt(cdfid,NCGLOBAL,'lonmax',NCFLOAT,1,lonmax,iret)
+      call ncapt(cdfid,NCGLOBAL,'latmin',NCFLOAT,1,latmin,iret)
+      call ncapt(cdfid,NCGLOBAL,'latmax',NCFLOAT,1,latmax,iret)
+      call ncapt(cdfid,NCGLOBAL,'dellon',NCFLOAT,1,dellon,iret)
+      call ncapt(cdfid,NCGLOBAL,'dellat',NCFLOAT,1,dellat,iret)
+* store data type and version
+      call ncapt(cdfid,NCGLOBAL,'dattyp',NCLONG,1,dattyp,iret)
+      call ncapt(cdfid,NCGLOBAL,'datver',NCLONG,1,datver,iret)
+      call ncapt(cdfid,NCGLOBAL,'cstver',NCLONG,1,cstver,iret)
+* store information of lm model vertical grid
+      call ncapt(cdfid,NCGLOBAL,'psref',NCFLOAT,1,psref,iret)
+      call ncapt(cdfid,NCGLOBAL,'tstar',NCFLOAT,1,tstar,iret)
+      call ncapt(cdfid,NCGLOBAL,'tbeta',NCFLOAT,1,tbeta,iret)
+      call ncapt(cdfid,NCGLOBAL,'pintf',NCFLOAT,1,pintf,iret)
+      call ncapt(cdfid,NCGLOBAL,'p0top',NCFLOAT,1,p0top,iret)
+* leave define mode
+      call ncendf(cdfid, iret)
+* store starty, etc
+      corner(1) = 1
+      edges(1) = 1
+      call ncvpt(cdfid, startyid, corner, edges, idate(1), iret)
+      call ncvpt(cdfid, startmid, corner, edges, idate(2), iret)
+      call ncvpt(cdfid, startdid, corner, edges, idate(3), iret)
+      call ncvpt(cdfid, starthid, corner, edges, idate(4), iret)
+      end
+      subroutine globcst(cdfnam,datar,aklev,bklev,aklay,bklay,stdate)
+C------------------------------------------------------------------------
+C+
+C NAME:
+C     subroutine globcst
+C
+C PURPOSE:
+C     instead of writing a constants-file (*_cst), the information
+C     is added to the netCDF file as global variables
+C     the data format is compatible with the one requested by
+C     the IVE ETH/MIT version, contact author about details
+C
+C CATEGORY:
+C     model,netCDF
+C
+C CALLING SEQUENCE:
+C     subroutine globcst(cdfnam,datar,aklev,bklev,aklay,bklay,stdate)
+C
+C INPUTS:
+C     cdfnam    name of netCDF file
+C               The file needs to exist, otherwise an ERROR occurs,
+C               i.e. nothing is done
+C     datar     array contains all required parameters to write file
+C               datar(1):       number of points along x
+C               datar(2):       number of points along y
+C               datar(3):       maximum latitude of data region (ymax)
+C               datar(4):       minimum longitude of data region (xmin)
+C               datar(5):       minimum latitude of data region (ymin)
+C               datar(6):       maximum longitude of data region (xmax)
+C               datar(7):       grid increment along x
+C               datar(8):       grid increment along y
+C               datar(9):       number of levels
+C               datar(10):      data type (forecast or analysis)
+C               datar(11):      data version
+C               datar(12):      constants file version
+C               datar(13):      longitude of pole of coordinate system
+C               datar(14):      latitude of pole of coordinate system
+C     aklev     array contains the aklev values
+C     bklev     array contains the bklev values
+C     aklay     array contains the aklay values
+C     bklay     array contains the bklay values
+C     stdate    array contains date (year,month,day,time,step) of first
+C               field on file (start-date), dimensionised as stdate(5)
+C     list    the griblist-ASCII-file
+C     varno   the GRIB code number
+C
+C OUTPUTS:
+C     Adds cdf-information to EXISTING netCDF-file
+C
+C MODIFICATION HISTORY:
+C
+C     June  93    Christoph Schaer (ETHZ) created
+C     Nov   93    Heini Wernli (ETHZ) wricst
+C     Nov   98    David N. Bresch (MIT) wricst to globcst
+C-
+C     Sun include statement.
+      include "netcdf.inc"
+      integer   nchar,maxlev
+      parameter (nchar=20,maxlev=32)
+      real      aklev(maxlev),bklev(maxlev)
+      real      aklay(maxlev),bklay(maxlev)
+      integer   datar(14)
+      integer   stdate(5)
+      character*80 cdfnam
+C     declarations for constants-variables
+      integer   nz
+      integer   dattyp, datver, cstver
+C     further declarations
+      integer   ierr                    ! error flag
+      integer   cdfid                   ! NetCDF id
+      integer   xid,yid,zid             ! dimension ids
+      integer   pollonid, pollatid,     ! variable ids
+     >          aklevid, bklevid, aklayid, bklayid,
+     >          lonminid, lonmaxid, latminid, latmaxid,
+     >          dellonid, dellatid,
+     >          startyid, startmid, startdid, starthid, startsid,
+     >          dattypid, datverid, cstverid
+      nz=datar(9)                       ! number of levels
+C     Set data-type and -version, version of cst-file-format
+      dattyp=datar(10)
+      datver=datar(11)
+      cstver=datar(12)
+C     Initially set error to false
+      ierr=0
+C     open the netCDF-file:
+      call cdfwopn(cdfnam,cdfid,ierr)
+      if (ierr.ne.0) then
+         print*,'ERROR opening netCDF-file ',cdfnam
+         return
+      endif
+C     Put file into define mode
+      call ncredf(cdfid,ierr)
+      if (ierr.ne.0) then
+         print*,'ERROR switching to netCDF redefine mode'
+         return
+      endif
+C     Define the dimensions
+      xid = ncddef (cdfid,'nx',datar(1),ierr)
+      yid = ncddef (cdfid,'ny',datar(2),ierr)
+      zid = ncddef (cdfid,'nz',datar(9),ierr)
+C     Define integer constants
+      pollonid = ncvdef(cdfid,'pollon', NCFLOAT,0,0,ierr)
+      pollatid = ncvdef(cdfid,'pollat', NCFLOAT,0,0,ierr)
+      aklevid = ncvdef (cdfid, 'aklev', NCFLOAT, 1, zid, ierr)
+      bklevid = ncvdef (cdfid, 'bklev', NCFLOAT, 1, zid, ierr)
+      aklayid = ncvdef (cdfid, 'aklay', NCFLOAT, 1, zid, ierr)
+      bklayid = ncvdef (cdfid, 'bklay', NCFLOAT, 1, zid, ierr)
+      lonminid = ncvdef (cdfid, 'lonmin', NCFLOAT, 0, 0, ierr)
+      lonmaxid = ncvdef (cdfid, 'lonmax', NCFLOAT, 0, 0, ierr)
+      latminid = ncvdef (cdfid, 'latmin', NCFLOAT, 0, 0, ierr)
+      latmaxid = ncvdef (cdfid, 'latmax', NCFLOAT, 0, 0, ierr)
+      dellonid = ncvdef (cdfid, 'dellon', NCFLOAT, 0, 0, ierr)
+      dellatid = ncvdef (cdfid, 'dellat', NCFLOAT, 0, 0, ierr)
+      startyid = ncvdef (cdfid, 'starty', NCLONG, 0, 0, ierr)
+      startmid = ncvdef (cdfid, 'startm', NCLONG, 0, 0, ierr)
+      startdid = ncvdef (cdfid, 'startd', NCLONG, 0, 0, ierr)
+      starthid = ncvdef (cdfid, 'starth', NCLONG, 0, 0, ierr)
+      startsid = ncvdef (cdfid, 'starts', NCLONG, 0, 0, ierr)
+      dattypid = ncvdef (cdfid, 'dattyp', NCLONG, 0, 0, ierr)
+      datverid = ncvdef (cdfid, 'datver', NCLONG, 0, 0, ierr)
+      cstverid = ncvdef (cdfid, 'cstver', NCLONG, 0, 0, ierr)
+C     Leave define mode
+      call ncendf(cdfid,ierr)
+      if (ierr.ne.0) then
+         print*,'ERROR exiting define mode'
+         return
+      endif
+C     Store levels
+      call ncvpt(cdfid, aklevid, 1, nz, aklev, ierr)
+      call ncvpt(cdfid, bklevid, 1, nz, bklev, ierr)
+      call ncvpt(cdfid, aklayid, 1, nz, aklay, ierr)
+      call ncvpt(cdfid, bklayid, 1, nz, bklay, ierr)
+C     Store position of pole (trivial for ECMWF data)
+      call ncvpt1(cdfid, pollonid, 1, real(datar(13))/1000., ierr)
+      call ncvpt1(cdfid, pollatid, 1, real(datar(14))/1000., ierr)
+C     Store horizontal data borders and grid increments
+      call ncvpt1(cdfid, lonminid, 1, real(datar(4))/1000., ierr)
+      call ncvpt1(cdfid, lonmaxid, 1, real(datar(6))/1000., ierr)
+      call ncvpt1(cdfid, latminid, 1, real(datar(5))/1000., ierr)
+      call ncvpt1(cdfid, latmaxid, 1, real(datar(3))/1000., ierr)
+      call ncvpt1(cdfid, dellonid, 1, real(datar(7))/1000., ierr)
+      call ncvpt1(cdfid, dellatid, 1, real(datar(8))/1000., ierr)
+C     Store date of first field on file (start-date)
+      call ncvpt1(cdfid, startyid, 1, stdate(1), ierr)
+      call ncvpt1(cdfid, startmid, 1, stdate(2), ierr)
+      call ncvpt1(cdfid, startdid, 1, stdate(3), ierr)
+      call ncvpt1(cdfid, starthid, 1, stdate(4), ierr)
+      call ncvpt1(cdfid, startsid, 1, stdate(5), ierr)
+C     Store datatype and version
+      call ncvpt1(cdfid, dattypid, 1, dattyp, ierr)
+      call ncvpt1(cdfid, datverid, 1, datver, ierr)
+C     Store version of the constants file format
+      call ncvpt1(cdfid, cstverid, 1, cstver, ierr)
+      if (ierr.ne.0) then
+         print*,'ERROR adding cst-date as global variables'
+         return
+      endif
+C     Store strings
+      call ncclos(cdfid,ierr)
+      if (ierr.ne.0) then
+         print*,'ERROR closing netCDF file'
+      endif
+      return
+      end
+      subroutine getsdat(cdfid,varnam,time,ix,iy,iz,sx,sy,sz,dat,error)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This routine is called to read the data within a selected
+c        domain of a variable from an IVE-NetCDF file.
+c        Prior to calling this routine, the file must be opened with
+c        a call to opncdf (for extension) or crecdf (for creation) or
+c        readcdf (for readonly).
+c     Arguments:
+c        cdfid   int   input   file-identifier
+c                              (must be obtained by calling routine
+c                              opncdf,readcdf  or crecdf)
+c        varnam  char  input   the user-supplied variable name
+c        time    real  input   the user-supplied time-level of the
+c                              data to be read from the file (the time-
+c                              levels stored in the file can be obtained
+c                              with a call to gettimes).
+c        ix/y/z  int   input   indices of lower left corner of selected
+c                              data volume.
+c        sx/y/z  int   input   size of selected data volume
+c        dat     real  output  data-array with dimensions (sx,sy,sz).
+c        error   int   output  indicates possible errors found in this
+c                              routine.
+c                              error = 0   no errors detected.
+c                              error = 1   the variable is not present on
+c                                          the file.
+c                              error = 2   the value of 'time' is not
+c                                          known.to the file.
+c                              error = 6,7,8   data volume too large
+c                              error =10   another error.
+c     History:
+c       June  93    Christoph Schaer (ETHZ)  Created getdat
+c       Nov   93    Heini Wernli (ETHZ)      Created getsdat
+c-----------------------------------------------------------------------
+      include "netcdf.inc"
+C     Declaration of local variables
+      character*(*) varnam
+      character*(20) chars
+      integer cdfid
+      integer     ix,iy,iz,sx,sy,sz
+      real        dat(sx,sy,sz)
+      real        misdat,varmin(4),varmax(4),stag(4)
+      real        time, timeval
+      integer     corner(4),edgeln(4),didtim,vardim(4),ndims
+      integer     error, ierr
+      integer     ntime
+      integer     idtime,idvar,iflag
+      integer     i
+      call ncpopt(NCVERBOS)
+c     access the variable
+      call getdef (cdfid, trim(varnam), ndims, misdat,
+     &                           vardim, varmin, varmax, stag, ierr)
+      if (ierr.ne.0) then
+        print *,'*ERROR* in getdef in getdat'
+        error=1
+        return
+      endif
+      idvar=ncvid(cdfid,trim(varnam),ierr)
+      if (ierr.ne.0) then
+        print *,'*ERROR* in ncvid in getsdat'
+        error=1
+        return
+      endif
+C     Get times-array
+      didtim=ncdid(cdfid,'time',ierr)
+      if (ierr.ne.0) then
+        print *,'*ERROR* didtim in getsdat'
+        error=10
+        return
+      endif
+      call ncdinq(cdfid,didtim,chars,ntime,ierr)
+      if (ierr.ne.0) then
+        print *,'*ERROR* in ncdinq in getsdat'
+        error=10
+        return
+      endif
+      idtime=ncvid(cdfid,'time',ierr)
+      if (ierr.ne.0) then
+        print *,'*ERROR* in ncvid for time in getsdat'
+        error=10
+        return
+      endif
+c     find appropriate time-index
+      iflag=0
+      do i=1,ntime
+        call ncvgt1(cdfid,idtime,i,timeval,ierr)
+        if (ierr.ne.0) print *,'*ERROR* in ncvgt1 in getsdat'
+        if (time.eq.timeval) iflag=i
+      enddo
+      if (iflag.eq.0) then
+        error=2
+        print *,'Error: Unknown time in getsdat'
+        print *,time,timeval
+        return
+      endif
+C     Define data volume to be written (index space)
+      corner(1)=ix
+      corner(2)=iy
+      corner(3)=iz
+      corner(4)=iflag
+      edgeln(1)=sx
+      edgeln(2)=sy
+      edgeln(3)=sz
+      edgeln(4)=1
+C     Check if data volume is within data domain
+      if (ix+sx-1.gt.vardim(1)) then
+        error=7
+        print *,'Error: data volume too large in x-direction'
+        print *,ix,sx,vardim(1)
+        return
+      endif
+      if (iy+sy-1.gt.vardim(2)) then
+        error=8
+        print *,'Error: data volume too large in y-direction'
+        return
+      endif
+      if (iz+sz-1.gt.vardim(3)) then
+        error=9
+        print *,'Error: data volume too large in z-direction'
+        return
+      endif
+C     Read data from NetCDF file
+      call ncvgt(cdfid,idvar,corner,edgeln,dat,error)
+      if (error.ne.0) then
+        print *, 'corner ',corner(1),corner(2),corner(3)
+        print *, 'edgeln ',edgeln(1),edgeln(2),edgeln(3)
+        print *, '*ERROR* in ncvgt in getsdat'
+        error=10
+      endif
+      end
+      subroutine getlevs(cstid,nlev,aklev,bklev,aklay,bklay,error)
+c-----------------------------------------------------------------------
+c     Purpose:
+c       This routine is called to get the level arrays aklev and
+c       bklev from a NetCDF constants file.
+c     Arguments:
+c       cstid     int   input   identifier for NetCDF constants file
+c       nlev      int   input   number of levels
+c       aklev     real  output  array contains all aklev values
+c       bklev     real  output  array contains all bklev values
+c       aklay     real  output  array contains all aklay values
+c       bklay     real  output  array contains all bklay values
+c       error     int   output  error flag
+c                               error = 0   no errors detected
+c                               error = 1   error detected
+c     History:
+c       Aug. 93   Heini Wernli          Created.
+c-----------------------------------------------------------------------
+      integer   error
+      integer   cstid
+      integer   ncdid,ncvid             ! NetCDF functions
+      integer   didz,idak,idbk,idaky,idbky
+      integer   nlev
+      real      aklev(nlev),bklev(nlev),aklay(nlev),bklay(nlev)
+      character*(20) dimnam
+      integer   i
+      didz      =ncdid(cstid,'nz',error)
+      if (error.ne.0) goto 920
+      idak      =ncvid(cstid,'aklev',error)
+      if (error.ne.0) goto 920
+      idbk      =ncvid(cstid,'bklev',error)
+      if (error.ne.0) goto 920
+      idaky     =ncvid(cstid,'aklay',error)
+      if (error.ne.0) goto 920
+      idbky     =ncvid(cstid,'bklay',error)
+      if (error.ne.0) goto 920
+      call ncdinq(cstid,didz,dimnam,nlev,error) ! read number of levels
+      if (error.ne.0) goto 920
+      do 10 i=1,nlev
+        call ncvgt1(cstid,idak,i,aklev(i),error)      ! get aklev
+        call ncvgt1(cstid,idbk,i,bklev(i),error)      ! get bklev
+        call ncvgt1(cstid,idaky,i,aklay(i),error)      ! get aklay
+        call ncvgt1(cstid,idbky,i,bklay(i),error)      ! get bklay
+        if (error.ne.0) goto 920
+continue
+      return
+c     Error exits.
+write(*,*)'*ERROR*: An error occured in subroutine getlevs'
+      return
+      end
+      subroutine getntim(cdfid,ntimes,ierr)
+C------------------------------------------------------------------------
+C     Purpose:
+C        Get number of times on the specified NetCDF file
+C     Arguments:
+C        cdfid  int  input   identifier for NetCDF file
+C        ntimes int  output  number of times on the file
+C        error  int  output  errorflag
+C     History:
+C        Heini Wernli, ETHZ
+C------------------------------------------------------------------------
+      include "netcdf.inc"
+      integer   ierr
+      integer didtim,ntimes
+      integer   cdfid,idtime
+      integer   ncopts
+      character*(20) dimnam
+c     Get current value of error options, and make sure netCDF-errors do
+c     not abort execution
+      call ncgopt (ncopts)
+      call ncpopt(NCVERBOS)
+      didtim=ncdid(cdfid,'time',ierr)   ! inquire id for time dimension
+      if (ierr.ne.0) goto 900
+      idtime=ncvid(cdfid,'time',ierr)   ! inquire id for time array
+      if (ierr.ne.0) goto 900
+      call ncdinq(cdfid,didtim,dimnam,ntimes,ierr)      ! inquire # of times
+      if (ierr.ne.0) goto 900
+c     normal exit
+      call ncpopt (ncopts)
+      return
+c     error exit
+ntimes=1
+      call ncpopt (ncopts)
+      end
+      subroutine getstart(cdfid,idate,ierr)
+C------------------------------------------------------------------------
+C     Purpose:
+C       Get start date for fields on specified NetCDF file
+C     Arguments:
+C       cdfid   int     input   identifier for NetCDF file
+C       idate   int     output  array contains date (year,month,day,time,step)
+C                               dimensioned as idate(5)
+C       ierr    int     output  error flag
+C------------------------------------------------------------------------
+      include "netcdf.inc"
+c     variable declarations
+      integer   ierr
+      integer   idate(5)
+      integer   cdfid,ncopts,idvar,nvars
+      integer   ndims,ngatts,recdim,i,vartyp,nvatts,vardim(4)
+      character*20 vnam(100)
+c     Get current value of error options, and make sure NetCDF-errors do
+c     not abort execution
+      call ncgopt (ncopts)
+      call ncpopt (NCVERBOS)
+      idvar=ncvid(cdfid,'starty',ierr)
+      if (ierr.ne.0) goto 930
+      call ncvgt1(cdfid,idvar,1,idate(1),ierr)
+      if (ierr.ne.0) goto 920
+      idvar=ncvid(cdfid,'startm',ierr)
+      if (ierr.ne.0) goto 920
+      call ncvgt1(cdfid,idvar,1,idate(2),ierr)
+      if (ierr.ne.0) goto 920
+      idvar=ncvid(cdfid,'startd',ierr)
+      if (ierr.ne.0) goto920
+      call ncvgt1(cdfid,idvar,1,idate(3),ierr)
+      if (ierr.ne.0) goto 920
+      idvar=ncvid(cdfid,'starth',ierr)
+      if (ierr.ne.0) goto 920
+      call ncvgt1(cdfid,idvar,1,idate(4),ierr)
+      if (ierr.ne.0) goto 920
+C     Starts is not defined on all files
+C     Only ask for it if it exists
+C     Inquire number of dimensions, variables and attributes
+      idate(5)=0
+      call ncinq(cdfid,ndims,nvars,ngatts,recdim,ierr)
+      do i=1,nvars
+        call ncvinq(cdfid,i,vnam(i),vartyp,ndims,vardim,nvatts,ierr)
+        if (vnam(i).eq.'starts') then
+          idvar=ncvid(cdfid,'starts',ierr)
+          call ncvgt1(cdfid,idvar,1,idate(5),ierr)
+          if (ierr.ne.0) goto 920
+        endif
+      enddo
+c     normal exit
+      call ncpopt (ncopts)
+      return
+c     error exit
+continue
+      write (6, *) 'ERROR: An error occurred while attempting to ',
+     &             'read the starting-time in subroutine putstart.'
+continue
+      call ncpopt (ncopts)
+      end
+      subroutine putstart(cdfid,idate,ierr)
+C----------------------------------------------------------------------
+C     Purpose:
+C        Puts the 'starting-time' on the specified NetCDF file.
+C     Arguments:
+C        cdfid   int     input   identifier for NetCDF file
+C        idate   int     input   array contains date (year,month,day,time,step)
+C                                dimensioned as idate(5)
+C        ierr    int     output  error flag
+C------------------------------------------------------------------------
+      include "netcdf.inc"
+c     variable declarations
+      integer   ierr,idate(5),startid(5),cdfid,ncopts,i
+c     Get current value of error options, and make sure NetCDF-errors do
+c     not abort execution
+      call ncgopt (ncopts)
+      call ncpopt (NCVERBOS)
+c     define variables
+      call ncredf(cdfid,ierr)
+      if (ierr.ne.0) goto 920
+      startid(1) = ncvdef (cdfid, 'starty', NCLONG, 0, 0, ierr)
+      if (ierr.ne.0) goto 920
+      startid(2) = ncvdef (cdfid, 'startm', NCLONG, 0, 0, ierr)
+      if (ierr.ne.0) goto 920
+      startid(3) = ncvdef (cdfid, 'startd', NCLONG, 0, 0, ierr)
+      if (ierr.ne.0) goto 920
+      startid(4) = ncvdef (cdfid, 'starth', NCLONG, 0, 0, ierr)
+      if (ierr.ne.0) goto 920
+      startid(5) = ncvdef (cdfid, 'starts', NCLONG, 0, 0, ierr)
+      if (ierr.ne.0) goto 920
+      call ncendf(cdfid, ierr)
+      if (ierr.ne.0) goto 920
+c     store variables
+      do i=1,5
+        call ncvpt1(cdfid,startid(i),1,idate(i),ierr)
+        if (ierr.ne.0) goto 920
+      enddo
+c     synchronyse output to disk, revert to previous error-mode, and exit
+      call ncsnc (cdfid,ierr)
+      call ncpopt (ncopts)
+      return
+c     error exit
+write (6, *) 'ERROR: An error occurred while attempting to ',
+     &             'write the starting-time in subroutine putstart.'
+      call ncpopt (ncopts)
+      call ncclos (cdfid, ierr)
+      end
+      subroutine getgrid(cdfid,dx,dy,ierr)
+C------------------------------------------------------------------------
+C     Purpose:
+C       Get grid increments for fields on specified NetCDF file
+C     Arguments:
+C       cdfid   int     input   identifier for NetCDF file
+C       dx      real    output  grid increment along latitude
+C       dy      real    output  grid increment along longitude
+C       ierr    int     output  error flag
+C------------------------------------------------------------------------
+      integer   ierr
+      integer   cdfid
+      integer   ncvid
+      integer   idilon,idilat
+      real      dx,dy
+      idilon    =ncvid(cdfid,'dellon',ierr)
+      if (ierr.ne.0) return
+      idilat    =ncvid(cdfid,'dellat',ierr)
+      if (ierr.ne.0) return
+      call ncvgt1(cdfid,idilon,1,dx,ierr)
+      if (ierr.ne.0) return
+      call ncvgt1(cdfid,idilat,1,dy,ierr)
+      if (ierr.ne.0) return
+      end
+      subroutine getdattyp(cdfid,typ,ierr)
+C------------------------------------------------------------------------
+C     Purpose:
+C       Get data type for specified NetCDF file
+C     Arguments:
+C       cdfid   int     input   identifier for NetCDF file
+C       typ     int     output  data type: 1 (52) for pressure (theta) coord
+C       ierr    int     output  error flag
+C------------------------------------------------------------------------
+      integer   ierr
+      integer   cdfid
+      integer   ncvid
+      integer   idtyp,typ
+      idtyp    =ncvid(cdfid,'dattyp',ierr)
+      if (ierr.ne.0) return
+      call ncvgt1(cdfid,idtyp,1,typ,ierr)
+      if (ierr.ne.0) return
+      end
+      subroutine getpole(cdfid,pollon,pollat,ierr)
+C------------------------------------------------------------------------
+C     Purpose:
+C       Get physical coordinates of pole of coordinate system
+C     Arguments:
+C       cdfid   int     input   identifier for NetCDF file
+C       pollon  real    output  longitude of pole
+C       pollat  real    output  latitude of pole
+C       ierr    int     output  error flag
+C------------------------------------------------------------------------
+      integer   ierr
+      integer   cdfid
+      integer   ncvid
+      integer   idplon,idplat
+      real      pollon,pollat
+      idplon    =ncvid(cdfid,'pollon',ierr)
+      if (ierr.ne.0) return
+      idplat    =ncvid(cdfid,'pollat',ierr)
+      if (ierr.ne.0) return
+      call ncvgt1(cdfid,idplon,1,pollon,ierr)
+      if (ierr.ne.0) return
+      call ncvgt1(cdfid,idplat,1,pollat,ierr)
+      if (ierr.ne.0) return
+      end
+      subroutine getmc2grid(cdfid,polx,poly,delx,shem,phi0,lam0,ierr)
+C------------------------------------------------------------------------
+C     Purpose:
+C       Get physical coordinates of pole of coordinate system
+C     Arguments:
+C       cdfid   int     input   identifier for NetCDF file
+C       ierr    int     output  error flag
+C------------------------------------------------------------------------
+      integer   ierr
+      integer   cdfid
+      integer   ncvid
+      integer   idpolx,idpoly,iddelx,idshem,idphi0,idlam0
+      real      polx,poly,delx,shem,phi0,lam0
+      idpolx    =ncvid(cdfid,'polx',ierr)
+      if (ierr.ne.0) return
+      idpoly    =ncvid(cdfid,'poly',ierr)
+      if (ierr.ne.0) return
+      iddelx    =ncvid(cdfid,'delx',ierr)
+      if (ierr.ne.0) return
+      idshem    =ncvid(cdfid,'shem',ierr)
+      if (ierr.ne.0) return
+      idphi0    =ncvid(cdfid,'phi0',ierr)
+      if (ierr.ne.0) return
+      idlam0    =ncvid(cdfid,'lam0',ierr)
+      if (ierr.ne.0) return
+      call ncvgt1(cdfid,idpolx,1,polx,ierr)
+      if (ierr.ne.0) return
+      call ncvgt1(cdfid,idpoly,1,poly,ierr)
+      if (ierr.ne.0) return
+      call ncvgt1(cdfid,iddelx,1,delx,ierr)
+      if (ierr.ne.0) return
+      call ncvgt1(cdfid,idshem,1,shem,ierr)
+      if (ierr.ne.0) return
+      call ncvgt1(cdfid,idphi0,1,phi0,ierr)
+      if (ierr.ne.0) return
+      call ncvgt1(cdfid,idlam0,1,lam0,ierr)
+      if (ierr.ne.0) return
+      end
+      subroutine getcfn(cdfid,cfn,ierr)
+C------------------------------------------------------------------------
+C     Purpose:
+C       Get name of constants file
+C     Arguments:
+C       cdfid   int     input   identifier for NetCDF file
+C       cfn     char    output  name of constants file
+C       ierr    int     output  error flag
+C------------------------------------------------------------------------
+      include 'netcdf.inc'
+      integer   ierr
+      integer   cdfid,lenstr
+      character*80 cfn
+      lenstr=80
+      call ncagtc(cdfid,NCGLOBAL,"constants_file_name",cfn,lenstr,ierr)
+      if (ierr.ne.0) write(*,*)'error in SR getcfn'
+      end
+      subroutine getdim (cdfid, varnam, nx, ny, nz, error)
+c-------------------------------------------------------------------------
+c     Purpose:
+c        This routine is called to get the dimensions of
+c        a variable from an IVE-NetCDF file for use with the IVE plotting
+c        package. Prior to calling this routine, the file must be opened
+c        with a call to opncdf.
+c     Arguments:
+c        cdfid   int   input   file-identifier
+c                              (can be obtained by calling routine
+c                              opncdf)
+c        varnam  char  input   the user-supplied variable name.
+c                              (can be obtained by calling routine
+c                              opncdf)
+c        nx      int   output  the zonal dimension of the variable.
+c        ny      int   output  the meridional dimension of the variable.
+c        nz      int   output  the vertical dimension of the variable.
+c
+c        error   int   output  indicates possible errors found in this
+c                              routine.
+c                              error = 0   no errors detected.
+c                              error = 1   the variable is not on the file.
+c                              error =10   other errors.
+c     History:
+c        March 2000    Heini Wernli (ETHZ)     Created.
+c-----------------------------------------------------------------------
+      include  "netcdf.inc"
+c     Argument declarations.
+      character *(*) varnam
+      integer        vardim(4), ndim, error, cdfid
+      integer        nx,ny,nz
+c     Local variable declarations.
+      character *(20) dimnam(MAXNCDIM),vnam
+      integer         id,i,k
+      integer         ndims,nvars,ngatts,recdim,dimsiz(MAXNCDIM)
+      integer         vartyp,nvatts, ncopts
+c     Get current value of error options.
+      call ncgopt (ncopts)
+c     make sure NetCDF-errors do not abort execution
+      call ncpopt(NCVERBOS)
+c     Initially set error to indicate no errors.
+      error = 0
+c     inquire for number of dimensions
+      call ncinq(cdfid,ndims,nvars,ngatts,recdim,error)
+      if (error.eq.1) goto 920
+c     read dimension-table
+      do i=1,ndims
+        call ncdinq(cdfid,i,dimnam(i),dimsiz(i),error)
+        if (error.gt.0) goto 920
+      enddo
+c     get id of the variable
+      id=ncvid(cdfid,varnam,error)
+      if (error.eq.1) goto 910
+c     inquire about variable
+      call ncvinq(cdfid,id,vnam,vartyp,ndim,vardim,nvatts,error)
+      if (vartyp.ne.NCFLOAT) error=1
+      if (error.gt.0) goto 920
+c     get dimensions from dimension-table
+      do k=1,ndim
+        vardim(k)=dimsiz(vardim(k))
+      enddo
+      nx=vardim(1)
+      ny=vardim(2)
+      nz=vardim(3)
+c     normal exit
+      call ncpopt (ncopts)
+      return
+c     Error exits.
+write (6, *) '*ERROR*: The selected variable could not be found ',
+     &             'in the file by getdim.'
+      call ncpopt (ncopts)
+      call ncclos (cdfid, error)
+      return
+write (6, *) '*ERROR*: An error occurred while attempting to ',
+     &             'read the data file in subroutine getcdf.'
+      call ncpopt (ncopts)
+      call ncclos (cdfid, error)
+      end
+      subroutine new_gettra(cdfid,varnam,ix,ntimes,array,ierr)
+C------------------------------------------------------------------------
+C
+C     Reads the time-evolution for one grid-point of the variable
+C     indicated by varnam.
+C
+C     cdfid     int     input   identifier for NetCDF file
+C     varnam    char    input   name of variable
+C     ix        int     input   index for trajectory to read
+C     ntimes    int     input   number of time-indices to read
+C     array     real    output  array contains the readed values
+C     ierr      int     output  error flag
+C------------------------------------------------------------------------
+C     Declaration of attributes
+      integer   cdfid
+      character*(*) varnam
+      integer   ix
+      integer   ntimes
+      real      array(ntimes)
+C     Declaration of local variables
+      integer   corner(4),edgeln(4)
+      integer   idvar,ierr
+      integer   ncvid
+      integer   strend
+      corner(1)=ix
+      corner(2)=1
+      corner(3)=1
+      corner(4)=1
+      edgeln(1)=1
+      edgeln(2)=1
+      edgeln(3)=1
+      edgeln(4)=ntimes
+      idvar =ncvid(cdfid,varnam(1:strend(varnam)),ierr)
+      call ncvgt(cdfid,idvar,corner,edgeln,array,ierr)
+      if (ierr.ne.0) goto 991
+      return
+stop 'Variable not found on NetCDF file in SR new_gettra'
+      end
+      subroutine getvars(cdfid,nvars,vnam,ierr)
+C------------------------------------------------------------------------
+C     Opens the NetCDF file 'filnam' and returns its identifier cdfid.
+C     filnam    char    input   name of NetCDF file to open
+C     nvars     int     output  number of variables on file
+C     vnam      char    output  array with variable names
+C     ierr      int     output  error flag
+C------------------------------------------------------------------------
+      include 'netcdf.inc'
+      integer   cdfid,ierr,nvars
+      character*(*) vnam(*)
+      integer   ndims,ngatts,recdim,i,vartyp,nvatts,vardim(4)
+      call ncpopt(NCVERBOS)
+C     Inquire number of dimensions, variables and attributes
+      call ncinq(cdfid,ndims,nvars,ngatts,recdim,ierr)
+C     Inquire variable names from NetCDF file
+      do i=1,nvars
+        call ncvinq(cdfid,i,vnam(i),vartyp,ndims,vardim,nvatts,ierr)
+      enddo
+      return
+      end

 /trunk/lib/libgm2em.f
 ,0 → 1,229
+      REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
+C
+C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
+C
+C**** PHTOPHS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
+C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
+C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
+C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
+C**   AUFRUF   :   PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM)
+C**   ENTRIES  :   KEINE
+C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF
+C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
+C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
+C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
+C**   VERSIONS-
+C**   DATUM    :   03.05.90
+C**
+C**   EXTERNALS:   KEINE
+C**   EINGABE-
+C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
+C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
+C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
+C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
+C**   AUSGABE-
+C**   PARAMETER:   ROTIERTE BREITE PHIS ALS WERT DER FUNKTION
+C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
+C**
+C**   COMMON-
+C**   BLOECKE  :   KEINE
+C**
+C**   FEHLERBE-
+C**   HANDLUNG :   KEINE
+C**   VERFASSER:   G. DE MORSIER
+      REAL        LAM,PHI,POLPHI,POLLAM
+      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
+      ZSINPOL = SIN(ZPIR18*POLPHI)
+      ZCOSPOL = COS(ZPIR18*POLPHI)
+      ZLAMPOL = ZPIR18*POLLAM
+      ZPHI    = ZPIR18*PHI
+      ZLAM    = LAM
+      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
+      ZLAM    = ZPIR18*ZLAM
+      ZARG    = ZCOSPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL) + ZSINPOL*SIN(ZPHI)
+      PHTOPHS = ZRPI18*ASIN(ZARG)
+      RETURN
+      END
+      REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
+C
+C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
+C
+C**** PHSTOPH  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
+C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
+C****                 ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
+C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
+C**   AUFRUF   :   PHI = PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
+C**   ENTRIES  :   KEINE
+C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
+C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
+C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
+C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
+C**   VERSIONS-
+C**   DATUM    :   03.05.90
+C**
+C**   EXTERNALS:   KEINE
+C**   EINGABE-
+C**   PARAMETER:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
+C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
+C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
+C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
+C**   AUSGABE-
+C**   PARAMETER:   WAHRE GEOGRAPHISCHE BREITE ALS WERT DER FUNKTION
+C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
+C**
+C**   COMMON-
+C**   BLOECKE  :   KEINE
+C**
+C**   FEHLERBE-
+C**   HANDLUNG :   KEINE
+C**   VERFASSER:   D.MAJEWSKI
+      REAL        LAMS,PHIS,POLPHI,POLLAM
+      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
+      SINPOL = SIN(ZPIR18*POLPHI)
+      COSPOL = COS(ZPIR18*POLPHI)
+      ZPHIS  = ZPIR18*PHIS
+      ZLAMS  = LAMS
+      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
+      ZLAMS  = ZPIR18*ZLAMS
+      ARG     = COSPOL*COS(ZPHIS)*COS(ZLAMS) + SINPOL*SIN(ZPHIS)
+      PHSTOPH = ZRPI18*ASIN(ARG)
+      RETURN
+      END
+      REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)
+C
+C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
+C
+C**** LMTOLMS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM
+C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
+C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
+C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
+C**   AUFRUF   :   LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM)
+C**   ENTRIES  :   KEINE
+C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF
+C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
+C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
+C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
+C**   VERSIONS-
+C**   DATUM    :   03.05.90
+C**
+C**   EXTERNALS:   KEINE
+C**   EINGABE-
+C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
+C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
+C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
+C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
+C**   AUSGABE-
+C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
+C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
+C**
+C**   COMMON-
+C**   BLOECKE  :   KEINE
+C**
+C**   FEHLERBE-
+C**   HANDLUNG :   KEINE
+C**   VERFASSER:   G. DE MORSIER
+      REAL        LAM,PHI,POLPHI,POLLAM
+      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
+      ZSINPOL = SIN(ZPIR18*POLPHI)
+      ZCOSPOL = COS(ZPIR18*POLPHI)
+      ZLAMPOL =     ZPIR18*POLLAM
+      ZPHI    =     ZPIR18*PHI
+      ZLAM    = LAM
+      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
+      ZLAM    = ZPIR18*ZLAM
+      ZARG1   = - SIN(ZLAM-ZLAMPOL)*COS(ZPHI)
+      ZARG2   = - ZSINPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL)+ZCOSPOL*SIN(ZPHI)
+      IF (ABS(ZARG2).LT.1.E-30) THEN
+        IF (ABS(ZARG1).LT.1.E-30) THEN
+          LMTOLMS =   0.0
+        ELSEIF (ZARG1.GT.0.) THEN
+              LMTOLMS =  90.0
+            ELSE
+              LMTOLMS = -90.0
+            ENDIF
+      ELSE
+        LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
+      ENDIF
+      RETURN
+      END
+      REAL FUNCTION LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
+C
+C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
+C
+C**** LMSTOLM  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
+C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
+C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
+C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
+C**   AUFRUF   :   LAM = LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
+C**   ENTRIES  :   KEINE
+C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
+C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
+C**                IM ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
+C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
+C**   VERSIONS-
+C**   DATUM    :   03.05.90
+C**
+C**   EXTERNALS:   KEINE
+C**   EINGABE-
+C**   PARAMETER:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
+C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
+C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
+C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
+C**   AUSGABE-
+C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
+C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
+C**
+C**   COMMON-
+C**   BLOECKE  :   KEINE
+C**
+C**   FEHLERBE-
+C**   HANDLUNG :   KEINE
+C**   VERFASSER:   D.MAJEWSKI
+      REAL        LAMS,PHIS,POLPHI,POLLAM
+      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
+      ZSINPOL = SIN(ZPIR18*POLPHI)
+      ZCOSPOL = COS(ZPIR18*POLPHI)
+      ZLAMPOL = ZPIR18*POLLAM
+      ZPHIS   = ZPIR18*PHIS
+      ZLAMS   = LAMS
+      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
+      ZLAMS   = ZPIR18*ZLAMS
+      ZARG1   = SIN(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
+ZCOSPOL*           SIN(ZPHIS)) -
+COS(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
+      ZARG2   = COS(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
+ZCOSPOL*           SIN(ZPHIS)) +
+SIN(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
+      IF (ABS(ZARG2).LT.1.E-30) THEN
+        IF (ABS(ZARG1).LT.1.E-30) THEN
+          LMSTOLM =   0.0
+        ELSEIF (ZARG1.GT.0.) THEN
+              LMSTOLAM =  90.0
+            ELSE
+              LMSTOLAM = -90.0
+            ENDIF
+      ELSE
+        LMSTOLM = ZRPI18*ATAN2(ZARG1,ZARG2)
+      ENDIF
+      RETURN
+      END

 /trunk/lib/libipo.f
 ,0 → 1,1368
+      real function int2d(ar,n1,n2,rid,rjd)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This subroutine interpolates a 2d-array to an arbitrary
+c        location within the grid.
+c     Arguments:
+c        ar        real  input   surface pressure, define as ar(n1,n2)
+c        n1,n2     int   input   dimensions of ar
+c        ri,rj     real  input   grid location to be interpolated to
+c     History:
+c-----------------------------------------------------------------------
+c     argument declarations
+      integer   n1,n2
+      real      ar(n1,n2), rid,rjd
+c     local declarations
+      integer   i,j,ip1,jp1,ih,jh
+      real      frac0i,frac0j,frac1i,frac1j,ri,rj
+c     do linear interpolation
+      ri=amax1(1.,amin1(float(n1),rid))
+      rj=amax1(1.,amin1(float(n2),rjd))
+      ih=nint(ri)
+      jh=nint(rj)
+c     Check for interpolation in i
+*     if (abs(float(ih)-ri).lt.1.e-3) then
+*       i  =ih
+*       ip1=ih
+*     else
+        i =min0(int(ri),n1-1)
+        ip1=i+1
+*     endif
+c     Check for interpolation in j
+*     if (abs(float(jh)-rj).lt.1.e-3) then
+*       j  =jh
+*       jp1=jh
+*     else
+        j =min0(int(rj),n2-1)
+        jp1=j+1
+*     endif
+      if ((i.eq.ip1).and.(j.eq.jp1)) then
+c        no interpolation at all
+         int2d=ar(i,j)
+      else
+         frac0i=ri-float(i)
+         frac0j=rj-float(j)
+         frac1i=1.-frac0i
+         frac1j=1.-frac0j
+         int2d = ar(i  ,j  ) * frac1i * frac1j
+     &         + ar(i  ,jp1) * frac1i * frac0j
+     &         + ar(ip1,j  ) * frac0i * frac1j
+     &         + ar(ip1,jp1) * frac0i * frac0j
+      endif
+      end
+      real function int2dm(ar,n1,n2,rid,rjd,misdat)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This subroutine interpolates a 2d-array to an arbitrary
+c        location within the grid. The interpolation includes the
+c        testing of the missing data flag 'misdat'.
+c     Arguments:
+c        ar        real  input   surface pressure, define as ar(n1,n2)
+c        n1,n2     int   input   dimensions of ar
+c        ri,rj     real  input   grid location to be interpolated to
+c        misdat    real  input   missing data flag (on if misdat<>0)
+c     Warning:
+c        This routine has not yet been seriously tested
+c     History:
+c-----------------------------------------------------------------------
+c     argument declarations
+      integer   n1,n2
+      real      ar(n1,n2), rid,rjd, misdat
+c     local declarations
+      integer   i,j,ip1,jp1,ih,jh
+      real      frac0i,frac0j,frac1i,frac1j,ri,rj,int2d
+c     check if routine without missing data checking can be called instead
+      if (misdat.eq.0.) then
+        int2dm=int2d(ar,n1,n2,rid,rjd)
+        return
+      endif
+c     do linear interpolation
+      ri=amax1(1.,amin1(float(n1),rid))
+      rj=amax1(1.,amin1(float(n2),rjd))
+      ih=nint(ri)
+      jh=nint(rj)
+c     Check for interpolation in i
+*     if (abs(float(ih)-ri).lt.1.e-3) then
+*       i  =ih
+*       ip1=ih
+*     else
+        i =min0(int(ri),n1-1)
+        ip1=i+1
+*     endif
+c     Check for interpolation in j
+*     if (abs(float(jh)-rj).lt.1.e-3) then
+*       j  =jh
+*       jp1=jh
+*     else
+        j =min0(int(rj),n2-1)
+        jp1=j+1
+*     endif
+      if ((i.eq.ip1).and.(j.eq.jp1)) then
+c        no interpolation at all
+         int2dm=ar(i,j)
+      else
+         if ((misdat.eq.ar(i  ,j  )).or.
+     &       (misdat.eq.ar(i  ,jp1)).or.
+     &       (misdat.eq.ar(ip1,j  )).or.
+     &       (misdat.eq.ar(ip1,jp1))) then
+           int2dm=misdat
+         else
+           frac0i=ri-float(i)
+           frac0j=rj-float(j)
+           frac1i=1.-frac0i
+           frac1j=1.-frac0j
+           int2dm = ar(i  ,j  ) * frac1i * frac1j
+     &            + ar(i  ,jp1) * frac1i * frac0j
+     &            + ar(ip1,j  ) * frac0i * frac1j
+     &            + ar(ip1,jp1) * frac0i * frac0j
+         endif
+      endif
+      end
+      real function int2dp(ar,n1,n2,rid,rjd)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This subroutine interpolates a 2d-array to an arbitrary
+c        location within the grid. The 2d-array is periodic in the
+c        i-direction: ar(0,.)=ar(n1,.) and ar(1,.)=ar(n1+1,.).
+c        Therefore rid can take values in the range 0,...,n1+1
+c     Arguments:
+c        ar        real  input   surface pressure, define as ar(n1,n2)
+c        n1,n2     int   input   dimensions of ar
+c        ri,rj     real  input   grid location to be interpolated to
+c     History:
+c-----------------------------------------------------------------------
+c     argument declarations
+      integer   n1,n2
+      real      ar(0:n1+1,n2), rid,rjd
+c     local declarations
+      integer   i,j,ip1,jp1,ih,jh
+      real      frac0i,frac0j,frac1i,frac1j,ri,rj
+c     do linear interpolation
+      ri=amax1(0.,amin1(float(n1+1),rid))
+      rj=amax1(1.,amin1(float(n2),rjd))
+      ih=nint(ri)
+      jh=nint(rj)
+c     Check for interpolation in i
+*     if (abs(float(ih)-ri).lt.1.e-5) then
+*       i  =ih
+*       ip1=ih
+*     else
+        i =min0(int(ri),n1-1)
+        ip1=i+1
+*     endif
+c     Check for interpolation in j
+*     if (abs(float(jh)-rj).lt.1.e-5) then
+*       j  =jh
+*       jp1=jh
+*     else
+        j =min0(int(rj),n2-1)
+        jp1=j+1
+*     endif
+      if ((i.eq.ip1).and.(j.eq.jp1)) then
+c        no interpolation at all
+         int2dp=ar(i,j)
+      else
+         frac0i=ri-float(i)
+         frac0j=rj-float(j)
+         frac1i=1.-frac0i
+         frac1j=1.-frac0j
+         int2dp = ar(i  ,j  ) * frac1i * frac1j
+     &         + ar(i  ,jp1) * frac1i * frac0j
+     &         + ar(ip1,j  ) * frac0i * frac1j
+     &         + ar(ip1,jp1) * frac0i * frac0j
+      endif
+      end
+      real function cint2d(ar,n1,n2,n3,rid,rjd,ikd)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This subroutine interpolates a 3d-array to an arbitrary
+c        location in the horizontal. ikd specifies the level (must
+c        be an integer). A bicubic method is applied (following
+c        the Numerical Recipes).
+c     Arguments:
+c        ar           real  input   field, define as ar(n1,n2,n3)
+c        n1,n2        int   input   dimensions of ar
+c        rid,rjd      real  input   grid location to be interpolated to
+c        ikd          int   input   level for interpolation
+c     History:
+c-----------------------------------------------------------------------
+c     argument declarations
+      integer   n1,n2,n3,ikd
+      real      ar(n1,n2,n3), rid,rjd
+c     local declarations
+      integer   i,j,k
+      real      y(4),y1(4),y2(4),y12(4)
+c     indices of lower left corner of interpolation grid
+      i=int(rid)
+      j=int(rjd)
+      k=ikd
+c     do not interpolate if i or j are at the boundaries
+      if ((i.eq.1).or.(j.eq.1).or.(i.eq.n1).or.(j.eq.n2)) then
+        cint2d=ar(i,j,k)
+        return
+      endif
+c     define the arrays y,y1,y2,y12 necessary for the bicubic
+c     interpolation (following the Numerical Recipes).
+      y(1)=ar(i,j,k)
+      y(2)=ar(i+1,j,k)
+      y(3)=ar(i+1,j+1,k)
+      y(4)=ar(i,j+1,k)
+      y1(1)=-(ar(i-1,j,k)-ar(i+1,j,k))/2.
+      y1(2)=-(ar(i,j,k)-ar(i+2,j,k))/2.
+      y1(3)=-(ar(i,j+1,k)-ar(i+2,j+1,k))/2.
+      y1(4)=-(ar(i-1,j+1,k)-ar(i+1,j+1,k))/2.
+      y2(1)=-(ar(i,j-1,k)-ar(i,j+1,k))/2.
+      y2(2)=-(ar(i+1,j-1,k)-ar(i+1,j+1,k))/2.
+      y2(3)=-(ar(i+1,j,k)-ar(i+1,j+2,k))/2.
+      y2(4)=-(ar(i,j,k)-ar(i,j+2,k))/2.
+      y12(1)=(ar(i+1,j+1,k)-ar(i+1,j-1,k)-ar(i-1,j+1,k)+ar(i-1,j-1,k))/4.
+      y12(2)=(ar(i+2,j+1,k)-ar(i+2,j-1,k)-ar(i,j+1,k)+ar(i,j-1,k))/4.
+      y12(3)=(ar(i+2,j+2,k)-ar(i+2,j,k)-ar(i,j+2,k)+ar(i,j,k))/4.
+      y12(4)=(ar(i+1,j+2,k)-ar(i+1,j,k)-ar(i-1,j+2,k)+ar(i-1,j,k))/4.
+      call bcuint(y,y1,y2,y12,i,j,rid,rjd,cint2d)
+      return
+      end
+      real function int3d(ar,n1,n2,n3,rid,rjd,rkd)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This subroutine interpolates a 3d-array to an arbitrary
+c        location within the grid.
+c     Arguments:
+c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
+c        n1,n2,n3  int   input   dimensions of ar
+c        ri,rj,rk  real  input   grid location to be interpolated to
+c     History:
+c-----------------------------------------------------------------------
+c     argument declarations
+      integer   n1,n2,n3
+      real      ar(n1,n2,n3), rid,rjd,rkd
+c     local declarations
+      integer   i,j,k,ip1,jp1,kp1,ih,jh,kh
+      real      frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk
+c     do linear interpolation
+      ri=amax1(1.,amin1(float(n1),rid))
+      rj=amax1(1.,amin1(float(n2),rjd))
+      rk=amax1(1.,amin1(float(n3),rkd))
+      ih=nint(ri)
+      jh=nint(rj)
+      kh=nint(rk)
+c     Check for interpolation in i
+*     if (abs(float(ih)-ri).lt.1.e-3) then
+*       i  =ih
+*       ip1=ih
+*     else
+        i =min0(int(ri),n1-1)
+        ip1=i+1
+*     endif
+c     Check for interpolation in j
+      if (abs(float(jh)-rj).lt.1.e-3) then
+        j  =jh
+        jp1=jh
+      else
+        j =min0(int(rj),n2-1)
+        jp1=j+1
+      endif
+c     Check for interpolation in k
+*     if (abs(float(kh)-rk).lt.1.e-3) then
+*       k  =kh
+*       kp1=kh
+*     else
+        k =min0(int(rk),n3-1)
+        kp1=k+1
+*     endif
+      if (k.eq.kp1) then
+c       no interpolation in k
+        if ((i.eq.ip1).and.(j.eq.jp1)) then
+c          no interpolation at all
+           int3d=ar(i,j,k)
+c          print *,'int3d 00: ',rid,rjd,rkd,int3d
+        else
+c          horizontal interpolation only
+           frac0i=ri-float(i)
+           frac0j=rj-float(j)
+           frac1i=1.-frac0i
+           frac1j=1.-frac0j
+           int3d = ar(i  ,j  ,k  ) * frac1i * frac1j
+     &           + ar(i  ,jp1,k  ) * frac1i * frac0j
+     &           + ar(ip1,j  ,k  ) * frac0i * frac1j
+     &           + ar(ip1,jp1,k  ) * frac0i * frac0j
+c          print *,'int3d 10: ',rid,rjd,rkd,int3d
+        endif
+      else
+        frac0k=rk-float(k)
+        frac1k=1.-frac0k
+        if ((i.eq.ip1).and.(j.eq.jp1)) then
+c          vertical interpolation only
+           int3d = ar(i  ,j  ,k  ) * frac1k
+     &           + ar(i  ,j  ,kp1) * frac0k
+c          print *,'int3d 01: ',rid,rjd,rkd,int3d
+        else
+c          full 3d interpolation
+           frac0i=ri-float(i)
+           frac0j=rj-float(j)
+           frac1i=1.-frac0i
+           frac1j=1.-frac0j
+           int3d = ar(i  ,j  ,k  ) * frac1i * frac1j * frac1k
+     &           + ar(i  ,jp1,k  ) * frac1i * frac0j * frac1k
+     &           + ar(ip1,j  ,k  ) * frac0i * frac1j * frac1k
+     &           + ar(ip1,jp1,k  ) * frac0i * frac0j * frac1k
+     &           + ar(i  ,j  ,kp1) * frac1i * frac1j * frac0k
+     &           + ar(i  ,jp1,kp1) * frac1i * frac0j * frac0k
+     &           + ar(ip1,j  ,kp1) * frac0i * frac1j * frac0k
+     &           + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
+c          print *,'int3d 11: ',rid,rjd,rkd,int3d
+        endif
+      endif
+      end
+      real function int3dm(ar,n1,n2,n3,rid,rjd,rkd,misdat)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This subroutine interpolates a 3d-array to an arbitrary
+c        location within the grid. The interpolation includes the
+c        testing of the missing data flag 'misdat'.
+c     Arguments:
+c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
+c        n1,n2,n3  int   input   dimensions of ar
+c        ri,rj,rk  real  input   grid location to be interpolated to
+c        misdat    real  input   missing data flag (on if misdat<>0)
+c     Warning:
+c        This routine has not yet been seriously tested
+c     History:
+c-----------------------------------------------------------------------
+c     argument declarations
+      integer   n1,n2,n3
+      real      ar(n1,n2,n3), rid,rjd,rkd, misdat
+c     local declarations
+      integer   i,j,k,ip1,jp1,kp1,ih,jh,kh
+      real      frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk,int3d
+c     check if routine without missing data checking can be called instead
+      if (misdat.eq.0.) then
+        int3dm=int3d(ar,n1,n2,n3,rid,rjd,rkd)
+        return
+      endif
+c     do linear interpolation
+      ri=amax1(1.,amin1(float(n1),rid))
+      rj=amax1(1.,amin1(float(n2),rjd))
+      rk=amax1(1.,amin1(float(n3),rkd))
+      ih=nint(ri)
+      jh=nint(rj)
+      kh=nint(rk)
+c     Check for interpolation in i
+*     if (abs(float(ih)-ri).lt.1.e-3) then
+*       i  =ih
+*       ip1=ih
+*     else
+        i =min0(int(ri),n1-1)
+        ip1=i+1
+*     endif
+c     Check for interpolation in j
+*     if (abs(float(jh)-rj).lt.1.e-3) then
+*       j  =jh
+*       jp1=jh
+*     else
+        j =min0(int(rj),n2-1)
+        jp1=j+1
+*     endif
+c     Check for interpolation in k
+*     if (abs(float(kh)-rk).lt.1.e-3) then
+*       k  =kh
+*       kp1=kh
+*     else
+        k =min0(int(rk),n3-1)
+        kp1=k+1
+*     endif
+      if (k.eq.kp1) then
+c       no interpolation in k
+        if ((i.eq.ip1).and.(j.eq.jp1)) then
+c          no interpolation at all
+           if (misdat.eq.ar(i,j,k)) then
+             int3dm=misdat
+           else
+             int3dm=ar(i,j,k)
+           endif
+c          print *,'int3dm 00: ',rid,rjd,rkd,int3dm
+        else
+c          horizontal interpolation only
+           if ((misdat.eq.ar(i  ,j  ,k  )).or.
+     &         (misdat.eq.ar(i  ,jp1,k  )).or.
+     &         (misdat.eq.ar(ip1,j  ,k  )).or.
+     &         (misdat.eq.ar(ip1,jp1,k  ))) then
+             int3dm=misdat
+           else
+             frac0i=ri-float(i)
+             frac0j=rj-float(j)
+             frac1i=1.-frac0i
+             frac1j=1.-frac0j
+             int3dm = ar(i  ,j  ,k  ) * frac1i * frac1j
+     &              + ar(i  ,jp1,k  ) * frac1i * frac0j
+     &              + ar(ip1,j  ,k  ) * frac0i * frac1j
+     &              + ar(ip1,jp1,k  ) * frac0i * frac0j
+c            print *,'int3dm 10: ',rid,rjd,rkd,int3dm
+           endif
+        endif
+      else
+        frac0k=rk-float(k)
+        frac1k=1.-frac0k
+        if ((i.eq.ip1).and.(j.eq.jp1)) then
+c          vertical interpolation only
+           if ((misdat.eq.ar(i  ,j  ,k  )).or.
+     &         (misdat.eq.ar(i  ,j  ,kp1))) then
+             int3dm=misdat
+           else
+             int3dm = ar(i  ,j  ,k  ) * frac1k
+     &              + ar(i  ,j  ,kp1) * frac0k
+c            print *,'int3dm 01: ',rid,rjd,rkd,int3dm
+           endif
+        else
+c          full 3d interpolation
+           if ((misdat.eq.ar(i  ,j  ,k  )).or.
+     &         (misdat.eq.ar(i  ,jp1,k  )).or.
+     &         (misdat.eq.ar(ip1,j  ,k  )).or.
+     &         (misdat.eq.ar(ip1,jp1,k  )).or.
+     &         (misdat.eq.ar(i  ,j  ,kp1)).or.
+     &         (misdat.eq.ar(i  ,jp1,kp1)).or.
+     &         (misdat.eq.ar(ip1,j  ,kp1)).or.
+     &         (misdat.eq.ar(ip1,jp1,kp1))) then
+             int3dm=misdat
+           else
+             frac0i=ri-float(i)
+             frac0j=rj-float(j)
+             frac1i=1.-frac0i
+             frac1j=1.-frac0j
+             int3dm = ar(i  ,j  ,k  ) * frac1i * frac1j * frac1k
+     &              + ar(i  ,jp1,k  ) * frac1i * frac0j * frac1k
+     &              + ar(ip1,j  ,k  ) * frac0i * frac1j * frac1k
+     &              + ar(ip1,jp1,k  ) * frac0i * frac0j * frac1k
+     &              + ar(i  ,j  ,kp1) * frac1i * frac1j * frac0k
+     &              + ar(i  ,jp1,kp1) * frac1i * frac0j * frac0k
+     &              + ar(ip1,j  ,kp1) * frac0i * frac1j * frac0k
+     &              + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
+c            print *,'int3dm 11: ',rid,rjd,rkd,int3dm
+           endif
+        endif
+      endif
+      end
+      real function int3dmlog(ar,n1,n2,n3,rid,rjd,rkd,misdat)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This subroutine interpolates a 3d-array to an arbitrary
+c        location within the grid. The interpolation includes the
+c        testing of the missing data flag 'misdat'.
+c        Prior to vertical interpolations the log is taken from the array.
+c     Arguments:
+c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
+c        n1,n2,n3  int   input   dimensions of ar
+c        ri,rj,rk  real  input   grid location to be interpolated to
+c        misdat    real  input   missing data flag (on if misdat<>0)
+c     Warning:
+c        This routine has not yet been seriously tested
+c     History:
+c-----------------------------------------------------------------------
+c     argument declarations
+      integer   n1,n2,n3
+      real      ar(n1,n2,n3), rid,rjd,rkd, misdat
+c     local declarations
+      integer   i,j,k,ip1,jp1,kp1,ih,jh,kh
+      real      frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk,int3d
+c     print*,'hallo in SR int3dmlog'
+c     check if routine without missing data checking can be called instead
+      if (misdat.eq.0.) then
+        int3dmlog=int3d(ar,n1,n2,n3,rid,rjd,rkd)
+        return
+      endif
+c     do linear interpolation
+      ri=amax1(1.,amin1(float(n1),rid))
+      rj=amax1(1.,amin1(float(n2),rjd))
+      rk=amax1(1.,amin1(float(n3),rkd))
+      ih=nint(ri)
+      jh=nint(rj)
+      kh=nint(rk)
+c     Check for interpolation in i
+*     if (abs(float(ih)-ri).lt.1.e-3) then
+*       i  =ih
+*       ip1=ih
+*     else
+        i =min0(int(ri),n1-1)
+        ip1=i+1
+*     endif
+c     Check for interpolation in j
+*     if (abs(float(jh)-rj).lt.1.e-3) then
+*       j  =jh
+*       jp1=jh
+*     else
+        j =min0(int(rj),n2-1)
+        jp1=j+1
+*     endif
+c     Check for interpolation in k
+*     if (abs(float(kh)-rk).lt.1.e-3) then
+*       k  =kh
+*       kp1=kh
+*     else
+        k =min0(int(rk),n3-1)
+        kp1=k+1
+*     endif
+      if (k.eq.kp1) then
+c       no interpolation in k
+        if ((i.eq.ip1).and.(j.eq.jp1)) then
+c          no interpolation at all
+           if (misdat.eq.ar(i,j,k)) then
+             int3dmlog=misdat
+           else
+             int3dmlog=ar(i,j,k)
+           endif
+c          print *,'int3dmlog 00: ',rid,rjd,rkd,int3dmlog
+        else
+c          horizontal interpolation only
+           if ((misdat.eq.ar(i  ,j  ,k  )).or.
+     &         (misdat.eq.ar(i  ,jp1,k  )).or.
+     &         (misdat.eq.ar(ip1,j  ,k  )).or.
+     &         (misdat.eq.ar(ip1,jp1,k  ))) then
+             int3dmlog=misdat
+           else
+             frac0i=ri-float(i)
+             frac0j=rj-float(j)
+             frac1i=1.-frac0i
+             frac1j=1.-frac0j
+             int3dmlog = ar(i  ,j  ,k  ) * frac1i * frac1j
+     &              + ar(i  ,jp1,k  ) * frac1i * frac0j
+     &              + ar(ip1,j  ,k  ) * frac0i * frac1j
+     &              + ar(ip1,jp1,k  ) * frac0i * frac0j
+c            print *,'int3dmlog 10: ',rid,rjd,rkd,int3dmlog
+           endif
+        endif
+      else
+        frac0k=rk-float(k)
+        frac1k=1.-frac0k
+        if ((i.eq.ip1).and.(j.eq.jp1)) then
+c          vertical interpolation only
+           if ((misdat.eq.ar(i  ,j  ,k  )).or.
+     &         (misdat.eq.ar(i  ,j  ,kp1))) then
+             int3dmlog=misdat
+           else
+             int3dmlog = log(ar(i  ,j  ,k  )) * frac1k
+     &                 + log(ar(i  ,j  ,kp1)) * frac0k
+             int3dmlog = exp(int3dmlog)
+c            print *,'int3dmlog 01: ',rid,rjd,rkd,int3dmlog
+           endif
+        else
+c          full 3d interpolation
+           if ((misdat.eq.ar(i  ,j  ,k  )).or.
+     &         (misdat.eq.ar(i  ,jp1,k  )).or.
+     &         (misdat.eq.ar(ip1,j  ,k  )).or.
+     &         (misdat.eq.ar(ip1,jp1,k  )).or.
+     &         (misdat.eq.ar(i  ,j  ,kp1)).or.
+     &         (misdat.eq.ar(i  ,jp1,kp1)).or.
+     &         (misdat.eq.ar(ip1,j  ,kp1)).or.
+     &         (misdat.eq.ar(ip1,jp1,kp1))) then
+             int3dmlog=misdat
+           else
+             frac0i=ri-float(i)
+             frac0j=rj-float(j)
+             frac1i=1.-frac0i
+             frac1j=1.-frac0j
+             int3dmlog = log(ar(i  ,j  ,k  )) * frac1i * frac1j * frac1k
+     &                 + log(ar(i  ,jp1,k  )) * frac1i * frac0j * frac1k
+     &                 + log(ar(ip1,j  ,k  )) * frac0i * frac1j * frac1k
+     &                 + log(ar(ip1,jp1,k  )) * frac0i * frac0j * frac1k
+     &                 + log(ar(i  ,j  ,kp1)) * frac1i * frac1j * frac0k
+     &                 + log(ar(i  ,jp1,kp1)) * frac1i * frac0j * frac0k
+     &                 + log(ar(ip1,j  ,kp1)) * frac0i * frac1j * frac0k
+     &                 + log(ar(ip1,jp1,kp1)) * frac0i * frac0j * frac0k
+             int3dmlog = exp(int3dmlog)
+c            print *,'int3dmlog 11: ',rid,rjd,rkd,int3dmlog
+           endif
+        endif
+      endif
+      end
+      real function int3dmvc(ar,n1,n2,n3,rid,rjd,rkd,misdat)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This subroutine interpolates a 3d-array to an arbitrary
+c        location within the grid. The interpolation includes the
+c        testing of the missing data flag 'misdat'.
+c        In the vertical a Lagrangian cubic interpolation is
+c        performed.
+c     Arguments:
+c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
+c        n1,n2,n3  int   input   dimensions of ar
+c        ri,rj,rk  real  input   grid location to be interpolated to
+c        misdat    real  input   missing data flag (on if misdat<>0)
+c     Warning:
+c        This routine has not yet been seriously tested
+c     History:
+c-----------------------------------------------------------------------
+c     argument declarations
+      integer   n1,n2,n3
+      real      ar(n1,n2,n3), rid,rjd,rkd, misdat
+c     local declarations
+      integer   i,j,k,ip1,jp1,kp1,ih,jh,kh,klow,n
+      real      frac0i,frac0j,frac1i,frac1j,ri,rj,rk
+      real      int2d(4)
+      real      int3dm
+c     if n3 < 4 then do call linear interpolation in the vertical
+      if (n3.lt.4) then
+        int3dmvc=int3dm(ar,n1,n2,n3,rid,rjd,rkd,misdat)
+        return
+      endif
+c     do linear interpolation in the horizontal
+      ri=amax1(1.,amin1(float(n1),rid))
+      rj=amax1(1.,amin1(float(n2),rjd))
+      rk=amax1(1.,amin1(float(n3),rkd))
+      ih=nint(ri)
+      jh=nint(rj)
+      kh=nint(rk)
+c     Check for interpolation in i
+*     if (abs(float(ih)-ri).lt.1.e-3) then
+*       i  =ih
+*       ip1=ih
+*     else
+        i =min0(int(ri),n1-1)
+        ip1=i+1
+*     endif
+c     Check for interpolation in j
+*     if (abs(float(jh)-rj).lt.1.e-3) then
+*       j  =jh
+*       jp1=jh
+*     else
+        j =min0(int(rj),n2-1)
+        jp1=j+1
+*     endif
+c     Check for interpolation in k
+*     if (abs(float(kh)-rk).lt.1.e-3) then
+*       k  =kh
+*       kp1=kh
+*     else
+        k =min0(int(rk),n3-1)
+        kp1=k+1
+*     endif
+      if (k.eq.kp1) then
+c       no interpolation in k
+        if ((i.eq.ip1).and.(j.eq.jp1)) then
+c          no interpolation at all
+           int3dmvc=ar(i,j,k)
+c          print *,'int3dmvc 00: ',rid,rjd,rkd,int3dmvc
+        else
+c          horizontal interpolation only
+           if ((misdat.eq.ar(i  ,j  ,k  )).or.
+     &         (misdat.eq.ar(i  ,jp1,k  )).or.
+     &         (misdat.eq.ar(ip1,j  ,k  )).or.
+     &         (misdat.eq.ar(ip1,jp1,k  ))) then
+             int3dmvc=misdat
+           else
+             frac0i=ri-float(i)
+             frac0j=rj-float(j)
+             frac1i=1.-frac0i
+             frac1j=1.-frac0j
+             int3dmvc = ar(i  ,j  ,k  ) * frac1i * frac1j
+     &              + ar(i  ,jp1,k  ) * frac1i * frac0j
+     &              + ar(ip1,j  ,k  ) * frac0i * frac1j
+     &              + ar(ip1,jp1,k  ) * frac0i * frac0j
+c            print *,'int3dmvc 10: ',rid,rjd,rkd,int3dmvc
+           endif
+        endif
+      else
+        if ((i.eq.ip1).and.(j.eq.jp1)) then
+c          vertical interpolation only
+           if ((misdat.eq.ar(i  ,j  ,k  )).or.
+     &         (misdat.eq.ar(i  ,j  ,kp1))) then
+             int3dmvc=misdat
+           else
+             if (k-1.lt.1) then
+               klow=1
+             else if (k+2.gt.n3) then
+               klow=n3-3
+             else
+               klow=k-1
+             endif
+             call cubint(ar(i,j,klow),ar(i,j,klow+1),ar(i,j,klow+2),
+     &                   ar(i,j,klow+3),klow,rk,int3dmvc)
+           endif
+        else
+c          full 3d interpolation
+           if ((misdat.eq.ar(i  ,j  ,k  )).or.
+     &         (misdat.eq.ar(i  ,jp1,k  )).or.
+     &         (misdat.eq.ar(ip1,j  ,k  )).or.
+     &         (misdat.eq.ar(ip1,jp1,k  )).or.
+     &         (misdat.eq.ar(i  ,j  ,kp1)).or.
+     &         (misdat.eq.ar(i  ,jp1,kp1)).or.
+     &         (misdat.eq.ar(ip1,j  ,kp1)).or.
+     &         (misdat.eq.ar(ip1,jp1,kp1))) then
+             int3dmvc=misdat
+           else
+             frac0i=ri-float(i)
+             frac0j=rj-float(j)
+             frac1i=1.-frac0i
+             frac1j=1.-frac0j
+             if (k-1.lt.1) then
+               klow=1
+             else if (k+2.gt.n3) then
+               klow=n3-3
+             else
+               klow=k-1
+             endif
+             do n=1,4
+               int2d(n) = ar(i  ,j  ,klow-1+n  ) * frac1i * frac1j
+     &                  + ar(i  ,jp1,klow-1+n  ) * frac1i * frac0j
+     &                  + ar(ip1,j  ,klow-1+n  ) * frac0i * frac1j
+     &                  + ar(ip1,jp1,klow-1+n  ) * frac0i * frac0j
+             enddo
+             call cubint(int2d(1),int2d(2),int2d(3),int2d(4),
+     &                   klow,rk,int3dmvc)
+           endif
+        endif
+      endif
+      end
+      real function int4d(ar,n1,n2,n3,n4,rid,rjd,rkd,rld)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This subroutine interpolates a 4d-array to an arbitrary
+c        location within the grid.
+c     Arguments:
+c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
+c        n1,..,n4  int   input   dimensions of ar
+c        ri,..,rl  real  input   grid location to be interpolated to
+c     History:
+c-----------------------------------------------------------------------
+c     argument declarations
+      integer   n1,n2,n3,n4
+      real      ar(n1,n2,n3,n4), rid,rjd,rkd,rld
+c     local declarations
+      integer   l,lp1,lh
+      real      frac0l,frac1l,rl,int3d
+c     do linear interpolation in l-direction
+      rl=amax1(1.,amin1(float(n4),rld))
+      lh=nint(rl)
+c     Check for interpolation in l
+*     if (abs(float(lh)-rl).lt.1.e-3) then
+*       l  =lh
+*       lp1=lh
+*     else
+        l =min0(int(rl),n4-1)
+        lp1=l+1
+*     endif
+      if (l.eq.lp1) then
+c       no interpolation in l
+        int4d=int3d(ar(1,1,1,l),n1,n2,n3,rid,rjd,rkd)
+      else
+c       interpolation in l
+        frac0l=rl-float(l)
+        frac1l=1.-frac0l
+        int4d = int3d(ar(1,1,1,l  ),n1,n2,n3,rid,rjd,rkd) * frac1l
+     &        + int3d(ar(1,1,1,lp1),n1,n2,n3,rid,rjd,rkd) * frac0l
+      endif
+      end
+      real function int4dm(ar,n1,n2,n3,n4,rid,rjd,rkd,rld,misdat)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This subroutine interpolates a 4d-array to an arbitrary
+c        location within the grid. The interpolation includes the
+c        testing of the missing data flag 'misdat'.
+c     Arguments:
+c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
+c        n1,..,n4  int   input   dimensions of ar
+c        ri,..,rl  real  input   grid location to be interpolated to
+c        misdat    real  input   missing data flag (on if misdat<>0)
+c     Warning:
+c        This routine has not yet been seriously tested.
+c     History:
+c-----------------------------------------------------------------------
+c     argument declarations
+      integer   n1,n2,n3,n4
+      real      ar(n1,n2,n3,n4), rid,rjd,rkd,rld, misdat
+c     local declarations
+      integer   l,lp1,lh
+      real      frac0l,frac1l,rl,rint0,rint1,int4d,int3dm
+c     check whether missing data checking is required
+      if (misdat.eq.0.) then
+        int4dm=int4d(ar,n1,n2,n3,n4,rid,rjd,rkd,rld)
+        return
+      endif
+c     do linear interpolation in l-direction
+      rl=amax1(1.,amin1(float(n4),rld))
+      lh=nint(rl)
+c     Check for interpolation in l
+*     if (abs(float(lh)-rl).lt.1.e-3) then
+*       l  =lh
+*       lp1=lh
+*     else
+        l =min0(int(rl),n4-1)
+        lp1=l+1
+*     endif
+      if (l.eq.lp1) then
+c       no interpolation in l
+        int4dm = int3dm(ar(1,1,1,l),n1,n2,n3,rid,rjd,rkd,misdat)
+      else
+c       interpolation in l
+        frac0l=rl-float(l)
+        frac1l=1.-frac0l
+        rint0 = int3dm(ar(1,1,1,l  ),n1,n2,n3,rid,rjd,rkd,misdat)
+        rint1 = int3dm(ar(1,1,1,lp1),n1,n2,n3,rid,rjd,rkd,misdat)
+        if ((rint0.eq.misdat).or.(rint1.eq.misdat)) then
+          int4dm = misdat
+        else
+          int4dm = rint0*frac1l + rint1*frac0l
+        endif
+      endif
+      end
+      real function int3dl(ar,n1,n2,n3,levels,rid,rjd,rkd)
+c-----------------------------------------------------------------------
+c     Purpose:
+c        This subroutine interpolates a 3d-array to an arbitrary
+c        location within the grid. The vertical interpolation is linear
+c        in log(pressure).
+c     Arguments:
+c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
+c        n1,n2,n3  int   input   dimensions of ar
+c        levels    real  input   array contains pressure levels for ar
+c        ri,rj,rk  real  input   grid location to be interpolated to
+c     History:
+c        Based on int3d                 July 93
+c-----------------------------------------------------------------------
+c     argument declarations
+      integer   n1,n2,n3
+      real      ar(n1,n2,n3), rid,rjd,rkd
+      real      levels(n3)
+c     local declarations
+      real      pval
+      integer   i,j,k,ip1,jp1,kp1,ih,jh,kh
+      real      frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk
+c     do linear interpolation
+      ri=amax1(1.,amin1(float(n1),rid))
+      rj=amax1(1.,amin1(float(n2),rjd))
+      rk=amax1(1.,amin1(float(n3),rkd))
+      ih=nint(ri)
+      jh=nint(rj)
+      kh=nint(rk)
+c     Check for interpolation in i
+*     if (abs(float(ih)-ri).lt.1.e-3) then
+*       i  =ih
+*       ip1=ih
+*     else
+        i =min0(int(ri),n1-1)
+        ip1=i+1
+*     endif
+c     Check for interpolation in j
+*     if (abs(float(jh)-rj).lt.1.e-3) then
+*       j  =jh
+*       jp1=jh
+*     else
+        j =min0(int(rj),n2-1)
+        jp1=j+1
+*     endif
+c     Check for interpolation in k
+*     if (abs(float(kh)-rk).lt.1.e-3) then
+*       k  =kh
+*       kp1=kh
+*     else
+        k =min0(int(rk),n3-1)
+        kp1=k+1
+*     endif
+      if (k.eq.kp1) then
+c       no interpolation in k
+        if ((i.eq.ip1).and.(j.eq.jp1)) then
+c          no interpolation at all
+           int3dl=ar(i,j,k)
+c          print *,'int3dl 00: ',rid,rjd,rkd,int3dl
+        else
+c          horizontal interpolation only
+           frac0i=ri-float(i)
+           frac0j=rj-float(j)
+           frac1i=1.-frac0i
+           frac1j=1.-frac0j
+           int3dl = ar(i  ,j  ,k  ) * frac1i * frac1j
+     &            + ar(i  ,jp1,k  ) * frac1i * frac0j
+     &            + ar(ip1,j  ,k  ) * frac0i * frac1j
+     &            + ar(ip1,jp1,k  ) * frac0i * frac0j
+c          print *,'int3dl 10: ',rid,rjd,rkd,int3dl
+        endif
+      else
+*       frac0k=rk-float(k)
+c       calculate the pressure value to be interpolated to
+        pval=levels(int(rk))
+     >            -(rk-aint(rk))*(levels(int(rk))-levels(int(rk)+1))
+        frac0k=log(levels(int(rk))/pval)
+     &        /log(levels(int(rk))/levels(int(rk)+1))
+        frac1k=1.-frac0k
+        if ((i.eq.ip1).and.(j.eq.jp1)) then
+c          vertical interpolation only
+           int3dl = ar(i  ,j  ,k  ) * frac1k
+     &            + ar(i  ,j  ,kp1) * frac0k
+c          print *,'int3dl 01: ',rid,rjd,rkd,int3dl
+        else
+c          full 3d interpolation
+           frac0i=ri-float(i)
+           frac0j=rj-float(j)
+           frac1i=1.-frac0i
+           frac1j=1.-frac0j
+           int3dl = ar(i  ,j  ,k  ) * frac1i * frac1j * frac1k
+     &            + ar(i  ,jp1,k  ) * frac1i * frac0j * frac1k
+     &            + ar(ip1,j  ,k  ) * frac0i * frac1j * frac1k
+     &            + ar(ip1,jp1,k  ) * frac0i * frac0j * frac1k
+     &            + ar(i  ,j  ,kp1) * frac1i * frac1j * frac0k
+     &            + ar(i  ,jp1,kp1) * frac1i * frac0j * frac0k
+     &            + ar(ip1,j  ,kp1) * frac0i * frac1j * frac0k
+     &            + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
+c          print *,'int3dl 11: ',rid,rjd,rkd,int3dl
+        endif
+      endif
+      end
+      subroutine bcucof(y,y1,y2,y12,c)
+c-----------------------------------------------------------------------
+c       Given arrays y,y1,y2 and y12, each of length 4, containing the
+c       function, gradients, and cross derivative at the four grid points
+c       of a rectangular grid cell (numbered counterclockwise from the
+c       lower left), and given d1 and d2, the length of the grid cell in
+c       the 1- and 2-directions, this routine returns the table c that is
+c       used by routine bcuint for biqubic interpolation.
+c     Source: Numerical Recipes, Fortran Version, p.99
+c-----------------------------------------------------------------------
+      real      c(4,4),y(4),y1(4),y2(4),y12(4),cl(16),x(16)
+      integer   wt(16,16)
+      data      wt/1,0,-3,2,4*0,-3,0,9,-6,2,0,-6,4,
+     >          8*0,3,0,-9,6,-2,0,6,-4,
+     >          10*0,9,-6,2*0,-6,4,2*0,3,-2,6*0,-9,6,2*0,6,-4,
+     >          4*0,1,0,-3,2,-2,0,6,-4,1,0,-3,2,8*0,-1,0,3,-2,1,0,-3,2,
+     >          10*0,-3,2,2*0,3,-2,6*0,3,-2,2*0,-6,4,2*0,3,-2,
+     >          0,1,-2,1,5*0,-3,6,-3,0,2,-4,2,9*0,3,-6,3,0,-2,4,-2,
+     >          10*0,-3,3,2*0,2,-2,2*0,-1,1,6*0,3,-3,2*0,-2,2,
+     >          5*0,1,-2,1,0,-2,4,-2,0,1,-2,1,9*0,-1,2,-1,0,1,-2,1,
+     >          10*0,1,-1,2*0,-1,1,6*0,-1,1,2*0,2,-2,2*0,-1,1/
+      real      xx
+      integer   i,j,k,l
+      do i=1,4                  ! pack a temporary vector x
+        x(i)=y(i)
+        x(i+4)=y1(i)
+        x(i+8)=y2(i)
+        x(i+12)=y12(i)
+      enddo
+      do i=1,16                 ! matrix multiply by the stored table
+        xx=0.
+        do k=1,16
+          xx=xx+wt(i,k)*x(k)
+        enddo
+        cl(i)=xx
+      enddo
+      l=0
+      do i=1,4                  ! unpack the result into the output table
+      do j=1,4
+        l=l+1
+        c(i,j)=cl(l)
+      enddo
+      enddo
+      return
+      end
+      subroutine bcuint(y,y1,y2,y12,x1l,x2l,x1,x2,ansy)
+c-----------------------------------------------------------------------
+c       Bicubic interpolation within a grid square. Input quantities are
+c       y,y1,y2,y12 (as described in bcucof); x1l and x1u, the lower and
+c       upper coordinates of the grid square in the 1-direction; x2l and
+c       x2u likewise for the 2-direction; and x1,x2, the coordinates of
+c       the desired point for the interpolation. The interplated function
+c       value is returned as ansy. This routine calls bcucof.
+c     Source: Numerical Recipes, Fortran Version, p.99/100
+c     !!! changed the proposed code !!!
+c-----------------------------------------------------------------------
+      real      y(4),y1(4),y2(4),y12(4),c(4,4)
+      real      ansy,x1,x2,t,u
+      integer   i,x1l,x2l
+      call bcucof(y,y1,y2,y12,c)
+      t=x1-real(x1l)
+      u=x2-real(x2l)
+      ansy=0.
+      do i=4,1,-1
+        ansy=t*ansy+((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)
+      enddo
+      return
+      end
+      subroutine cubint(ya1,ya2,ya3,ya4,k,x,y)
+c-----------------------------------------------------------------------
+c       Interface routine for SR polint for special case of cubic
+c       interpolation in index space, with xa=k,k+1,k+2,k+3
+c-----------------------------------------------------------------------
+      integer   k
+      real      ya1,ya2,ya3,ya4,x,y
+      integer   n
+      real      xa(4),ya(4),dy
+      do n=1,4
+        xa(1)=real(k)
+        xa(2)=real(k+1)
+        xa(3)=real(k+2)
+        xa(4)=real(k+3)
+        ya(1)=ya1
+        ya(2)=ya2
+        ya(3)=ya3
+        ya(4)=ya4
+      enddo
+      call polint(xa,ya,4,x,y,dy)
+      return
+      end
+      subroutine polint(xa,ya,n,x,y,dy)
+c-----------------------------------------------------------------------
+c       Given arrays xa and ya, each of length n, and given a value x, this
+c       routine returns a value y, and an error estimate dy. If P(x) is the
+c       polynomial of degree n-1 such that p(xa(i))=ya(i),i=1,...,n, then
+c       the returned value y=p(x)
+c     Source: Numerical Recipes, Fortran Version, p.82
+c-----------------------------------------------------------------------
+      integer   nmax,n
+      parameter(nmax=10)
+      real      xa(n),ya(n),x,y,dy
+      integer   i,m,ns
+      real      c(nmax),d(nmax),den,dif,dift,ho,hp,w
+      ns=1
+      dif=abs(x-xa(1))
+      do i=1,n
+        dift=abs(x-xa(i))
+        if (dift.lt.dif) then
+          ns=i
+          dif=dift
+        endif
+        c(i)=ya(i)
+        d(i)=ya(i)
+      enddo
+      y=ya(ns)
+      ns=ns-1
+      do m=1,n-1
+        do i=1,n-m
+          ho=xa(i)-x
+          hp=xa(i+m)-x
+          w=c(i+1)-d(i)
+          den=ho-hp
+          den=w/den
+          d(i)=hp*den
+          c(i)=ho*den
+        enddo
+        if (2*ns.lt.n-m) then
+          dy=c(ns+1)
+        else
+          dy=d(ns)
+          ns=ns-1
+        endif
+        y=y+dy
+      enddo
+      return
+      end
+      subroutine filt2d (a,af,f1,f2,nx,ny,fil,misdat,
+     &                                     iperx,ipery,ispol,inpol)
+c     =============================================================
+c     Apply a conservative diffusion operator onto the 2d field a,
+c     with full missing data checking.
+c
+c     a     real   inp  array to be filtered, dimensioned (nx,ny)
+c     af    real   out  filtered array, dimensioned (nx,ny), can be
+c                       equivalenced with array a in the calling routine
+c     f1    real        workarray, dimensioned (nx+1,ny)
+c     f2    real        workarray, dimensioned (nx,ny+1)
+c     fil   real   inp  filter-coeff., 0<afil<=1. Maximum filtering with afil=1
+c                       corresponds to one application of the linear filter.
+c     misdat real  inp  missing-data value, a(i,j)=misdat indicates that
+c                       the corresponding value is not available. The
+c                       misdat-checking can be switched off with with misdat=0.
+c     iperx int    inp  periodic boundaries in the x-direction (1=yes,0=no)
+c     ipery int    inp  periodic boundaries in the y-direction (1=yes,0=no)
+c     inpol int    inp  northpole at j=ny  (1=yes,0=no)
+c     ispol int    inp  southpole at j=1   (1=yes,0=no)
+c
+c     Christoph Schaer, 1993
+c     argument declaration
+      integer     nx,ny
+      real        a(nx,ny),af(nx,ny),f1(nx+1,ny),f2(nx,ny+1),fil,misdat
+      integer     iperx,ipery,inpol,ispol
+c     local variable declaration
+      integer     i,j,is
+      real        fh
+c     compute constant fh
+      fh=0.125*fil
+c     compute fluxes in x-direction
+      if (misdat.eq.0.) then
+        do j=1,ny
+        do i=2,nx
+          f1(i,j)=a(i-1,j)-a(i,j)
+        enddo
+        enddo
+      else
+        do j=1,ny
+        do i=2,nx
+          if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then
+            f1(i,j)=0.
+          else
+            f1(i,j)=a(i-1,j)-a(i,j)
+          endif
+        enddo
+        enddo
+      endif
+      if (iperx.eq.1) then
+c       do periodic boundaries in the x-direction
+        do j=1,ny
+          f1(1,j)=f1(nx,j)
+          f1(nx+1,j)=f1(2,j)
+        enddo
+      else
+c       set boundary-fluxes to zero
+        do j=1,ny
+          f1(1,j)=0.
+          f1(nx+1,j)=0.
+        enddo
+      endif
+c     compute fluxes in y-direction
+      if (misdat.eq.0.) then
+        do j=2,ny
+        do i=1,nx
+          f2(i,j)=a(i,j-1)-a(i,j)
+        enddo
+        enddo
+      else
+        do j=2,ny
+        do i=1,nx
+          if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then
+            f2(i,j)=0.
+          else
+            f2(i,j)=a(i,j-1)-a(i,j)
+          endif
+        enddo
+        enddo
+      endif
+c     set boundary-fluxes to zero
+      do i=1,nx
+        f2(i,1)=0.
+        f2(i,ny+1)=0.
+      enddo
+      if (ipery.eq.1) then
+c       do periodic boundaries in the x-direction
+        do i=1,nx
+          f2(i,1)=f2(i,ny)
+          f2(i,ny+1)=f2(i,2)
+        enddo
+      endif
+      if (iperx.eq.1) then
+        if (ispol.eq.1) then
+c         do south-pole
+          is=(nx-1)/2
+          do i=1,nx
+            f2(i,1)=-f2(mod(i-1+is,nx)+1,2)
+          enddo
+        endif
+        if (inpol.eq.1) then
+c         do north-pole
+          is=(nx-1)/2
+          do i=1,nx
+            f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny)
+          enddo
+        endif
+      endif
+c     compute flux-convergence -> filter
+      if (misdat.eq.0.) then
+        do j=1,ny
+        do i=1,nx
+            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
+        enddo
+        enddo
+      else
+        do j=1,ny
+        do i=1,nx
+          if (a(i,j).eq.misdat) then
+            af(i,j)=misdat
+          else
+            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
+          endif
+        enddo
+        enddo
+      endif
+      end
+      subroutine pipo(var3d,p3d,lev,var,nx,ny,nz,mdv,mode)
+C     ====================================================
+C     Interpolates the 3d variable var3d on the pressure surface
+C     defined by lev. The interpolated field is returned as var.
+C     p3d denotes the 3d pressure array.
+C     mode determines the way of vertical interpolation:
+C       mode=0 is for linear interpolation
+C       mode=1 is for logarithmic interpolation
+      integer   nx,ny,nz,mode
+      real      lev,mdv
+      real      var3d(nx,ny,nz),p3d(nx,ny,nz),var(nx,ny)
+      integer   i,j,k
+      real      kind
+      real      int3dm
+      do i=1,nx
+      do j=1,ny
+        kind=0.
+        do k=1,nz-1
+          if ((p3d(i,j,k).ge.lev).and.(p3d(i,j,k+1).le.lev)) then
+            kind=float(k)+(p3d(i,j,k)-lev)/
+     >                   (p3d(i,j,k)-p3d(i,j,k+1))
+            goto 100
+          endif
+        enddo
+continue
+        if (kind.eq.0.) then
+          var(i,j)=mdv
+        else
+          var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)
+        endif
+      enddo
+      enddo
+      return
+      end
+      subroutine thipo(var3d,th3d,lev,var,nx,ny,nz,mdv,mode)
+C     ======================================================
+C     Interpolates the 3d variable var3d on the isentropic surface
+C     defined by lev. The interpolated field is returned as var.
+C     th3d denotes the 3d theta array.
+C     mode determines the way of vertical interpolation:
+C       mode=0 is for linear interpolation
+C       mode=1 is for logarithmic interpolation
+      integer   nx,ny,nz,mode
+      real      lev,mdv
+      real      var3d(nx,ny,nz),th3d(nx,ny,nz),var(nx,ny)
+      integer   i,j,k
+      real      kind
+      real      int3dm
+      do i=1,nx
+      do j=1,ny
+        kind=0
+        do k=1,nz-1
+          if ((th3d(i,j,k).le.lev).and.(th3d(i,j,k+1).ge.lev)) then
+            kind=float(k)+(th3d(i,j,k)-lev)/
+     >                   (th3d(i,j,k)-th3d(i,j,k+1))
+            goto 100
+          endif
+        enddo
+continue
+        if (kind.eq.0) then
+          var(i,j)=mdv
+        else
+          var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)
+        endif
+      enddo
+      enddo
+      return
+      end

Subversion Repositories pvinversion.ecmwf

Compare Revisions

Ignore whitespace Rev 2 → Rev 3