WebSVN – lagranto.um – Diff – /tags/2.1/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. 
 
      use netcdf
 
      implicit none
 
c     Declaration of subroutine parameters
      integer      fid              ! File identifier
      character*80 filename         ! Filename
 
c     Declaration of auxiliary variables
      integer      ierr
 
c     Open netcdf file
      ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
 
c     Exception handling
      return
      
 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).
 
      use netcdf
 
      implicit none
 
c     Declaration of subroutine parameters 
      integer      fid                 ! File identifier
      real         xmin,xmax,ymin,ymax ! Domain size
      real         dx,dy               ! Horizontal resolution
      integer      nx,ny,nz            ! Grid dimensions
      real         pollon,pollat       ! Longitude and latitude of pole
      real         p3(nx,ny,nz)        ! Staggered levels
      real         ps(nx,ny)           ! Surface pressure
      real         time                ! Time of the grid information
      real         ak(nz),bk(nz)       ! Ak and Bk for layers or levels
      real         stagz               ! Vertical staggering (0 or -0.5)
      character*80 fieldname           ! Variable from which to take grid info
      character*80 timecheck           ! Either 'yes' or 'no'
      
c     Numerical and physical parameters
      real          eps                 ! Numerical epsilon
      parameter    (eps=0.001)
 
c     Netcdf variables
      integer      vardim(4)
      real         varmin(4),varmax(4)
      real         mdv
      real         stag(4)
      integer      ndim
      character*80 cstfile
      integer      cstid
      real         times
      integer      ntimes
      real         aklay(nz),bklay(nz),aklev(nz),bklev(nz)
      integer      nvars
      character*80 vars(100)
      integer        dimids (nf90_max_var_dims)
      character*80   dimname(nf90_max_var_dims)
 
c     Auxiliary varaibles
      integer      ierr       
      integer      i,j,k
      integer      isok
      real         tmp(200)
      character*80 varname
      real         rtime
      integer      varid
      integer      cdfid
      real         rmax,rmin
 
c     Set file identifier
      if (fid.lt.0) then
        cdfid = -fid
      else 
        cdfid = fid
      endif
      
c     Get varid
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
         
c     Get numebr of dimensions -> ndim      
      ierr = nf90_inquire_variable(cdfid, varid, ndims  = ndim)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
         
c     Get dim IDs         
      ierr = nf90_inquire_variable(cdfid, varid, 
     >                                   dimids = dimids(1:ndim))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 
c     Get dimensions -> vardim(1:ndim)         
      do i=1,ndim
           ierr = nf90_inquire_dimension(cdfid, dimids(i), 
     >                               name = dimname(i) )
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
           ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      enddo
         
c     Get domain min -> varmin(1:3)        
      ierr  = nf90_get_att(cdfid, varid, "xmin", varmin(1))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr  = nf90_get_att(cdfid, varid, "ymin", varmin(2))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr  = nf90_get_att(cdfid, varid, "zmin", varmin(3))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
c     Get domain max -> varmax(1:3)     
      ierr  = nf90_get_att(cdfid, varid, "xmax", varmax(1))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr  = nf90_get_att(cdfid, varid, "ymax", varmax(2))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr  = nf90_get_att(cdfid, varid, "zmax", varmax(3))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
c     Get vertical staggering -> stagz           
      ierr  = nf90_get_att(cdfid, varid, "zstag", stagz)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
c     Get missing data value -> mdv         
      ierr  = nf90_get_att(cdfid, varid, "missing_value", mdv)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 
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 name of constants file -> cstfile
      ierr  = nf90_get_att(cdfid,nf90_global,
     >                    "constants_file_name ", cstfile )
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 
c     Get pole position from constants file
      ierr = NF90_OPEN(TRIM(cstfile),nf90_nowrite, cstid)
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
      varname='pollon'
      ierr = NF90_INQ_VARID(cstid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cstid,varid,pollon)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname='pollat'
      ierr = NF90_INQ_VARID(cstid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cstid,varid,pollat)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = NF90_CLOSE(cstid)
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
c     Skip the rest if fid < 0: 
      if ( fid.lt.0 ) goto 900      
      
c     Get ak and bk from constants file
      ierr = NF90_OPEN(TRIM(cstfile),nf90_nowrite, cstid)
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
      varname='aklev'
      ierr = NF90_INQ_VARID(cstid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cstid,varid,aklev)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname='bklev'
      ierr = NF90_INQ_VARID(cstid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cstid,varid,bklev)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname='aklay'
      ierr = NF90_INQ_VARID(cstid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cstid,varid,aklay)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname='bklay'
      ierr = NF90_INQ_VARID(cstid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cstid,varid,bklay)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = NF90_CLOSE(cstid)
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
c     print*,'AKLEV : '
c     print*, (aklev(i),i=1,20)
c     print*,'BKLEV : '
c     print*,(bklev(i),i=1,20)
c     print*,'AKLAY : '
c     print*,(aklay(i),i=1,20)
c     print*,'BKLAY : '
c     print*,(bklay(i),i=1,20)
      
