WebSVN – lagranto.ecmwf – Diff – /trunk/lib/ioinp.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. 
 
      implicit none
 
c     Declaration of subroutine parameters
      integer      fid              ! File identifier
      character*80 filename         ! Filename
 
c     Declaration of auxiliary variables
      integer      ierr
 
c     Open IVE netcdf file
      call cdfopn(filename,fid,ierr)
      if (ierr.ne.0) goto 900
      
c     Exception handling
      return
      
 900  print*,'Cannot open input file ',trim(filename)
      stop
 
      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).
 
      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
      real         times(10)
      integer      ntimes
      real         aklay(nz),bklay(nz),aklev(nz),bklev(nz)
      integer      nvars
      character*80 vars(100)
 
c     Auxiliary varaibles
      integer      ierr       
      integer      i,j,k
      integer      isok
      real         tmp(200)
      character*80 varname
      real         rtime
      integer      is2d
      integer      plev
 
c     Init the flag for 2D variables - assume a 3D field
      is2d = 0
 
c     Init the flag for pressure levels (PS is not needed)
      plev = 0
 
c     Inquire dimensions and grid constants if <fid> is negative
      if (fid.lt.0) then
 
c        Get grid info for the selected variable
         if ( fieldname.eq.'PLEV' ) then
            varname = 'PS'
            stagz   = 0.
            call getdef(-fid,varname,ndim,mdv,vardim,
     >                  varmin,varmax,stag,ierr)
            if (ierr.ne.0) goto 900
            call getcfn(-fid,cstfile,ierr)
            if (ierr.ne.0) goto 903
            call cdfopn(cstfile,cstid,ierr)
            if (ierr.ne.0) goto 903
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
            if (ierr.ne.0) goto 903
            call clscdf(cstid,ierr)
            if (ierr.ne.0) goto 903
            
         elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then       
            varname = 'PS'
            stagz   = -0.5
            call getdef(-fid,varname,ndim,mdv,vardim,
     >                  varmin,varmax,stag,ierr)
            if (ierr.ne.0) goto 900
            call getcfn(-fid,cstfile,ierr)
            if (ierr.ne.0) goto 903
            call cdfopn(cstfile,cstid,ierr)
            if (ierr.ne.0) goto 903
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
            if (ierr.ne.0) goto 903
            call clscdf(cstid,ierr)
            if (ierr.ne.0) goto 903
 
         else
            varname = fieldname
            call getdef(-fid,varname,ndim,mdv,vardim,
     >                  varmin,varmax,stag,ierr)
            if (ierr.ne.0) goto 900
            
         endif
 
c        Set the grid dimensions and constants - vardim(3) is taken from constants file
         nx   = vardim(1)
         ny   = vardim(2)
         nz   = vardim(3)
         xmin = varmin(1)
         ymin = varmin(2)
         xmax = varmax(1)
         ymax = varmax(2)
         dx   = (xmax-xmin)/real(nx-1)
         dy   = (ymax-ymin)/real(ny-1)
 
c        Get pole position - if no constants file available, set default pole
         call getcfn(-fid,cstfile,ierr)
         if (ierr.eq.0) then  
            call cdfopn(cstfile,cstid,ierr)
            if (ierr.ne.0) goto 903
            call getpole(cstid,pollon,pollat,ierr)
            if (ierr.ne.0) goto 903
            call clscdf(cstid,ierr)
            if (ierr.ne.0) goto 903
         else
            pollon = 0.
            pollat = 90.
         endif
 
c     Get non-constant grid parameters (surface pressure and vertical grid)
      else
         
c        Special handling for fieldname 'P.ML' -> in this case the fields
c        P and PS are available on the data file and can be read in. There
c        is no need to reconstruct it from PS,AK and BK. This mode is
c        used for model-level (2D) trajectories
         if ( fieldname.eq.'P.ML' ) then
 
c           Get the right time to read
            call gettimes(fid,times,ntimes,ierr)
            if (ierr.ne.0) goto 901
            isok=0
            do i=1,ntimes
               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
 
c           Read surface pressure and 3D pressure
            varname='PS'
            call getdat(fid,varname,rtime,0,ps,ierr)
            if (ierr.ne.0) goto 904
            varname='P'
            call getdat(fid,varname,rtime,0,p3,ierr)
            if (ierr.ne.0) goto 904
 
c           Set MDV to 1050. - otherwise interpolation routines don't work
            do i=1,nx
              do j=1,ny
                do k=1,nz
                   if ( p3(i,j,k).lt.0. ) p3(i,j,k) = 1050.
                enddo
              enddo
            enddo
 
c           Don't care about other stuff - finish subroutine
            goto 800
 
         endif
 
c        Get full grid info - in particular staggering flag; set flag for 2D tracing
         if ( fieldname.eq.'PLEV' ) then
            varname = 'PS'
            stagz   = 0.
            call getdef(fid,varname,ndim,mdv,vardim,
     >                  varmin,varmax,stag,ierr)
            if (ierr.ne.0) goto 900
            call getcfn(fid,cstfile,ierr)
            if (ierr.ne.0) goto 903
            call cdfopn(cstfile,cstid,ierr)
            if (ierr.ne.0) goto 903
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
            if (ierr.ne.0) goto 903
            call clscdf(cstid,ierr)
            if (ierr.ne.0) goto 903
            
         elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then       
            varname = 'PS'
            stagz   = -0.5
            call getdef(fid,varname,ndim,mdv,vardim,
     >                  varmin,varmax,stag,ierr)
            if (ierr.ne.0) goto 900
            call getcfn(fid,cstfile,ierr)
            if (ierr.ne.0) goto 903
            call cdfopn(cstfile,cstid,ierr)
            if (ierr.ne.0) goto 903
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
            if (ierr.ne.0) goto 903
            call clscdf(cstid,ierr)
            if (ierr.ne.0) goto 903
 
         else         
            varname=fieldname
            call getdef(fid,varname,ndim,mdv,vardim,
     >                  varmin,varmax,stag,ierr)
            if (ierr.ne.0) goto 900
            if (vardim(3).eq.1) is2d = 1
         endif
 
