WebSVN – lagranto.ecmwf – Diff – /trunk/lib/ioinp_cdo.f/


c     ************************************************************
c     * This package provides input routines to read the wind    *
c     * and other fields from IVE necdf files. The routines are  *
c     *                                                          *
c     * 1) input_open  : to open a data file                     *
c     * 2) input_grid  : to read the grid information, including *
c     *                  the vertical levels                     *
c     * 3) input_wind  : to read the wind components             *
c     * 4) input_close : to close an input file                  *
c     *                                                          *
c     * The file is characterised by an filename <filename> and  *
c     * a file identifier <fid>. The horizontal grid is given by *
c     * <xmin,xmax,ymin,ymax,dx,dy,nx,ny> where the pole of the  *
c     * rotated grid is given by <pollon,pollat>. The vertical   *
c     * grid is characterised by the surface pressure <ps> and   *
c     * the pressure at staggered <slev> and unstaggered <ulev>  *
c     * levels. The number of levels is given by <nz>. Finally,  *
c     * the retrieval of the wind <field> with name <fieldname>  *
c     * is characterised by a <time> and a missing data value    *
c     * <mdv>.                                                   *
c     *                                                          *
c     * Author: Michael Sprenger, Autumn 2008                    *
c     ************************************************************
 
c     ------------------------------------------------------------
c     Open input file
c     ------------------------------------------------------------
 
      subroutine input_open (fid,filename)
 
c     Open the input file with filename <filename> and return the
c     file identifier <fid> for further reference. 
 
      use netcdf
 
      implicit none
 
c     Declaration of subroutine parameters
      integer      fid              ! File identifier
      character*80 filename         ! Filename
 
c     Declaration of auxiliary variables
      integer      ierr
 
c     Open netcdf file
      ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
 
c     Exception handling
      return
 
      end
      
c     ------------------------------------------------------------
c     Read information about the grid
c     ------------------------------------------------------------
      
      subroutine input_grid 
     >                   (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >                    time,pollon,pollat,p3,ps,nz,ak,bk,stagz,
     >                    timecheck)
 
c     Read grid information at <time> from file with identifier <fid>. 
c     The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>
c     with pole position at <pollon,pollat> and grid dimension <nx,ny>.
c     The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either
c     on the staggered or unstaggered grid (with <stagz> as the flag).
c     The surface pressure is given in <ps(nx,ny)>. If <fid> is negative, 
c     only the grid dimensions and grid parameters (xmin...pollat,nz) are 
c     determined and returned (this is needed for dynamical allocation of 
c     memory).
 
      use netcdf
 
      implicit none
 
c     Declaration of subroutine parameters 
      integer      fid                 ! File identifier
      real         xmin,xmax,ymin,ymax ! Domain size
      real         dx,dy               ! Horizontal resolution
      integer      nx,ny,nz            ! Grid dimensions
      real         pollon,pollat       ! Longitude and latitude of pole
      real         p3(nx,ny,nz)        ! Staggered levels
      real         ps(nx,ny)           ! Surface pressure
      real         time                ! Time of the grid information
      real         ak(nz),bk(nz)       ! Ak and Bk for layers or levels
      real         stagz               ! Vertical staggering (0 or -0.5)
      character*80 fieldname           ! Variable from which to take grid info
      character*80 timecheck           ! Either 'yes' or 'no'
      
c     Numerical and physical parameters
      real          eps                 ! Numerical epsilon
      parameter    (eps=0.001)
 
c     Netcdf variables
      integer      vardim(4)
      real         varmin(4),varmax(4)
      real         mdv
      real         stag(4)
      integer      ndim
      character*80 cstfile
      integer      cstid
 
      integer      nvars
      character*80 vars(100)
      integer        dimids (nf90_max_var_dims)
      character*80   dimname(nf90_max_var_dims)
      real,allocatable, dimension (:)     :: lon,lat,lev
      real,allocatable, dimension (:)     :: times
      real,allocatable, dimension (:,:)   :: tmp2
      real,allocatable, dimension (:,:,:) :: tmp3
      real,allocatable, dimension (:)     :: aktmp,bktmp
      character*80  units
 
c     Auxiliary variables
      integer      ierr       
      integer      i,j,k
      integer      isok
      real         tmp(200)
      character*80 varname
      real         rtime
      integer      varid
      integer      cdfid
      integer      stat
      real         delta
      integer      closear
      real         maxps,minps
 
c     ------ Set file identifier --------------------------------------
      if (fid.lt.0) then
        cdfid = -fid
      else 
        cdfid = fid
      endif
 
c     Special handling if 3D pressure is
      if ( fieldname.eq.'P' ) then
         fieldname = 'U'
      endif
 
c     Get number of dimensions of variable -> ndim
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_variable(cdfid, varid, ndims  = ndim)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      if ( ndim.ne.4 ) then
          print*,' ERROR: netCDF variables need to be 4D'
          print*,'      ',trim(fieldname)
          stop
      endif
 
c     Get dimensions -> vardim(1:ndim),dimname(1:ndim)
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_variable(cdfid, varid, 
     >                                   dimids = dimids(1:ndim))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      do i=1,ndim
           ierr = nf90_inquire_dimension(cdfid, dimids(i), 
     >                               name = dimname(i) )
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
           ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      enddo
 
c     Check whether the list of dimensions is OK
      if ( ( dimname(1).ne.'lon'  ).or.
     >     ( dimname(2).ne.'lat'  ).or. 
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
     >     ( dimname(4).ne.'time' ) )
     >then
        print*,' ERROR: the dimensions of the variable are not correct'
        print*,'        expected -> lon / lat / lev / time'
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
        stop
      endif
 