c     Get time information (check if time is correct)
      varname = 'time'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,times)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      isok=0
      if (abs(time-times).lt.eps) then
               isok = 1
               rtime = times
      elseif (timecheck.eq.'no') then
               isok = 1
               rtime = times
      endif
      if ( isok.eq.0 ) then
         print*,' ERROR: time ',rtime,' not found on netCDF file' 
         stop
      endif
      
c     print*,'TIME : ',rtime
 
c     Read surface pressure
      varname='PS'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,ps)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 
c     Check that PS is in hPa - quick and dirty
      rmax = -1.e9
      rmin = +1.e9
      do i=1,nx
       do j=1,ny
        if ( (ps(i,j).gt.rmax).and.abs(ps(i,j)-mdv).gt.eps)  then
          rmax = ps(i,j)
        endif
        if ( (ps(i,j).lt.rmin).and.abs(ps(i,j)-mdv).gt.eps)  then
          rmin = ps(i,j)
        endif
       enddo
      enddo
      if ( (rmin.lt.200.).or.(rmax.gt.1500.) ) then
        print*,' ERROR: PS must be in hPa... Stop'
        stop
      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
      
c     Exit point
 900  continue    
   
      return
      
      end
 
c     ------------------------------------------------------------
c     Read wind information
c     ------------------------------------------------------------
 
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
     >                       timecheck)
 
c     Read the wind component <fieldname> from the file with identifier
c     <fid> and save it in the 3d array <field>. The vertical staggering 
c     information is provided in <stagz> and gives the reference to either
c     the layer or level field from <input_grid>. A consistency check is
c     performed to have an agreement with the grid specified by <xmin,xmax,
c     ymin,ymax,dx,dy,nx,ny,nz>.
 
      use netcdf
 
      implicit none
 
c     Declaration of variables and parameters
      integer      fid                 ! File identifier
      character*80 fieldname           ! Name of the wind field
      integer      nx,ny,nz            ! Dimension of fields
      real         field(nx,ny,nz)     ! 3d wind field
      real         stagz               ! Staggering in the z direction
      real         mdv                 ! Missing data flag
      real         xmin,xmax,ymin,ymax ! Domain size
      real         dx,dy               ! Horizontal resolution
      real         time                ! Time
      character*80 timecheck           ! Either 'yes' or 'no'
 
c     Numerical and physical parameters
      real        eps                 ! Numerical epsilon
      parameter  (eps=0.001)
      real        notimecheck         ! 'Flag' for no time check
      parameter  (notimecheck=7.26537)
 
c     Netcdf variables
      integer      ierr
      character*80 varname
      integer      vardim(4)
      real         varmin(4),varmax(4)
      real         stag(4)
      integer      ndim
      real         times(10)
      integer      ntimes
      character*80 cstfile
      integer      cstid
      real         aklay(200),bklay(200),aklev(200),bklev(200)
      real         ps(nx,ny)
      integer      dimids (nf90_max_var_dims)
      character*80 dimname(nf90_max_var_dims)
      integer      varid
      integer      cdfid
 
c     Auxiliary variables
      integer      isok
      integer      i,j,k
      integer      nz1
      real         rtime
 
c     Get varid
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
         
c     Get number of vertical levels -> vardim(3)    
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_inquire_variable(fid, varid, 
     >                                   dimids = dimids(1:ndim))
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      do i=1,ndim
           ierr = nf90_inquire_dimension(fid, dimids(i), 
     >                               name = dimname(i) )
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      enddo
 
c     Read data 
      varname=fieldname
      ierr = nf90_get_var(fid,varid,field)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
  
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
      if ( vardim(3).eq.1 ) then
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  field(i,j,k) = field(i,j,1)
               enddo
            enddo
         enddo
      endif
         
c     Exit point
      return
 
      end
 
c     ------------------------------------------------------------
c     Close input file
c     ------------------------------------------------------------
 
      subroutine input_close(fid)
 
c     Close the input file with file identifier <fid>.
 
      use netcdf
 
      implicit none
 
c     Declaration of subroutine parameters
      integer fid
 
c     Auxiliary variables
      integer ierr
 
c     Close file
      ierr = NF90_CLOSE(fid)
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
      end

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

Subversion Repositories lagranto.um

(root)/tags/2.1/lib/ioinp.f – Rev 11 → 12