WebSVN – lagranto.um – Diff – /tags/1.0/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
 
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
         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        We want the longitudes within -180 ... + 180
         if ( xmin.lt.-180.) then
            xmin = xmin + 360.
            xmax = xmax + 360.
         else if ( xmax.gt.360. ) then
            xmin = xmin - 360.
            xmax = xmax - 360.
         endif
 
c        Get pole position and number of vertical levels
         call getcfn(-fid,cstfile,ierr)
         if (ierr.ne.0) goto 903
         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
 
c     Get non-constant grid parameters (surface pressure and vertical grid)
      else
         
c        Get full grid info - in particular staggering flag
         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        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        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        Set the grid dimensions and constants
         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        We want the longitudes within -180 ... + 180
         if ( xmin.lt.-180.) then
            xmin = xmin + 360.
            xmax = xmax + 360.
         else if ( xmax.gt.360. ) then
            xmin = xmin - 360.
            xmax = xmax - 360.
         endif
 
c        Read surface pressure
         varname='PS'
         call getdat(fid,varname,rtime,0,ps,ierr)
         if (ierr.ne.0) goto 904
 
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     Exception handling
      return
      
 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

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

Subversion Repositories lagranto.um

(root)/tags/1.0/lib/ioinp.f – Rev 3 → 4