c        Get time information (check if time is correct)
         call gettimes(fid,times,ntimes,ierr)
         if (ierr.ne.0) goto 901
         isok=0
         do i=1,ntimes
            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) goto 905
 
c        If 2D tracing requested: take dummay values for PS, AKLEV,BKLEV,AKLAY,BKLAY
         if ( is2d.eq.1 ) then
            
            do i=1,nx
               do j=1,ny
                  ps(i,j) = 1050.
               enddo
            enddo
            
            do k=1,nz
               aklev(k) = 0.
               bklev(k) = real(nz-k)/real(nz-1) 
               aklay(k) = 0.
               bklay(k) = real(nz-k)/real(nz-1) 
            enddo
 
c        3D tracing - read PS, AKLEV,BKLEV,AKLAY;BKLAY
         else
 
c           Read the level coefficients from the constants file
            call getcfn(fid,cstfile,ierr)
            if (ierr.ne.0) goto 903
            call cdfopn(cstfile,cstid,ierr)
            if (ierr.ne.0) goto 903
            call getlevs(cstid,vardim(3),aklev,bklev,aklay,bklay,ierr)
            if (ierr.ne.0) goto 903
            call clscdf(cstid,ierr)
            if (ierr.ne.0) goto 903
 
c           Check whether PS is needed to get the 3d pressure field
            plev = 1
            do i=1,nz
              if ( (abs(stagz).lt.eps).and.(abs(bklev(i)).gt.eps) ) then
                plev = 0
              endif
              if ( (abs(stagz).gt.eps).and.(abs(bklay(i)).gt.eps) ) then
                plev = 0
              endif
            enddo
 
c           Read surface pressure if needed
            if ( plev.ne.1 ) then
              varname='PS'
              call getdat(fid,varname,rtime,0,ps,ierr)
              if (ierr.ne.0) goto 904
            else
              do i=1,nx
                do j=1,ny
                    ps(i,j) = 1000.
                enddo
              enddo
            endif
 
         endif
 
c        Calculate layer and level pressures
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  if ( abs(stagz).lt.eps ) then
                     p3(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
                  else
                     p3(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
                  endif
               enddo
            enddo
         enddo
 
c        Set the ak and bk for the vertical grid
         do k=1,nz
            if ( abs(stagz).lt.eps ) then
               ak(k)=aklev(k)
               bk(k)=bklev(k)
            else
               ak(k)=aklay(k)
               bk(k)=bklay(k)
            endif
         enddo
 
      endif
      
c     Exit point for subroutine
 800  continue
      return
      
c     Exception handling
 900  print*,'Cannot retrieve grid dimension from ',fid
      stop
 901  print*,'Cannot retrieve grid parameters from ',fid
      stop
 902  print*,'Grid inconsistency detected for ',fid
      stop
 903  print*,'Problem with level coefficients from ',trim(cstfile)
      stop
 904  print*,'Cannot read surface pressure from ',trim(cstfile)
      stop
 905  print*,'Cannot find time ',time,' on ',fid
      stop
 906  print*,'Unable to get grid info ',fid
      stop
      
      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>.
 
      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)
 
c     Auxiliary variables
      integer      isok
      integer      i,j,k
      integer      nz1
      real         rtime
 
c     Read variable definition - for P, PLEV and PLAY: load also ak,bk
      if ( ( fieldname.eq.'PLEV' ).or.
     >     ( fieldname.eq.'PLAY' ).or.
     >     ( fieldname.eq.'P'    ) )
     >then
         call getcfn(fid,cstfile,ierr)
         if (ierr.ne.0) goto 905
         call cdfopn(cstfile,cstid,ierr)
         if (ierr.ne.0) goto 905
         call getlevs(cstid,nz1,aklev,bklev,aklay,bklay,ierr)
         if (ierr.ne.0) goto 905
         call clscdf(cstid,ierr)
         if (ierr.ne.0) goto 905
         varname = 'PS'
         call getdef(fid,varname,ndim,mdv,vardim,
     >               varmin,varmax,stag,ierr)
         vardim(3) = nz1
         if (ierr.ne.0) goto 900
 
      else
         varname = fieldname
         call getdef(fid,varname,ndim,mdv,vardim,
     >               varmin,varmax,stag,ierr)
         if (ierr.ne.0) goto 900
         stagz=stag(3)
      endif
 
c     Get time information (set time to first one in the file)
      call gettimes(fid,times,ntimes,ierr)
      if (ierr.ne.0) goto 902
      isok=0
      do i=1,ntimes
         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 ) goto 904
 