c     Allocate memory for reading arrays
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
      allocate(lon(vardim(1)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lon     ***' 
      allocate(lat(vardim(2)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lat     ***' 
      allocate(times(vardim(4)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array times   ***'
      allocate(aktmp(vardim(3)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
      allocate(bktmp(vardim(3)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
 
c     Get domain longitudes and latitudes
      varname = dimname(1)
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,lon)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname = dimname(2)
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,lat)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
c     Get ak and bk
      varname='hyam'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,aktmp)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname='hybm'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,bktmp)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 
c     Check that unit of ak is in hPa - if necessary correct it
      varname='hyam'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_att(cdfid, varid, "units", units)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      if ( units.eq.'Pa' ) then
 
         do k=1,vardim(3)
            aktmp(k) = 0.01 * aktmp(k)
         enddo
      endif
 
c     Get time information (check if time is correct)
      varname = 'time'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,times)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      isok=0
      do i=1,vardim(4)
        if (abs(time-times(i)).lt.eps) then
               isok = 1
               rtime = times(i)
        elseif (timecheck.eq.'no') then
               isok = 1
               rtime = times(1)
        endif
      enddo
      if ( isok.eq.0 ) then
         print*,' ERROR: time ',rtime,' not found on netCDF file' 
         stop
      endif
 
c     Read surface pressure
      varname='PS'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,tmp2)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
    
c     Check that surface pressure is in hPa
      maxps = -1.e39
      minps =  1.e39
      do i=1,vardim(1)
        do j=1,vardim(2)
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
             if (tmp2(i,j).lt.minps) minps = tmp2(i,j)
        enddo
      enddo
      if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then
         print*,' ERROR: surface pressre PS must be in hPa'
         print*,'       ',maxps,minps
         stop
      endif
 
c     Calculate layer and level pressures
      do i=1,vardim(1)
         do j=1,vardim(2)
               do k=1,vardim(3)
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
               enddo
         enddo
      enddo
 
c     Set the grid dimensions and constants
      nx      = vardim(1)
      ny      = vardim(2)
      nz      = vardim(3)
      xmin    = lon(1)
      ymin    = lat(1)
      xmax    = lon(nx)
      ymax    = lat(ny)
      dx      = (xmax-xmin)/real(nx-1)
      dy      = (ymax-ymin)/real(ny-1)
      pollon  = 0.
      pollat  = 90.
      stagz   = 0.
      delta   = xmax-xmin-360.
      if (abs(delta+dx).lt.eps) then
          xmax    = xmax + dx
          nx      = nx + 1
          closear = 1
      else
          closear = 0
      endif
 
c     Save the output arrays (if fid>0) - close arrays on request
      if ( fid.gt.0 ) then
 
         do j=1,vardim(2)
           do i=1,vardim(1)
             ps(i,j) = tmp2(i,j)
           enddo
           if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)
         enddo
 
         do j=1,vardim(2)
           do k=1,vardim(3)
             do i=1,vardim(1)
               p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)
             enddo
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
           enddo
         enddo
 
         do k=1,vardim(3)
            ak(k) = aktmp(vardim(3)-k+1)
            bk(k) = bktmp(vardim(3)-k+1)
         enddo
 
      endif
 
      return
      
      end
 
c     ------------------------------------------------------------
c     Read wind information
c     ------------------------------------------------------------
 
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
     >                       timecheck)
 
c     Read the wind component <fieldname> from the file with identifier
c     <fid> and save it in the 3d array <field>. The vertical staggering 
c     information is provided in <stagz> and gives the reference to either
c     the layer or level field from <input_grid>. A consistency check is
c     performed to have an agreement with the grid specified by <xmin,xmax,
c     ymin,ymax,dx,dy,nx,ny,nz>.
 
      use netcdf
 
      implicit none
 
c     Declaration of variables and parameters
      integer      fid                 ! File identifier
      character*80 fieldname           ! Name of the wind field
      integer      nx,ny,nz            ! Dimension of fields
      real         field(nx,ny,nz)     ! 3d wind field
      real         stagz               ! Staggering in the z direction
      real         mdv                 ! Missing data flag
      real         xmin,xmax,ymin,ymax ! Domain size
      real         dx,dy               ! Horizontal resolution
      real         time                ! Time
      character*80 timecheck           ! Either 'yes' or 'no'
 
c     Numerical and physical parameters
      real        eps                 ! Numerical epsilon
      parameter  (eps=0.001)
      real        notimecheck         ! 'Flag' for no time check
      parameter  (notimecheck=7.26537)
 
c     Netcdf variables
      integer      ierr
      character*80 varname
      integer      vardim(4)
      real         varmin(4),varmax(4)
      real         stag(4)
      integer      ndim
      real         times(10)
      integer      ntimes
      character*80 cstfile
      integer      cstid
      real         aklay(200),bklay(200),aklev(200),bklev(200)
      real         ps(nx,ny)
      integer      dimids (nf90_max_var_dims)
      character*80 dimname(nf90_max_var_dims)
      integer      varid
      integer      cdfid
      real,allocatable, dimension (:)     :: lon,lat
      real,allocatable, dimension (:,:)   :: tmp2
      real,allocatable, dimension (:,:,:) :: tmp3
 
c     Auxiliary variables
      integer      isok
      integer      i,j,k
      integer      nz1
      real         rtime
      integer      closear
      integer      stat
      real         delta
 
c     Get the number of dimensions -> ndim
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 
c     Get the dimensions of the arrays -> varid(1:ndim)
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_variable(fid, varid, 
     >                                   dimids = dimids(1:ndim))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      do i=1,ndim
           ierr = nf90_inquire_dimension(fid, dimids(i), 
     >                               name = dimname(i) )
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      enddo
 
c     Check whether the list of dimensions is OK
      if ( ( dimname(1).ne.'lon'  ).or.
     >     ( dimname(2).ne.'lat'  ).or.
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
     >     ( dimname(4).ne.'time' ) )
     >then
        print*,' ERROR: the dimensions of the variable are not correct'
        print*,'        expected -> lon / lat / lev / time'
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
        stop
      endif
 
c     Allocate memory for reading arrays - depending on <closear>
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
      allocate(lon(vardim(1)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
      allocate(lat(vardim(2)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
 
c     Get domain boundaries
      varname = dimname(1)
      ierr = NF90_INQ_VARID(fid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(fid,varid,lon)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname = dimname(2)
      ierr = NF90_INQ_VARID(fid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(fid,varid,lat)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 
c     Read data 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(fid,varid,tmp3)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
  
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
      if ( vardim(3).eq.1 ) then
         do i=1,vardim(1)
            do j=1,vardim(2)
               do k=1,vardim(3)
                  tmp3(i,j,k) = tmp3(i,j,1)
               enddo
            enddo
         enddo
      endif
 
c     Save the ouput array - close on request
 
      delta = varmax(1)-varmin(1)-360.
      if (abs(delta+dx).lt.eps) then
          closear = 1
      else
          closear = 0
      endif
 
      do j=1,vardim(2)
        do k=1,vardim(3)
          do i=1,vardim(1)
             field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
          enddo
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
        enddo
      enddo
         
c     Exit point
      return
 
      end
 
c     ------------------------------------------------------------
c     Close input file
c     ------------------------------------------------------------
 
      subroutine input_close(fid)
 
c     Close the input file with file identifier <fid>.
 
      use netcdf
 
      implicit none
 
c     Declaration of subroutine parameters
      integer fid
 
c     Auxiliary variables
      integer ierr
 
c     Close file
      ierr = NF90_CLOSE(fid)
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
      end
      
c     ------------------------------------------------------------
c     Get a list of variables on netCDF file
c     ------------------------------------------------------------
 
      subroutine input_getvars(fid,vnam,nvars)
 
c     List of variables on netCDF file
 
      use netcdf
 
      implicit none
 
c     Declaration of subroutine parameters
      integer      fid
      integer      nvars
      character*80 vnam(200)
 
c     Auxiliary variables
      integer ierr
      integer i
      integer nDims,nGlobalAtts,unlimdimid
 
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
      do i=1,nVars
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
      enddo
 
      end

Rev 21	Rev 23
1	`c ************************************************************`	1	`c ************************************************************`
2	`c * This package provides input routines to read the wind *`	2	`c * This package provides input routines to read the wind *`
3	`c * and other fields from IVE necdf files. The routines are *`	3	`c * and other fields from IVE necdf files. The routines are *`
4	`c * *`	4	`c * *`
5	`c * 1) input_open : to open a data file *`	5	`c * 1) input_open : to open a data file *`
6	`c * 2) input_grid : to read the grid information, including *`	6	`c * 2) input_grid : to read the grid information, including *`
7	`c * the vertical levels *`	7	`c * the vertical levels *`
8	`c * 3) input_wind : to read the wind components *`	8	`c * 3) input_wind : to read the wind components *`
9	`c * 4) input_close : to close an input file *`	9	`c * 4) input_close : to close an input file *`
10	`c * *`	10	`c * *`
11	`c * The file is characterised by an filename <filename> and *`	11	`c * The file is characterised by an filename <filename> and *`
12	`c * a file identifier <fid>. The horizontal grid is given by *`	12	`c * a file identifier <fid>. The horizontal grid is given by *`
13	`c * <xmin,xmax,ymin,ymax,dx,dy,nx,ny> where the pole of the *`	13	`c * <xmin,xmax,ymin,ymax,dx,dy,nx,ny> where the pole of the *`
14	`c * rotated grid is given by <pollon,pollat>. The vertical *`	14	`c * rotated grid is given by <pollon,pollat>. The vertical *`
15	`c * grid is characterised by the surface pressure <ps> and *`	15	`c * grid is characterised by the surface pressure <ps> and *`
16	`c * the pressure at staggered <slev> and unstaggered <ulev> *`	16	`c * the pressure at staggered <slev> and unstaggered <ulev> *`
17	`c * levels. The number of levels is given by <nz>. Finally, *`	17	`c * levels. The number of levels is given by <nz>. Finally, *`
18	`c * the retrieval of the wind <field> with name <fieldname> *`	18	`c * the retrieval of the wind <field> with name <fieldname> *`
19	`c * is characterised by a <time> and a missing data value *`	19	`c * is characterised by a <time> and a missing data value *`
20	`c * <mdv>. *`	20	`c * <mdv>. *`
21	`c * *`	21	`c * *`
22	`c * Author: Michael Sprenger, Autumn 2008 *`	22	`c * Author: Michael Sprenger, Autumn 2008 *`
23	`c ************************************************************`	23	`c ************************************************************`
24		24
25	`c ------------------------------------------------------------`	25	`c ------------------------------------------------------------`
26	`c Open input file`	26	`c Open input file`
27	`c ------------------------------------------------------------`	27	`c ------------------------------------------------------------`
28		28
29	`subroutine input_open (fid,filename)`	29	`subroutine input_open (fid,filename)`
30		30
31	`c Open the input file with filename <filename> and return the`	31	`c Open the input file with filename <filename> and return the`
32	`c file identifier <fid> for further reference.`	32	`c file identifier <fid> for further reference.`
33		33
34	`use netcdf`	34	`use netcdf`
35		35
36	`implicit none`	36	`implicit none`
37		37
38	`c Declaration of subroutine parameters`	38	`c Declaration of subroutine parameters`
39	`integer fid ! File identifier`	39	`integer fid ! File identifier`
40	`character*80 filename ! Filename`	40	`character*80 filename ! Filename`
41		41
42	`c Declaration of auxiliary variables`	42	`c Declaration of auxiliary variables`
43	`integer ierr`	43	`integer ierr`
44		44
45	`c Open netcdf file`	45	`c Open netcdf file`
46	`ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)`	46	`ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)`
47	`IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)`	47	`IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)`
48		48
49	`c Exception handling`	49	`c Exception handling`
50	`return`	50	`return`
51		51
52	`end`	52	`end`
53		53
54	`c ------------------------------------------------------------`	54	`c ------------------------------------------------------------`
55	`c Read information about the grid`	55	`c Read information about the grid`
56	`c ------------------------------------------------------------`	56	`c ------------------------------------------------------------`
57		57
58	`subroutine input_grid`	58	`subroutine input_grid`
59	`> (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`	59	`> (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`
60	`> time,pollon,pollat,p3,ps,nz,ak,bk,stagz,`	60	`> time,pollon,pollat,p3,ps,nz,ak,bk,stagz,`
61	`> timecheck)`	61	`> timecheck)`
62		62
63	`c Read grid information at <time> from file with identifier <fid>.`	63	`c Read grid information at <time> from file with identifier <fid>.`
64	`c The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>`	64	`c The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>`
65	`c with pole position at <pollon,pollat> and grid dimension <nx,ny>.`	65	`c with pole position at <pollon,pollat> and grid dimension <nx,ny>.`
66	`c The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either`	66	`c The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either`
67	`c on the staggered or unstaggered grid (with <stagz> as the flag).`	67	`c on the staggered or unstaggered grid (with <stagz> as the flag).`
68	`c The surface pressure is given in <ps(nx,ny)>. If <fid> is negative,`	68	`c The surface pressure is given in <ps(nx,ny)>. If <fid> is negative,`
69	`c only the grid dimensions and grid parameters (xmin...pollat,nz) are`	69	`c only the grid dimensions and grid parameters (xmin...pollat,nz) are`
70	`c determined and returned (this is needed for dynamical allocation of`	70	`c determined and returned (this is needed for dynamical allocation of`
71	`c memory).`	71	`c memory).`
72		72
73	`use netcdf`	73	`use netcdf`
74		74
75	`implicit none`	75	`implicit none`
76		76
77	`c Declaration of subroutine parameters`	77	`c Declaration of subroutine parameters`
78	`integer fid ! File identifier`	78	`integer fid ! File identifier`
79	`real xmin,xmax,ymin,ymax ! Domain size`	79	`real xmin,xmax,ymin,ymax ! Domain size`
80	`real dx,dy ! Horizontal resolution`	80	`real dx,dy ! Horizontal resolution`
81	`integer nx,ny,nz ! Grid dimensions`	81	`integer nx,ny,nz ! Grid dimensions`
82	`real pollon,pollat ! Longitude and latitude of pole`	82	`real pollon,pollat ! Longitude and latitude of pole`
83	`real p3(nx,ny,nz) ! Staggered levels`	83	`real p3(nx,ny,nz) ! Staggered levels`
84	`real ps(nx,ny) ! Surface pressure`	84	`real ps(nx,ny) ! Surface pressure`
85	`real time ! Time of the grid information`	85	`real time ! Time of the grid information`
86	`real ak(nz),bk(nz) ! Ak and Bk for layers or levels`	86	`real ak(nz),bk(nz) ! Ak and Bk for layers or levels`
87	`real stagz ! Vertical staggering (0 or -0.5)`	87	`real stagz ! Vertical staggering (0 or -0.5)`
88	`character*80 fieldname ! Variable from which to take grid info`	88	`character*80 fieldname ! Variable from which to take grid info`
89	`character*80 timecheck ! Either 'yes' or 'no'`	89	`character*80 timecheck ! Either 'yes' or 'no'`
90		90
91	`c Numerical and physical parameters`	91	`c Numerical and physical parameters`
92	`real eps ! Numerical epsilon`	92	`real eps ! Numerical epsilon`
93	`parameter (eps=0.001)`	93	`parameter (eps=0.001)`
94		94
95	`c Netcdf variables`	95	`c Netcdf variables`
96	`integer vardim(4)`	96	`integer vardim(4)`
97	`real varmin(4),varmax(4)`	97	`real varmin(4),varmax(4)`
98	`real mdv`	98	`real mdv`
99	`real stag(4)`	99	`real stag(4)`
100	`integer ndim`	100	`integer ndim`
101	`character*80 cstfile`	101	`character*80 cstfile`
102	`integer cstid`	102	`integer cstid`
103	`real aklay(nz),bklay(nz)`	-
104	`integer nvars`	103	`integer nvars`
105	`character*80 vars(100)`	104	`character*80 vars(100)`
106	`integer dimids (nf90_max_var_dims)`	105	`integer dimids (nf90_max_var_dims)`
107	`character*80 dimname(nf90_max_var_dims)`	106	`character*80 dimname(nf90_max_var_dims)`
108	`real,allocatable, dimension (:) :: lon,lat,lev`	107	`real,allocatable, dimension (:) :: lon,lat,lev`
109	`real,allocatable, dimension (:) :: times`	108	`real,allocatable, dimension (:) :: times`
110	`real,allocatable, dimension (:,:) :: tmp2`	109	`real,allocatable, dimension (:,:) :: tmp2`
111	`real,allocatable, dimension (:,:,:) :: tmp3`	110	`real,allocatable, dimension (:,:,:) :: tmp3`
-		111	`real,allocatable, dimension (:) :: aktmp,bktmp`
112	`character*80 units`	112	`character*80 units`
113		113
114	`c Auxiliary variables`	114	`c Auxiliary variables`
115	`integer ierr`	115	`integer ierr`
116	`integer i,j,k`	116	`integer i,j,k`
117	`integer isok`	117	`integer isok`
118	`real tmp(200)`	118	`real tmp(200)`
119	`character*80 varname`	119	`character*80 varname`
120	`real rtime`	120	`real rtime`
121	`integer varid`	121	`integer varid`
122	`integer cdfid`	122	`integer cdfid`
123	`integer stat`	123	`integer stat`
124	`real delta`	124	`real delta`
125	`integer closear`	125	`integer closear`
126	`real maxps,minps`	126	`real maxps,minps`
127		127
128	`c ------ Set file identifier --------------------------------------`	128	`c ------ Set file identifier --------------------------------------`
129	`if (fid.lt.0) then`	129	`if (fid.lt.0) then`
130	`cdfid = -fid`	130	`cdfid = -fid`
131	`else`	131	`else`
132	`cdfid = fid`	132	`cdfid = fid`
133	`endif`	133	`endif`
134		134
135	`c Special handling if 3D pressure is`	135	`c Special handling if 3D pressure is`
136	`if ( fieldname.eq.'P' ) then`	136	`if ( fieldname.eq.'P' ) then`
137	`fieldname = 'U'`	137	`fieldname = 'U'`
138	`endif`	138	`endif`
139		139
140	`c Get number of dimensions of variable -> ndim`	140	`c Get number of dimensions of variable -> ndim`
141	`ierr = NF90_INQ_VARID(cdfid,fieldname,varid)`	141	`ierr = NF90_INQ_VARID(cdfid,fieldname,varid)`
142	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	142	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
143	`ierr = nf90_inquire_variable(cdfid, varid, ndims = ndim)`	143	`ierr = nf90_inquire_variable(cdfid, varid, ndims = ndim)`
144	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	144	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
145	`if ( ndim.ne.4 ) then`	145	`if ( ndim.ne.4 ) then`
146	`print*,' ERROR: netCDF variables need to be 4D'`	146	`print*,' ERROR: netCDF variables need to be 4D'`
147	`print*,' ',trim(fieldname)`	147	`print*,' ',trim(fieldname)`
148	`stop`	148	`stop`
149	`endif`	149	`endif`
150		150
151	`c Get dimensions -> vardim(1:ndim),dimname(1:ndim)`	151	`c Get dimensions -> vardim(1:ndim),dimname(1:ndim)`
152	`ierr = NF90_INQ_VARID(cdfid,fieldname,varid)`	152	`ierr = NF90_INQ_VARID(cdfid,fieldname,varid)`
153	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	153	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
154	`ierr = nf90_inquire_variable(cdfid, varid,`	154	`ierr = nf90_inquire_variable(cdfid, varid,`
155	`> dimids = dimids(1:ndim))`	155	`> dimids = dimids(1:ndim))`
156	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	156	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
157	`do i=1,ndim`	157	`do i=1,ndim`
158	`ierr = nf90_inquire_dimension(cdfid, dimids(i),`	158	`ierr = nf90_inquire_dimension(cdfid, dimids(i),`
159	`> name = dimname(i) )`	159	`> name = dimname(i) )`
160	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	160	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
161	`ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))`	161	`ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))`
162	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	162	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
163	`enddo`	163	`enddo`
164		164
165	`c Check whether the list of dimensions is OK`	165	`c Check whether the list of dimensions is OK`
166	`if ( ( dimname(1).ne.'lon' ).or.`	166	`if ( ( dimname(1).ne.'lon' ).or.`
167	`> ( dimname(2).ne.'lat' ).or.`	167	`> ( dimname(2).ne.'lat' ).or.`
168	`> ( dimname(3).ne.'lev' ).and.( dimname(3).ne.'lev_2' ).or.`	168	`> ( dimname(3).ne.'lev' ).and.( dimname(3).ne.'lev_2' ).or.`
169	`> ( dimname(4).ne.'time' ) )`	169	`> ( dimname(4).ne.'time' ) )`
170	`>then`	170	`>then`
171	`print*,' ERROR: the dimensions of the variable are not correct'`	171	`print*,' ERROR: the dimensions of the variable are not correct'`
172	`print*,' expected -> lon / lat / lev / time'`	172	`print*,' expected -> lon / lat / lev / time'`
173	`print*, ( trim(dimname(i))//' / ',i=1,ndim )`	173	`print*, ( trim(dimname(i))//' / ',i=1,ndim )`
174	`stop`	174	`stop`
175	`endif`	175	`endif`
176		176
177	`c Allocate memory for reading arrays`	177	`c Allocate memory for reading arrays`
178	`allocate(tmp2(vardim(1),vardim(2)),stat=stat)`	178	`allocate(tmp2(vardim(1),vardim(2)),stat=stat)`
179	`if (stat.ne.0) print,' error allocating array tmp2 *'`	179	`if (stat.ne.0) print,' error allocating array tmp2 *'`
180	`allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)`	180	`allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)`
181	`if (stat.ne.0) print,' error allocating array tmp3 *'`	181	`if (stat.ne.0) print,' error allocating array tmp3 *'`
182	`allocate(lon(vardim(1)),stat=stat)`	182	`allocate(lon(vardim(1)),stat=stat)`
183	`if (stat.ne.0) print,' error allocating array lon *'`	183	`if (stat.ne.0) print,' error allocating array lon *'`
184	`allocate(lat(vardim(2)),stat=stat)`	184	`allocate(lat(vardim(2)),stat=stat)`
185	`if (stat.ne.0) print,' error allocating array lat *'`	185	`if (stat.ne.0) print,' error allocating array lat *'`
186	`allocate(times(vardim(4)),stat=stat)`	186	`allocate(times(vardim(4)),stat=stat)`
187	`if (stat.ne.0) print,' error allocating array times *'`	187	`if (stat.ne.0) print,' error allocating array times *'`
-		188	`allocate(aktmp(vardim(3)),stat=stat)`
-		189	`if (stat.ne.0) print,' error allocating array aktmp *'`
-		190	`allocate(bktmp(vardim(3)),stat=stat)`
-		191	`if (stat.ne.0) print,' error allocating array bktmp *'`
188		192
189	`c Get domain longitudes and latitudes`	193	`c Get domain longitudes and latitudes`
190	`varname = dimname(1)`	194	`varname = dimname(1)`
191	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	195	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
192	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	196	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
193	`ierr = nf90_get_var(cdfid,varid,lon)`	197	`ierr = nf90_get_var(cdfid,varid,lon)`
194	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	198	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
195	`varname = dimname(2)`	199	`varname = dimname(2)`
196	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	200	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
197	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	201	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
198	`ierr = nf90_get_var(cdfid,varid,lat)`	202	`ierr = nf90_get_var(cdfid,varid,lat)`
199	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	203	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
200		204
201	`c Get ak and bk`	205	`c Get ak and bk`
202	`varname='hyam'`	206	`varname='hyam'`
203	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	207	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
204	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	208	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
205	`ierr = nf90_get_var(cdfid,varid,ak)`	209	`ierr = nf90_get_var(cdfid,varid,aktmp)`
206	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	210	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
207	`varname='hybm'`	211	`varname='hybm'`
208	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	212	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
209	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	213	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
210	`ierr = nf90_get_var(cdfid,varid,bk)`	214	`ierr = nf90_get_var(cdfid,varid,bktmp)`
211	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	215	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
212		216
213	`c Check that unit of ak is in hPa - if necessary correct it`	217	`c Check that unit of ak is in hPa - if necessary correct it`
214	`varname='hyam'`	218	`varname='hyam'`
215	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	219	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
216	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	220	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
217	`ierr = nf90_get_att(cdfid, varid, "units", units)`	221	`ierr = nf90_get_att(cdfid, varid, "units", units)`
218	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	222	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
219	`if ( units.eq.'Pa' ) then`	223	`if ( units.eq.'Pa' ) then`
220	`print*,'Hallo'`	-
221	`do i=1,vardim(3)`	224	`do k=1,vardim(3)`
222	`ak(k) = 0.01 * ak(k)`	225	`aktmp(k) = 0.01 * aktmp(k)`
223	`enddo`	226	`enddo`
224	`endif`	227	`endif`
225		228
226	`c Get time information (check if time is correct)`	229	`c Get time information (check if time is correct)`
227	`varname = 'time'`	230	`varname = 'time'`
228	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	231	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
229	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	232	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
230	`ierr = nf90_get_var(cdfid,varid,times)`	233	`ierr = nf90_get_var(cdfid,varid,times)`
231	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	234	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
232	`isok=0`	235	`isok=0`
233	`do i=1,vardim(4)`	236	`do i=1,vardim(4)`
234	`if (abs(time-times(i)).lt.eps) then`	237	`if (abs(time-times(i)).lt.eps) then`
235	`isok = 1`	238	`isok = 1`
236	`rtime = times(i)`	239	`rtime = times(i)`
237	`elseif (timecheck.eq.'no') then`	240	`elseif (timecheck.eq.'no') then`
238	`isok = 1`	241	`isok = 1`
239	`rtime = times(1)`	242	`rtime = times(1)`
240	`endif`	243	`endif`
241	`enddo`	244	`enddo`
242	`if ( isok.eq.0 ) then`	245	`if ( isok.eq.0 ) then`
243	`print*,' ERROR: time ',rtime,' not found on netCDF file'`	246	`print*,' ERROR: time ',rtime,' not found on netCDF file'`
244	`stop`	247	`stop`
245	`endif`	248	`endif`
246		249
247	`c Read surface pressure`	250	`c Read surface pressure`
248	`varname='PS'`	251	`varname='PS'`
249	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	252	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
250	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	253	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
251	`ierr = nf90_get_var(cdfid,varid,tmp2)`	254	`ierr = nf90_get_var(cdfid,varid,tmp2)`
252	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	255	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
253		256
254	`c Check that surface pressure is in hPa`	257	`c Check that surface pressure is in hPa`
255	`maxps = -1.e39`	258	`maxps = -1.e39`
256	`minps = 1.e39`	259	`minps = 1.e39`
257	`do i=1,vardim(1)`	260	`do i=1,vardim(1)`
258	`do j=1,vardim(2)`	261	`do j=1,vardim(2)`
259	`if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)`	262	`if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)`
260	`if (tmp2(i,j).lt.minps) minps = tmp2(i,j)`	263	`if (tmp2(i,j).lt.minps) minps = tmp2(i,j)`
261	`enddo`	264	`enddo`
262	`enddo`	265	`enddo`
263	`if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then`	266	`if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then`
264	`print*,' ERROR: surface pressre PS must be in hPa'`	267	`print*,' ERROR: surface pressre PS must be in hPa'`
265	`print*,' ',maxps,minps`	268	`print*,' ',maxps,minps`
266	`stop`	269	`stop`
267	`endif`	270	`endif`
268		271
269	`c Calculate layer and level pressures`	272	`c Calculate layer and level pressures`
270	`do i=1,vardim(1)`	273	`do i=1,vardim(1)`
271	`do j=1,vardim(2)`	274	`do j=1,vardim(2)`
272	`do k=1,vardim(3)`	275	`do k=1,vardim(3)`
273	`tmp3(i,j,k)=ak(k)+bk(k)*tmp2(i,j)`	276	`tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)`
274	`enddo`	277	`enddo`
275	`enddo`	278	`enddo`
276	`enddo`	279	`enddo`
277		280
278	`c Set the grid dimensions and constants`	281	`c Set the grid dimensions and constants`
279	`nx = vardim(1)`	282	`nx = vardim(1)`
280	`ny = vardim(2)`	283	`ny = vardim(2)`
281	`nz = vardim(3)`	284	`nz = vardim(3)`
282	`xmin = lon(1)`	285	`xmin = lon(1)`
283	`ymin = lat(1)`	286	`ymin = lat(1)`
284	`xmax = lon(nx)`	287	`xmax = lon(nx)`
285	`ymax = lat(ny)`	288	`ymax = lat(ny)`
286	`dx = (xmax-xmin)/real(nx-1)`	289	`dx = (xmax-xmin)/real(nx-1)`
287	`dy = (ymax-ymin)/real(ny-1)`	290	`dy = (ymax-ymin)/real(ny-1)`
288	`pollon = 0.`	291	`pollon = 0.`
289	`pollat = 90.`	292	`pollat = 90.`
290	`stagz = 0.`	293	`stagz = 0.`
291	`delta = xmax-xmin-360.`	294	`delta = xmax-xmin-360.`
292	`if (abs(delta+dx).lt.eps) then`	295	`if (abs(delta+dx).lt.eps) then`
293	`xmax = xmax + dx`	296	`xmax = xmax + dx`
294	`nx = nx + 1`	297	`nx = nx + 1`
295	`closear = 1`	298	`closear = 1`
296	`else`	299	`else`
297	`closear = 0`	300	`closear = 0`
298	`endif`	301	`endif`
299		302
300	`c Save the output arrays (if fid>0) - close arrays on request`	303	`c Save the output arrays (if fid>0) - close arrays on request`
301	`if ( fid.gt.0 ) then`	304	`if ( fid.gt.0 ) then`
302		305
303	`do j=1,vardim(2)`	306	`do j=1,vardim(2)`
304	`do i=1,vardim(1)`	307	`do i=1,vardim(1)`
305	`ps(i,j) = tmp2(i,j)`	308	`ps(i,j) = tmp2(i,j)`
306	`enddo`	309	`enddo`
307	`if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)`	310	`if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)`
308	`enddo`	311	`enddo`
309		312
310	`do j=1,vardim(2)`	313	`do j=1,vardim(2)`
311	`do k=1,vardim(3)`	314	`do k=1,vardim(3)`
312	`do i=1,vardim(1)`	315	`do i=1,vardim(1)`
313	`p3(i,j,k) = tmp3(i,j,k)`	316	`p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)`
314	`enddo`	317	`enddo`
315	`if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)`	318	`if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)`
316	`enddo`	319	`enddo`
317	`enddo`	320	`enddo`
-		321
-		322	`do k=1,vardim(3)`
-		323	`ak(k) = aktmp(vardim(3)-k+1)`
-		324	`bk(k) = bktmp(vardim(3)-k+1)`
-		325	`enddo`
-		326
318	`endif`	327	`endif`
319		328
320	`return`	329	`return`
321		330
322	`end`	331	`end`
323		332
324	`c ------------------------------------------------------------`	333	`c ------------------------------------------------------------`
325	`c Read wind information`	334	`c Read wind information`
326	`c ------------------------------------------------------------`	335	`c ------------------------------------------------------------`
327		336
328	`subroutine input_wind (fid,fieldname,field,time,stagz,mdv,`	337	`subroutine input_wind (fid,fieldname,field,time,stagz,mdv,`
329	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,`	338	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,`
330	`> timecheck)`	339	`> timecheck)`
331		340
332	`c Read the wind component <fieldname> from the file with identifier`	341	`c Read the wind component <fieldname> from the file with identifier`
333	`c <fid> and save it in the 3d array <field>. The vertical staggering`	342	`c <fid> and save it in the 3d array <field>. The vertical staggering`
334	`c information is provided in <stagz> and gives the reference to either`	343	`c information is provided in <stagz> and gives the reference to either`
335	`c the layer or level field from <input_grid>. A consistency check is`	344	`c the layer or level field from <input_grid>. A consistency check is`
336	`c performed to have an agreement with the grid specified by <xmin,xmax,`	345	`c performed to have an agreement with the grid specified by <xmin,xmax,`
337	`c ymin,ymax,dx,dy,nx,ny,nz>.`	346	`c ymin,ymax,dx,dy,nx,ny,nz>.`
338		347
339	`use netcdf`	348	`use netcdf`
340		349
341	`implicit none`	350	`implicit none`
342		351
343	`c Declaration of variables and parameters`	352	`c Declaration of variables and parameters`
344	`integer fid ! File identifier`	353	`integer fid ! File identifier`
345	`character*80 fieldname ! Name of the wind field`	354	`character*80 fieldname ! Name of the wind field`
346	`integer nx,ny,nz ! Dimension of fields`	355	`integer nx,ny,nz ! Dimension of fields`
347	`real field(nx,ny,nz) ! 3d wind field`	356	`real field(nx,ny,nz) ! 3d wind field`
348	`real stagz ! Staggering in the z direction`	357	`real stagz ! Staggering in the z direction`
349	`real mdv ! Missing data flag`	358	`real mdv ! Missing data flag`
350	`real xmin,xmax,ymin,ymax ! Domain size`	359	`real xmin,xmax,ymin,ymax ! Domain size`
351	`real dx,dy ! Horizontal resolution`	360	`real dx,dy ! Horizontal resolution`
352	`real time ! Time`	361	`real time ! Time`
353	`character*80 timecheck ! Either 'yes' or 'no'`	362	`character*80 timecheck ! Either 'yes' or 'no'`
354		363
355	`c Numerical and physical parameters`	364	`c Numerical and physical parameters`
356	`real eps ! Numerical epsilon`	365	`real eps ! Numerical epsilon`
357	`parameter (eps=0.001)`	366	`parameter (eps=0.001)`
358	`real notimecheck ! 'Flag' for no time check`	367	`real notimecheck ! 'Flag' for no time check`
359	`parameter (notimecheck=7.26537)`	368	`parameter (notimecheck=7.26537)`
360		369
361	`c Netcdf variables`	370	`c Netcdf variables`
362	`integer ierr`	371	`integer ierr`
363	`character*80 varname`	372	`character*80 varname`
364	`integer vardim(4)`	373	`integer vardim(4)`
365	`real varmin(4),varmax(4)`	374	`real varmin(4),varmax(4)`
366	`real stag(4)`	375	`real stag(4)`
367	`integer ndim`	376	`integer ndim`
368	`real times(10)`	377	`real times(10)`
369	`integer ntimes`	378	`integer ntimes`
370	`character*80 cstfile`	379	`character*80 cstfile`
371	`integer cstid`	380	`integer cstid`
372	`real aklay(200),bklay(200),aklev(200),bklev(200)`	381	`real aklay(200),bklay(200),aklev(200),bklev(200)`
373	`real ps(nx,ny)`	382	`real ps(nx,ny)`
374	`integer dimids (nf90_max_var_dims)`	383	`integer dimids (nf90_max_var_dims)`
375	`character*80 dimname(nf90_max_var_dims)`	384	`character*80 dimname(nf90_max_var_dims)`
376	`integer varid`	385	`integer varid`
377	`integer cdfid`	386	`integer cdfid`
378	`real,allocatable, dimension (:) :: lon,lat`	387	`real,allocatable, dimension (:) :: lon,lat`
379	`real,allocatable, dimension (:,:) :: tmp2`	388	`real,allocatable, dimension (:,:) :: tmp2`
380	`real,allocatable, dimension (:,:,:) :: tmp3`	389	`real,allocatable, dimension (:,:,:) :: tmp3`
381		390
382	`c Auxiliary variables`	391	`c Auxiliary variables`
383	`integer isok`	392	`integer isok`
384	`integer i,j,k`	393	`integer i,j,k`
385	`integer nz1`	394	`integer nz1`
386	`real rtime`	395	`real rtime`
387	`integer closear`	396	`integer closear`
388	`integer stat`	397	`integer stat`
389	`real delta`	398	`real delta`
390		399
391	`c Get the number of dimensions -> ndim`	400	`c Get the number of dimensions -> ndim`
392	`ierr = NF90_INQ_VARID(fid,fieldname,varid)`	401	`ierr = NF90_INQ_VARID(fid,fieldname,varid)`
393	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	402	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
394	`ierr = nf90_inquire_variable(fid, varid, ndims = ndim)`	403	`ierr = nf90_inquire_variable(fid, varid, ndims = ndim)`
395	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	404	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
396		405
397	`c Get the dimensions of the arrays -> varid(1:ndim)`	406	`c Get the dimensions of the arrays -> varid(1:ndim)`
398	`ierr = NF90_INQ_VARID(fid,fieldname,varid)`	407	`ierr = NF90_INQ_VARID(fid,fieldname,varid)`
399	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	408	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
400	`ierr = nf90_inquire_variable(fid, varid,`	409	`ierr = nf90_inquire_variable(fid, varid,`
401	`> dimids = dimids(1:ndim))`	410	`> dimids = dimids(1:ndim))`
402	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	411	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
403	`do i=1,ndim`	412	`do i=1,ndim`
404	`ierr = nf90_inquire_dimension(fid, dimids(i),`	413	`ierr = nf90_inquire_dimension(fid, dimids(i),`
405	`> name = dimname(i) )`	414	`> name = dimname(i) )`
406	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	415	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
407	`ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))`	416	`ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))`
408	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	417	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
409	`enddo`	418	`enddo`
410		419
411	`c Check whether the list of dimensions is OK`	420	`c Check whether the list of dimensions is OK`
412	`if ( ( dimname(1).ne.'lon' ).or.`	421	`if ( ( dimname(1).ne.'lon' ).or.`
413	`> ( dimname(2).ne.'lat' ).or.`	422	`> ( dimname(2).ne.'lat' ).or.`
414	`> ( dimname(3).ne.'lev' ).and.( dimname(3).ne.'lev_2' ).or.`	423	`> ( dimname(3).ne.'lev' ).and.( dimname(3).ne.'lev_2' ).or.`
415	`> ( dimname(4).ne.'time' ) )`	424	`> ( dimname(4).ne.'time' ) )`
416	`>then`	425	`>then`
417	`print*,' ERROR: the dimensions of the variable are not correct'`	426	`print*,' ERROR: the dimensions of the variable are not correct'`
418	`print*,' expected -> lon / lat / lev / time'`	427	`print*,' expected -> lon / lat / lev / time'`
419	`print*, ( trim(dimname(i))//' / ',i=1,ndim )`	428	`print*, ( trim(dimname(i))//' / ',i=1,ndim )`
420	`stop`	429	`stop`
421	`endif`	430	`endif`
422		431
423	`c Allocate memory for reading arrays - depending on <closear>`	432	`c Allocate memory for reading arrays - depending on <closear>`
424	`allocate(tmp2(vardim(1),vardim(2)),stat=stat)`	433	`allocate(tmp2(vardim(1),vardim(2)),stat=stat)`
425	`if (stat.ne.0) print,' error allocating array tmp2 *'`	434	`if (stat.ne.0) print,' error allocating array tmp2 *'`
426	`allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)`	435	`allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)`
427	`if (stat.ne.0) print,' error allocating array tmp3 *'`	436	`if (stat.ne.0) print,' error allocating array tmp3 *'`
428	`allocate(lon(vardim(1)),stat=stat)`	437	`allocate(lon(vardim(1)),stat=stat)`
429	`if (stat.ne.0) print,' error allocating array lon *'`	438	`if (stat.ne.0) print,' error allocating array lon *'`
430	`allocate(lat(vardim(2)),stat=stat)`	439	`allocate(lat(vardim(2)),stat=stat)`
431	`if (stat.ne.0) print,' error allocating array lat *'`	440	`if (stat.ne.0) print,' error allocating array lat *'`
432		441
433	`c Get domain boundaries`	442	`c Get domain boundaries`
434	`varname = dimname(1)`	443	`varname = dimname(1)`
435	`ierr = NF90_INQ_VARID(fid,varname,varid)`	444	`ierr = NF90_INQ_VARID(fid,varname,varid)`
436	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	445	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
437	`ierr = nf90_get_var(fid,varid,lon)`	446	`ierr = nf90_get_var(fid,varid,lon)`
438	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	447	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
439	`varname = dimname(2)`	448	`varname = dimname(2)`
440	`ierr = NF90_INQ_VARID(fid,varname,varid)`	449	`ierr = NF90_INQ_VARID(fid,varname,varid)`
441	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	450	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
442	`ierr = nf90_get_var(fid,varid,lat)`	451	`ierr = nf90_get_var(fid,varid,lat)`
443	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	452	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
444		453
445	`c Read data`	454	`c Read data`
446	`ierr = NF90_INQ_VARID(fid,fieldname,varid)`	455	`ierr = NF90_INQ_VARID(fid,fieldname,varid)`
447	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	456	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
448	`ierr = nf90_get_var(fid,varid,tmp3)`	457	`ierr = nf90_get_var(fid,varid,tmp3)`
449	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	458	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
450		459
451	`c If the field is 2D, expand it to 3D - simple handling of 2D tracing`	460	`c If the field is 2D, expand it to 3D - simple handling of 2D tracing`
452	`if ( vardim(3).eq.1 ) then`	461	`if ( vardim(3).eq.1 ) then`
453	`do i=1,vardim(1)`	462	`do i=1,vardim(1)`
454	`do j=1,vardim(2)`	463	`do j=1,vardim(2)`
455	`do k=1,vardim(3)`	464	`do k=1,vardim(3)`
456	`tmp3(i,j,k) = tmp3(i,j,1)`	465	`tmp3(i,j,k) = tmp3(i,j,1)`
457	`enddo`	466	`enddo`
458	`enddo`	467	`enddo`
459	`enddo`	468	`enddo`
460	`endif`	469	`endif`
461		470
462	`c Save the ouput array - close on request`	471	`c Save the ouput array - close on request`
463	`dx = (varmax(1)-varmin(1))/real(vardim(1)-1)`	-
464	`delta = varmax(1)-varmin(1)-360.`	472	`delta = varmax(1)-varmin(1)-360.`
465	`if (abs(delta+dx).lt.eps) then`	473	`if (abs(delta+dx).lt.eps) then`
466	`closear = 1`	474	`closear = 1`
467	`else`	475	`else`
468	`closear = 0`	476	`closear = 0`
469	`endif`	477	`endif`
470		478
471	`do j=1,vardim(2)`	479	`do j=1,vardim(2)`
472	`do k=1,vardim(3)`	480	`do k=1,vardim(3)`
473	`do i=1,vardim(1)`	481	`do i=1,vardim(1)`
474	`field(i,j,k) = tmp3(i,j,k)`	482	`field(i,j,k) = tmp3(i,j,vardim(3)-k+1)`
475	`enddo`	483	`enddo`
476	`if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)`	484	`if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)`
477	`enddo`	485	`enddo`
478	`enddo`	486	`enddo`
479		487
480	`c Exit point`	488	`c Exit point`
481	`return`	489	`return`
482		490
483	`end`	491	`end`
484		492
485	`c ------------------------------------------------------------`	493	`c ------------------------------------------------------------`
486	`c Close input file`	494	`c Close input file`
487	`c ------------------------------------------------------------`	495	`c ------------------------------------------------------------`
488		496
489	`subroutine input_close(fid)`	497	`subroutine input_close(fid)`

Subversion Repositories lagranto.ecmwf

(root)/trunk/lib/ioinp_cdo.f/ – Rev 21 → 23