c     Read  field
      if ( ( fieldname.eq.'P' ).or.(fieldname.eq.'PLAY') )  then       ! P, PLAY
         stagz   = -0.5
         varname = 'PS'
         call getdat(fid,varname,rtime,0,ps,ierr)
         if (ierr.ne.0) goto 903
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  field(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
               enddo
            enddo
         enddo
         
      elseif ( fieldname.eq.'PLEV' )  then                             ! PLEV
         stagz   = 0.
         varname = 'PS'
         call getdat(fid,varname,rtime,0,ps,ierr)
         if (ierr.ne.0) goto 903
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  field(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
               enddo
            enddo
         enddo
 
      else                                                             ! Other fields
         varname=fieldname
         call getdat(fid,varname,rtime,0,field,ierr)
         if (ierr.ne.0) goto 903
 
      endif
 
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
      if ( vardim(3).eq.1 ) then
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  field(i,j,k) = field(i,j,1)
               enddo
            enddo
         enddo
      endif
         
 
 
c     Exception handling
      return
      
 900  print*,'Cannot retrieve definition for ',trim(varname),'  ',fid
      stop
 901  print*,'Grid inconsistency detected for ',trim(varname),'  ',fid
      stop
 902  print*,'Cannot retrieve time for ',trim(varname),'  ',fid
      stop
 903  print*,'Cannot load wind component ',trim(varname),'  ',fid
      stop
 904  print*,'Cannot load time ',time,' for ',trim(varname),'  ',fid
      stop
 905  print*,'Cannot load time vertical grid AK, BK from file  ',fid
      stop
      
      end
 
c     ------------------------------------------------------------
c     Close input file
c     ------------------------------------------------------------
 
      subroutine input_close(fid)
 
c     Close the input file with file identifier <fid>.
 
      implicit none
 
c     Declaration of subroutine parameters
      integer fid
 
c     Auxiliary variables
      integer ierr
 
c     Close file
      call clscdf(fid,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
 
      implicit none
 
c     Declaration of subroutine parameters
      integer      fid
      integer      nvars
      character*80 vnam(200)
 
c     Auxiliary variables
      integer ierr
      
c     Get list and save      
      call getvars(fid,nvars,vnam,ierr)
 
      end

Rev 44	Rev 46
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	`implicit none`	34	`implicit none`
35		35
36	`c Declaration of subroutine parameters`	36	`c Declaration of subroutine parameters`
37	`integer fid ! File identifier`	37	`integer fid ! File identifier`
38	`character*80 filename ! Filename`	38	`character*80 filename ! Filename`
39		39
40	`c Declaration of auxiliary variables`	40	`c Declaration of auxiliary variables`
41	`integer ierr`	41	`integer ierr`
42		42
43	`c Open IVE netcdf file`	43	`c Open IVE netcdf file`
44	`call cdfopn(filename,fid,ierr)`	44	`call cdfopn(filename,fid,ierr)`
45	`if (ierr.ne.0) goto 900`	45	`if (ierr.ne.0) goto 900`
46		46
47	`c Exception handling`	47	`c Exception handling`
48	`return`	48	`return`
49		49
50	`900 print*,'Cannot open input file ',trim(filename)`	50	`900 print*,'Cannot open input file ',trim(filename)`
51	`stop`	51	`stop`
52		52
53	`end`	53	`end`
54		54
55		55
56	`c ------------------------------------------------------------`	56	`c ------------------------------------------------------------`
57	`c Read information about the grid`	57	`c Read information about the grid`
58	`c ------------------------------------------------------------`	58	`c ------------------------------------------------------------`
59		59
60	`subroutine input_grid`	60	`subroutine input_grid`
61	`> (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`	61	`> (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`
62	`> time,pollon,pollat,p3,ps,nz,ak,bk,stagz,`	62	`> time,pollon,pollat,p3,ps,nz,ak,bk,stagz,`
63	`> timecheck)`	63	`> timecheck)`
64		64
65	`c Read grid information at <time> from file with identifier <fid>.`	65	`c Read grid information at <time> from file with identifier <fid>.`
66	`c The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>`	66	`c The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>`
67	`c with pole position at <pollon,pollat> and grid dimension <nx,ny>.`	67	`c with pole position at <pollon,pollat> and grid dimension <nx,ny>.`
68	`c The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either`	68	`c The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either`
69	`c on the staggered or unstaggered grid (with <stagz> as the flag).`	69	`c on the staggered or unstaggered grid (with <stagz> as the flag).`
70	`c The surface pressure is given in <ps(nx,ny)>. If <fid> is negative,`	70	`c The surface pressure is given in <ps(nx,ny)>. If <fid> is negative,`
71	`c only the grid dimensions and grid parameters (xmin...pollat,nz) are`	71	`c only the grid dimensions and grid parameters (xmin...pollat,nz) are`
72	`c determined and returned (this is needed for dynamical allocation of`	72	`c determined and returned (this is needed for dynamical allocation of`
73	`c memory).`	73	`c memory).`
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 times(10)`	103	`real times(10)`
104	`integer ntimes`	104	`integer ntimes`
105	`real aklay(nz),bklay(nz),aklev(nz),bklev(nz)`	105	`real aklay(nz),bklay(nz),aklev(nz),bklev(nz)`
106	`integer nvars`	106	`integer nvars`
107	`character*80 vars(100)`	107	`character*80 vars(100)`
108		108
109	`c Auxiliary varaibles`	109	`c Auxiliary varaibles`
110	`integer ierr`	110	`integer ierr`
111	`integer i,j,k`	111	`integer i,j,k`
112	`integer isok`	112	`integer isok`
113	`real tmp(200)`	113	`real tmp(200)`
114	`character*80 varname`	114	`character*80 varname`
115	`real rtime`	115	`real rtime`
116	`integer is2d`	116	`integer is2d`
117	`integer plev`	117	`integer plev`
118		118
119	`c Init the flag for 2D variables - assume a 3D field`	119	`c Init the flag for 2D variables - assume a 3D field`
120	`is2d = 0`	120	`is2d = 0`
121		121
122	`c Init the flag for pressure levels (PS is not needed)`	122	`c Init the flag for pressure levels (PS is not needed)`
123	`plev = 0`	123	`plev = 0`
124		124
125	`c Inquire dimensions and grid constants if <fid> is negative`	125	`c Inquire dimensions and grid constants if <fid> is negative`
126	`if (fid.lt.0) then`	126	`if (fid.lt.0) then`
127		127
128	`c Get grid info for the selected variable`	128	`c Get grid info for the selected variable`
129	`if ( fieldname.eq.'PLEV' ) then`	129	`if ( fieldname.eq.'PLEV' ) then`
130	`varname = 'PS'`	130	`varname = 'PS'`
131	`stagz = 0.`	131	`stagz = 0.`
132	`call getdef(-fid,varname,ndim,mdv,vardim,`	132	`call getdef(-fid,varname,ndim,mdv,vardim,`
133	`> varmin,varmax,stag,ierr)`	133	`> varmin,varmax,stag,ierr)`
134	`if (ierr.ne.0) goto 900`	134	`if (ierr.ne.0) goto 900`
135	`call getcfn(-fid,cstfile,ierr)`	135	`call getcfn(-fid,cstfile,ierr)`
136	`if (ierr.ne.0) goto 903`	136	`if (ierr.ne.0) goto 903`
137	`call cdfopn(cstfile,cstid,ierr)`	137	`call cdfopn(cstfile,cstid,ierr)`
138	`if (ierr.ne.0) goto 903`	138	`if (ierr.ne.0) goto 903`
139	`call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)`	139	`call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)`
140	`if (ierr.ne.0) goto 903`	140	`if (ierr.ne.0) goto 903`
141	`call clscdf(cstid,ierr)`	141	`call clscdf(cstid,ierr)`
142	`if (ierr.ne.0) goto 903`	142	`if (ierr.ne.0) goto 903`
143		143
144	`elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then`	144	`elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then`
145	`varname = 'PS'`	145	`varname = 'PS'`
146	`stagz = -0.5`	146	`stagz = -0.5`
147	`call getdef(-fid,varname,ndim,mdv,vardim,`	147	`call getdef(-fid,varname,ndim,mdv,vardim,`
148	`> varmin,varmax,stag,ierr)`	148	`> varmin,varmax,stag,ierr)`
149	`if (ierr.ne.0) goto 900`	149	`if (ierr.ne.0) goto 900`
150	`call getcfn(-fid,cstfile,ierr)`	150	`call getcfn(-fid,cstfile,ierr)`
151	`if (ierr.ne.0) goto 903`	151	`if (ierr.ne.0) goto 903`
152	`call cdfopn(cstfile,cstid,ierr)`	152	`call cdfopn(cstfile,cstid,ierr)`
153	`if (ierr.ne.0) goto 903`	153	`if (ierr.ne.0) goto 903`
154	`call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)`	154	`call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)`
155	`if (ierr.ne.0) goto 903`	155	`if (ierr.ne.0) goto 903`
156	`call clscdf(cstid,ierr)`	156	`call clscdf(cstid,ierr)`
157	`if (ierr.ne.0) goto 903`	157	`if (ierr.ne.0) goto 903`
158		158
159	`else`	159	`else`
160	`varname = fieldname`	160	`varname = fieldname`
161	`call getdef(-fid,varname,ndim,mdv,vardim,`	161	`call getdef(-fid,varname,ndim,mdv,vardim,`
162	`> varmin,varmax,stag,ierr)`	162	`> varmin,varmax,stag,ierr)`
163	`if (ierr.ne.0) goto 900`	163	`if (ierr.ne.0) goto 900`
164		164
165	`endif`	165	`endif`
166		166
167	`c Set the grid dimensions and constants - vardim(3) is taken from constants file`	167	`c Set the grid dimensions and constants - vardim(3) is taken from constants file`
168	`nx = vardim(1)`	168	`nx = vardim(1)`
169	`ny = vardim(2)`	169	`ny = vardim(2)`
170	`nz = vardim(3)`	170	`nz = vardim(3)`
171	`xmin = varmin(1)`	171	`xmin = varmin(1)`
172	`ymin = varmin(2)`	172	`ymin = varmin(2)`
173	`xmax = varmax(1)`	173	`xmax = varmax(1)`
174	`ymax = varmax(2)`	174	`ymax = varmax(2)`
175	`dx = (xmax-xmin)/real(nx-1)`	175	`dx = (xmax-xmin)/real(nx-1)`
176	`dy = (ymax-ymin)/real(ny-1)`	176	`dy = (ymax-ymin)/real(ny-1)`
177		177
178	`c Get pole position - if no constants file available, set default pole`	178	`c Get pole position - if no constants file available, set default pole`
179	`call getcfn(-fid,cstfile,ierr)`	179	`call getcfn(-fid,cstfile,ierr)`
180	`if (ierr.eq.0) then`	180	`if (ierr.eq.0) then`
181	`call cdfopn(cstfile,cstid,ierr)`	181	`call cdfopn(cstfile,cstid,ierr)`
182	`if (ierr.ne.0) goto 903`	182	`if (ierr.ne.0) goto 903`
183	`call getpole(cstid,pollon,pollat,ierr)`	183	`call getpole(cstid,pollon,pollat,ierr)`
184	`if (ierr.ne.0) goto 903`	184	`if (ierr.ne.0) goto 903`
185	`call clscdf(cstid,ierr)`	185	`call clscdf(cstid,ierr)`
186	`if (ierr.ne.0) goto 903`	186	`if (ierr.ne.0) goto 903`
187	`else`	187	`else`
188	`pollon = 0.`	188	`pollon = 0.`
189	`pollat = 90.`	189	`pollat = 90.`
190	`endif`	190	`endif`
191		191
192	`c Get non-constant grid parameters (surface pressure and vertical grid)`	192	`c Get non-constant grid parameters (surface pressure and vertical grid)`
193	`else`	193	`else`
194		194
195	`c Special handling for fieldname 'P.ML' -> in this case the fields`	195	`c Special handling for fieldname 'P.ML' -> in this case the fields`
196	`c P and PS are available on the data file and can be read in. There`	196	`c P and PS are available on the data file and can be read in. There`
197	`c is no need to reconstruct it from PS,AK and BK. This mode is`	197	`c is no need to reconstruct it from PS,AK and BK. This mode is`
198	`c used for model-level (2D) trajectories`	198	`c used for model-level (2D) trajectories`
199	`if ( fieldname.eq.'P.ML' ) then`	199	`if ( fieldname.eq.'P.ML' ) then`
200		200
201	`c Get the right time to read`	201	`c Get the right time to read`
202	`call gettimes(fid,times,ntimes,ierr)`	202	`call gettimes(fid,times,ntimes,ierr)`
203	`if (ierr.ne.0) goto 901`	203	`if (ierr.ne.0) goto 901`
204	`isok=0`	204	`isok=0`
205	`do i=1,ntimes`	205	`do i=1,ntimes`
206	`if (abs(time-times(i)).lt.eps) then`	206	`if (abs(time-times(i)).lt.eps) then`
207	`isok = 1`	207	`isok = 1`
208	`rtime = times(i)`	208	`rtime = times(i)`
209	`elseif (timecheck.eq.'no') then`	209	`elseif (timecheck.eq.'no') then`
210	`isok = 1`	210	`isok = 1`
211	`rtime = times(1)`	211	`rtime = times(1)`
212	`endif`	212	`endif`
213	`enddo`	213	`enddo`
214		214
215	`c Read surface pressure and 3D pressure`	215	`c Read surface pressure and 3D pressure`
216	`varname='PS'`	216	`varname='PS'`
217	`call getdat(fid,varname,rtime,0,ps,ierr)`	217	`call getdat(fid,varname,rtime,0,ps,ierr)`
218	`if (ierr.ne.0) goto 904`	218	`if (ierr.ne.0) goto 904`
219	`varname='P'`	219	`varname='P'`
220	`call getdat(fid,varname,rtime,0,p3,ierr)`	220	`call getdat(fid,varname,rtime,0,p3,ierr)`
221	`if (ierr.ne.0) goto 904`	221	`if (ierr.ne.0) goto 904`
222		222
223	`c Set MDV to 1050. - otherwise interpolation routines don't work`	223	`c Set MDV to 1050. - otherwise interpolation routines don't work`
224	`do i=1,nx`	224	`do i=1,nx`
225	`do j=1,ny`	225	`do j=1,ny`
226	`do k=1,nz`	226	`do k=1,nz`
227	`if ( p3(i,j,k).lt.0. ) p3(i,j,k) = 1050.`	227	`if ( p3(i,j,k).lt.0. ) p3(i,j,k) = 1050.`
228	`enddo`	228	`enddo`
229	`enddo`	229	`enddo`
230	`enddo`	230	`enddo`
231		231
232	`c Don't care about other stuff - finish subroutine`	232	`c Don't care about other stuff - finish subroutine`
233	`goto 800`	233	`goto 800`
234		234
235	`endif`	235	`endif`
236		236
237	`c Get full grid info - in particular staggering flag; set flag for 2D tracing`	237	`c Get full grid info - in particular staggering flag; set flag for 2D tracing`
238	`if ( fieldname.eq.'PLEV' ) then`	238	`if ( fieldname.eq.'PLEV' ) then`
239	`varname = 'PS'`	239	`varname = 'PS'`
240	`stagz = 0.`	240	`stagz = 0.`
241	`call getdef(fid,varname,ndim,mdv,vardim,`	241	`call getdef(fid,varname,ndim,mdv,vardim,`
242	`> varmin,varmax,stag,ierr)`	242	`> varmin,varmax,stag,ierr)`
243	`if (ierr.ne.0) goto 900`	243	`if (ierr.ne.0) goto 900`
244	`call getcfn(fid,cstfile,ierr)`	244	`call getcfn(fid,cstfile,ierr)`
245	`if (ierr.ne.0) goto 903`	245	`if (ierr.ne.0) goto 903`
246	`call cdfopn(cstfile,cstid,ierr)`	246	`call cdfopn(cstfile,cstid,ierr)`
247	`if (ierr.ne.0) goto 903`	247	`if (ierr.ne.0) goto 903`
248	`call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)`	248	`call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)`
249	`if (ierr.ne.0) goto 903`	249	`if (ierr.ne.0) goto 903`
250	`call clscdf(cstid,ierr)`	250	`call clscdf(cstid,ierr)`
251	`if (ierr.ne.0) goto 903`	251	`if (ierr.ne.0) goto 903`
252		252
253	`elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then`	253	`elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then`
254	`varname = 'PS'`	254	`varname = 'PS'`
255	`stagz = -0.5`	255	`stagz = -0.5`
256	`call getdef(fid,varname,ndim,mdv,vardim,`	256	`call getdef(fid,varname,ndim,mdv,vardim,`
257	`> varmin,varmax,stag,ierr)`	257	`> varmin,varmax,stag,ierr)`
258	`if (ierr.ne.0) goto 900`	258	`if (ierr.ne.0) goto 900`
259	`call getcfn(fid,cstfile,ierr)`	259	`call getcfn(fid,cstfile,ierr)`
260	`if (ierr.ne.0) goto 903`	260	`if (ierr.ne.0) goto 903`
261	`call cdfopn(cstfile,cstid,ierr)`	261	`call cdfopn(cstfile,cstid,ierr)`
262	`if (ierr.ne.0) goto 903`	262	`if (ierr.ne.0) goto 903`
263	`call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)`	263	`call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)`
264	`if (ierr.ne.0) goto 903`	264	`if (ierr.ne.0) goto 903`
265	`call clscdf(cstid,ierr)`	265	`call clscdf(cstid,ierr)`
266	`if (ierr.ne.0) goto 903`	266	`if (ierr.ne.0) goto 903`
267		267
268	`else`	268	`else`
269	`varname=fieldname`	269	`varname=fieldname`
270	`call getdef(fid,varname,ndim,mdv,vardim,`	270	`call getdef(fid,varname,ndim,mdv,vardim,`
271	`> varmin,varmax,stag,ierr)`	271	`> varmin,varmax,stag,ierr)`
272	`if (ierr.ne.0) goto 900`	272	`if (ierr.ne.0) goto 900`
273	`if (vardim(3).eq.1) is2d = 1`	273	`if (vardim(3).eq.1) is2d = 1`
274	`endif`	274	`endif`
275		275
276	`c Get time information (check if time is correct)`	276	`c Get time information (check if time is correct)`
277	`call gettimes(fid,times,ntimes,ierr)`	277	`call gettimes(fid,times,ntimes,ierr)`
278	`if (ierr.ne.0) goto 901`	278	`if (ierr.ne.0) goto 901`
279	`isok=0`	279	`isok=0`
280	`do i=1,ntimes`	280	`do i=1,ntimes`
281	`if (abs(time-times(i)).lt.eps) then`	281	`if (abs(time-times(i)).lt.eps) then`
282	`isok = 1`	282	`isok = 1`
283	`rtime = times(i)`	283	`rtime = times(i)`
284	`elseif (timecheck.eq.'no') then`	284	`elseif (timecheck.eq.'no') then`
285	`isok = 1`	285	`isok = 1`
286	`rtime = times(1)`	286	`rtime = times(1)`
287	`endif`	287	`endif`
288	`enddo`	288	`enddo`
289	`if ( isok.eq.0) goto 905`	289	`if ( isok.eq.0) goto 905`
290		290
291	`c If 2D tracing requested: take dummay values for PS, AKLEV,BKLEV,AKLAY,BKLAY`	291	`c If 2D tracing requested: take dummay values for PS, AKLEV,BKLEV,AKLAY,BKLAY`
292	`if ( is2d.eq.1 ) then`	292	`if ( is2d.eq.1 ) then`
293		293
294	`do i=1,nx`	294	`do i=1,nx`
295	`do j=1,ny`	295	`do j=1,ny`
296	`ps(i,j) = 1050.`	296	`ps(i,j) = 1050.`
297	`enddo`	297	`enddo`
298	`enddo`	298	`enddo`
299		299
300	`do k=1,nz`	300	`do k=1,nz`
301	`aklev(k) = 0.`	301	`aklev(k) = 0.`
302	`bklev(k) = real(nz-k)/real(nz-1)`	302	`bklev(k) = real(nz-k)/real(nz-1)`
303	`aklay(k) = 0.`	303	`aklay(k) = 0.`
304	`bklay(k) = real(nz-k)/real(nz-1)`	304	`bklay(k) = real(nz-k)/real(nz-1)`
305	`enddo`	305	`enddo`
306		306
307	`c 3D tracing - read PS, AKLEV,BKLEV,AKLAY;BKLAY`	307	`c 3D tracing - read PS, AKLEV,BKLEV,AKLAY;BKLAY`
308	`else`	308	`else`
309		309
310	`c Read the level coefficients from the constants file`	310	`c Read the level coefficients from the constants file`
311	`call getcfn(fid,cstfile,ierr)`	311	`call getcfn(fid,cstfile,ierr)`
312	`if (ierr.ne.0) goto 903`	312	`if (ierr.ne.0) goto 903`
313	`call cdfopn(cstfile,cstid,ierr)`	313	`call cdfopn(cstfile,cstid,ierr)`
314	`if (ierr.ne.0) goto 903`	314	`if (ierr.ne.0) goto 903`
315	`call getlevs(cstid,vardim(3),aklev,bklev,aklay,bklay,ierr)`	315	`call getlevs(cstid,vardim(3),aklev,bklev,aklay,bklay,ierr)`
316	`if (ierr.ne.0) goto 903`	316	`if (ierr.ne.0) goto 903`
317	`call clscdf(cstid,ierr)`	317	`call clscdf(cstid,ierr)`
318	`if (ierr.ne.0) goto 903`	318	`if (ierr.ne.0) goto 903`
319		319
320	`c Check whether PS is needed to get the 3d pressure field`	320	`c Check whether PS is needed to get the 3d pressure field`
321	`plev = 1`	321	`plev = 1`
322	`do i=1,nz`	322	`do i=1,nz`
323	`if ( (abs(stagz).lt.eps).and.(abs(bklev(i)).gt.eps) ) then`	323	`if ( (abs(stagz).lt.eps).and.(abs(bklev(i)).gt.eps) ) then`
324	`plev = 0`	324	`plev = 0`
325	`endif`	325	`endif`
326	`if ( (abs(stagz).gt.eps).and.(abs(bklay(i)).gt.eps) ) then`	326	`if ( (abs(stagz).gt.eps).and.(abs(bklay(i)).gt.eps) ) then`
327	`plev = 0`	327	`plev = 0`
328	`endif`	328	`endif`
329	`enddo`	329	`enddo`
330		330
331	`c Read surface pressure if needed`	331	`c Read surface pressure if needed`
332	`if ( plev.ne.1 ) then`	332	`if ( plev.ne.1 ) then`
333	`varname='PS'`	333	`varname='PS'`
334	`call getdat(fid,varname,rtime,0,ps,ierr)`	334	`call getdat(fid,varname,rtime,0,ps,ierr)`
335	`if (ierr.ne.0) goto 904`	335	`if (ierr.ne.0) goto 904`
336	`else`	336	`else`
337	`do i=1,nx`	337	`do i=1,nx`
338	`do j=1,ny`	338	`do j=1,ny`
339	`ps(i,j) = 1000.`	339	`ps(i,j) = 1000.`
340	`enddo`	340	`enddo`
341	`enddo`	341	`enddo`
342	`endif`	342	`endif`
343		343
344	`endif`	344	`endif`
345		345
346	`c Calculate layer and level pressures`	346	`c Calculate layer and level pressures`
347	`do i=1,nx`	347	`do i=1,nx`
348	`do j=1,ny`	348	`do j=1,ny`
349	`do k=1,nz`	349	`do k=1,nz`
350	`if ( abs(stagz).lt.eps ) then`	350	`if ( abs(stagz).lt.eps ) then`
351	`p3(i,j,k)=aklev(k)+bklev(k)*ps(i,j)`	351	`p3(i,j,k)=aklev(k)+bklev(k)*ps(i,j)`
352	`else`	352	`else`
353	`p3(i,j,k)=aklay(k)+bklay(k)*ps(i,j)`	353	`p3(i,j,k)=aklay(k)+bklay(k)*ps(i,j)`
354	`endif`	354	`endif`
355	`enddo`	355	`enddo`
356	`enddo`	356	`enddo`
357	`enddo`	357	`enddo`
358		358
359	`c Set the ak and bk for the vertical grid`	359	`c Set the ak and bk for the vertical grid`
360	`do k=1,nz`	360	`do k=1,nz`
361	`if ( abs(stagz).lt.eps ) then`	361	`if ( abs(stagz).lt.eps ) then`
362	`ak(k)=aklev(k)`	362	`ak(k)=aklev(k)`
363	`bk(k)=bklev(k)`	363	`bk(k)=bklev(k)`
364	`else`	364	`else`
365	`ak(k)=aklay(k)`	365	`ak(k)=aklay(k)`
366	`bk(k)=bklay(k)`	366	`bk(k)=bklay(k)`
367	`endif`	367	`endif`
368	`enddo`	368	`enddo`
369		369
370	`endif`	370	`endif`
371		371
372	`c Exit point for subroutine`	372	`c Exit point for subroutine`
373	`800 continue`	373	`800 continue`
374	`return`	374	`return`
375		375
376	`c Exception handling`	376	`c Exception handling`
377	`900 print*,'Cannot retrieve grid dimension from ',fid`	377	`900 print*,'Cannot retrieve grid dimension from ',fid`
378	`stop`	378	`stop`
379	`901 print*,'Cannot retrieve grid parameters from ',fid`	379	`901 print*,'Cannot retrieve grid parameters from ',fid`
380	`stop`	380	`stop`
381	`902 print*,'Grid inconsistency detected for ',fid`	381	`902 print*,'Grid inconsistency detected for ',fid`
382	`stop`	382	`stop`
383	`903 print*,'Problem with level coefficients from ',trim(cstfile)`	383	`903 print*,'Problem with level coefficients from ',trim(cstfile)`
384	`stop`	384	`stop`
385	`904 print*,'Cannot read surface pressure from ',trim(cstfile)`	385	`904 print*,'Cannot read surface pressure from ',trim(cstfile)`
386	`stop`	386	`stop`
387	`905 print*,'Cannot find time ',time,' on ',fid`	387	`905 print*,'Cannot find time ',time,' on ',fid`
388	`stop`	388	`stop`
389	`906 print*,'Unable to get grid info ',fid`	389	`906 print*,'Unable to get grid info ',fid`
390	`stop`	390	`stop`
391		391
392	`end`	392	`end`
393		393
394	`c ------------------------------------------------------------`	394	`c ------------------------------------------------------------`
395	`c Read wind information`	395	`c Read wind information`
396	`c ------------------------------------------------------------`	396	`c ------------------------------------------------------------`
397		397
398	`subroutine input_wind (fid,fieldname,field,time,stagz,mdv,`	398	`subroutine input_wind (fid,fieldname,field,time,stagz,mdv,`
399	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,`	399	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,`
400	`> timecheck)`	400	`> timecheck)`
401		401
402	`c Read the wind component <fieldname> from the file with identifier`	402	`c Read the wind component <fieldname> from the file with identifier`
403	`c <fid> and save it in the 3d array <field>. The vertical staggering`	403	`c <fid> and save it in the 3d array <field>. The vertical staggering`
404	`c information is provided in <stagz> and gives the reference to either`	404	`c information is provided in <stagz> and gives the reference to either`
405	`c the layer or level field from <input_grid>. A consistency check is`	405	`c the layer or level field from <input_grid>. A consistency check is`
406	`c performed to have an agreement with the grid specified by <xmin,xmax,`	406	`c performed to have an agreement with the grid specified by <xmin,xmax,`
407	`c ymin,ymax,dx,dy,nx,ny,nz>.`	407	`c ymin,ymax,dx,dy,nx,ny,nz>.`
408		408
409	`implicit none`	409	`implicit none`
410		410
411	`c Declaration of variables and parameters`	411	`c Declaration of variables and parameters`
412	`integer fid ! File identifier`	412	`integer fid ! File identifier`
413	`character*80 fieldname ! Name of the wind field`	413	`character*80 fieldname ! Name of the wind field`
414	`integer nx,ny,nz ! Dimension of fields`	414	`integer nx,ny,nz ! Dimension of fields`
415	`real field(nx,ny,nz) ! 3d wind field`	415	`real field(nx,ny,nz) ! 3d wind field`
416	`real stagz ! Staggering in the z direction`	416	`real stagz ! Staggering in the z direction`
417	`real mdv ! Missing data flag`	417	`real mdv ! Missing data flag`
418	`real xmin,xmax,ymin,ymax ! Domain size`	418	`real xmin,xmax,ymin,ymax ! Domain size`
419	`real dx,dy ! Horizontal resolution`	419	`real dx,dy ! Horizontal resolution`
420	`real time ! Time`	420	`real time ! Time`
421	`character*80 timecheck ! Either 'yes' or 'no'`	421	`character*80 timecheck ! Either 'yes' or 'no'`
422		422
423	`c Numerical and physical parameters`	423	`c Numerical and physical parameters`
424	`real eps ! Numerical epsilon`	424	`real eps ! Numerical epsilon`
425	`parameter (eps=0.001)`	425	`parameter (eps=0.001)`
426	`real notimecheck ! 'Flag' for no time check`	426	`real notimecheck ! 'Flag' for no time check`
427	`parameter (notimecheck=7.26537)`	427	`parameter (notimecheck=7.26537)`
428		428
429	`c Netcdf variables`	429	`c Netcdf variables`
430	`integer ierr`	430	`integer ierr`
431	`character*80 varname`	431	`character*80 varname`
432	`integer vardim(4)`	432	`integer vardim(4)`
433	`real varmin(4),varmax(4)`	433	`real varmin(4),varmax(4)`
434	`real stag(4)`	434	`real stag(4)`
435	`integer ndim`	435	`integer ndim`
436	`real times(10)`	436	`real times(10)`
437	`integer ntimes`	437	`integer ntimes`
438	`character*80 cstfile`	438	`character*80 cstfile`
439	`integer cstid`	439	`integer cstid`
440	`real aklay(200),bklay(200),aklev(200),bklev(200)`	440	`real aklay(200),bklay(200),aklev(200),bklev(200)`
441	`real ps(nx,ny)`	441	`real ps(nx,ny)`
442		442
443	`c Auxiliary variables`	443	`c Auxiliary variables`
444	`integer isok`	444	`integer isok`
445	`integer i,j,k`	445	`integer i,j,k`
446	`integer nz1`	446	`integer nz1`
447	`real rtime`	447	`real rtime`
448		448
449	`c Read variable definition - for P, PLEV and PLAY: load also ak,bk`	449	`c Read variable definition - for P, PLEV and PLAY: load also ak,bk`
450	`if ( ( fieldname.eq.'PLEV' ).or.`	450	`if ( ( fieldname.eq.'PLEV' ).or.`
451	`> ( fieldname.eq.'PLAY' ).or.`	451	`> ( fieldname.eq.'PLAY' ).or.`
452	`> ( fieldname.eq.'P' ) )`	452	`> ( fieldname.eq.'P' ) )`
453	`>then`	453	`>then`
454	`call getcfn(fid,cstfile,ierr)`	454	`call getcfn(fid,cstfile,ierr)`
455	`if (ierr.ne.0) goto 905`	455	`if (ierr.ne.0) goto 905`
456	`call cdfopn(cstfile,cstid,ierr)`	456	`call cdfopn(cstfile,cstid,ierr)`
457	`if (ierr.ne.0) goto 905`	457	`if (ierr.ne.0) goto 905`
458	`call getlevs(cstid,nz1,aklev,bklev,aklay,bklay,ierr)`	458	`call getlevs(cstid,nz1,aklev,bklev,aklay,bklay,ierr)`
459	`if (ierr.ne.0) goto 905`	459	`if (ierr.ne.0) goto 905`
460	`call clscdf(cstid,ierr)`	460	`call clscdf(cstid,ierr)`
461	`if (ierr.ne.0) goto 905`	461	`if (ierr.ne.0) goto 905`
462	`varname = 'PS'`	462	`varname = 'PS'`
463	`call getdef(fid,varname,ndim,mdv,vardim,`	463	`call getdef(fid,varname,ndim,mdv,vardim,`
464	`> varmin,varmax,stag,ierr)`	464	`> varmin,varmax,stag,ierr)`
465	`vardim(3) = nz1`	465	`vardim(3) = nz1`
466	`if (ierr.ne.0) goto 900`	466	`if (ierr.ne.0) goto 900`
467		467
468	`else`	468	`else`
469	`varname = fieldname`	469	`varname = fieldname`
470	`call getdef(fid,varname,ndim,mdv,vardim,`	470	`call getdef(fid,varname,ndim,mdv,vardim,`
471	`> varmin,varmax,stag,ierr)`	471	`> varmin,varmax,stag,ierr)`
472	`if (ierr.ne.0) goto 900`	472	`if (ierr.ne.0) goto 900`
473	`stagz=stag(3)`	473	`stagz=stag(3)`
474	`endif`	474	`endif`
475		475
476	`c Get time information (set time to first one in the file)`	476	`c Get time information (set time to first one in the file)`
477	`call gettimes(fid,times,ntimes,ierr)`	477	`call gettimes(fid,times,ntimes,ierr)`
478	`if (ierr.ne.0) goto 902`	478	`if (ierr.ne.0) goto 902`
479	`isok=0`	479	`isok=0`
480	`do i=1,ntimes`	480	`do i=1,ntimes`
481	`if (abs(time-times(i)).lt.eps) then`	481	`if (abs(time-times(i)).lt.eps) then`
482	`isok = 1`	482	`isok = 1`
483	`rtime = times(i)`	483	`rtime = times(i)`
484	`elseif (timecheck.eq.'no') then`	484	`elseif (timecheck.eq.'no') then`
485	`isok = 1`	485	`isok = 1`
486	`rtime = times(1)`	486	`rtime = times(1)`
487	`endif`	487	`endif`
488	`enddo`	488	`enddo`
489	`if ( isok.eq.0 ) goto 904`	489	`if ( isok.eq.0 ) goto 904`
490		490
491	`c Read field`	491	`c Read field`
492	`if ( ( fieldname.eq.'P' ).or.(fieldname.eq.'PLAY') ) then ! P, PLAY`	492	`if ( ( fieldname.eq.'P' ).or.(fieldname.eq.'PLAY') ) then ! P, PLAY`
493	`stagz = -0.5`	493	`stagz = -0.5`
494	`varname = 'PS'`	494	`varname = 'PS'`
495	`call getdat(fid,varname,rtime,0,ps,ierr)`	495	`call getdat(fid,varname,rtime,0,ps,ierr)`
496	`if (ierr.ne.0) goto 903`	496	`if (ierr.ne.0) goto 903`
497	`do i=1,nx`	497	`do i=1,nx`
498	`do j=1,ny`	498	`do j=1,ny`
499	`do k=1,nz`	499	`do k=1,nz`
500	`field(i,j,k)=aklay(k)+bklay(k)*ps(i,j)`	500	`field(i,j,k)=aklay(k)+bklay(k)*ps(i,j)`

Subversion Repositories lagranto.ecmwf

(root)/trunk/lib/ioinp.f – Rev 44 → 46