Subversion Repositories lagranto.ecmwf

Rev

Rev 23 | Rev 27 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 23 Rev 25
1 
c     ************************************************************
 1 
c     ************************************************************
2 
c     * This package provides input routines to read the wind    *
 2 
c     * This package provides input routines to read the wind    *
3 
c     * and other fields from IVE necdf files. The routines are  *
 3 
c     * and other fields from IVE necdf files. The routines are  *
4 
c     *                                                          *
 4 
c     *                                                          *
5 
c     * 1) input_open  : to open a data file                     *
 5 
c     * 1) input_open  : to open a data file                     *
6 
c     * 2) input_grid  : to read the grid information, including *
 6 
c     * 2) input_grid  : to read the grid information, including *
7 
c     *                  the vertical levels                     *
 7 
c     *                  the vertical levels                     *
8 
c     * 3) input_wind  : to read the wind components             *
 8 
c     * 3) input_wind  : to read the wind components             *
9 
c     * 4) input_close : to close an input file                  *
 9 
c     * 4) input_close : to close an input file                  *
10 
c     *                                                          *
 10 
c     *                                                          *
11 
c     * The file is characterised by an filename <filename> and  *
 11 
c     * The file is characterised by an filename <filename> and  *
12 
c     * a file identifier <fid>. The horizontal grid is given by *
 12 
c     * a file identifier <fid>. The horizontal grid is given by *
13 
c     * <xmin,xmax,ymin,ymax,dx,dy,nx,ny> where the pole of the  *
 13 
c     * <xmin,xmax,ymin,ymax,dx,dy,nx,ny> where the pole of the  *
14 
c     * rotated grid is given by <pollon,pollat>. The vertical   *
 14 
c     * rotated grid is given by <pollon,pollat>. The vertical   *
15 
c     * grid is characterised by the surface pressure <ps> and   *
 15 
c     * grid is characterised by the surface pressure <ps> and   *
16 
c     * the pressure at staggered <slev> and unstaggered <ulev>  *
 16 
c     * the pressure at staggered <slev> and unstaggered <ulev>  *
17 
c     * levels. The number of levels is given by <nz>. Finally,  *
 17 
c     * levels. The number of levels is given by <nz>. Finally,  *
18 
c     * the retrieval of the wind <field> with name <fieldname>  *
 18 
c     * the retrieval of the wind <field> with name <fieldname>  *
19 
c     * is characterised by a <time> and a missing data value    *
 19 
c     * is characterised by a <time> and a missing data value    *
20 
c     * <mdv>.                                                   *
 20 
c     * <mdv>.                                                   *
21 
c     *                                                          *
 21 
c     *                                                          *
22 
c     * Author: Michael Sprenger, Autumn 2008                    *
 22 
c     * Author: Michael Sprenger, Autumn 2008                    *
23 
c     ************************************************************
 23 
c     ************************************************************
24 
 
 24 
 
25 
c     ------------------------------------------------------------
 25 
c     ------------------------------------------------------------
26 
c     Open input file
 26 
c     Open input file
27 
c     ------------------------------------------------------------
 27 
c     ------------------------------------------------------------
28 
 
 28 
 
29 
      subroutine input_open (fid,filename)
 29 
      subroutine input_open (fid,filename)
30 
 
 30 
 
31 
c     Open the input file with filename <filename> and return the
 31 
c     Open the input file with filename <filename> and return the
32 
c     file identifier <fid> for further reference.
32 
c     file identifier <fid> for further reference.
33 
 
 33 
 
34 
      use netcdf
 34 
      use netcdf
35 
 
 35 
 
36 
      implicit none
 36 
      implicit none
37 
 
 37 
 
38 
c     Declaration of subroutine parameters
 38 
c     Declaration of subroutine parameters
39 
      integer      fid              ! File identifier
 39 
      integer      fid              ! File identifier
40 
      character*80 filename         ! Filename
 40 
      character*80 filename         ! Filename
41 
 
 41 
 
42 
c     Declaration of auxiliary variables
 42 
c     Declaration of auxiliary variables
43 
      integer      ierr
 43 
      integer      ierr
44 
 
 44 
 
45 
c     Open netcdf file
 45 
c     Open netcdf file
46 
      ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
 46 
      ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
47 
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
 47 
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
48 
 
 48 
 
49 
c     Exception handling
 49 
c     Exception handling
50 
      return
 50 
      return
51 
 
 51 
 
52 
      end
 52 
      end
53
53
54 
c     ------------------------------------------------------------
 54 
c     ------------------------------------------------------------
55 
c     Read information about the grid
 55 
c     Read information about the grid
56 
c     ------------------------------------------------------------
 56 
c     ------------------------------------------------------------
57
57
58 
      subroutine input_grid
58 
      subroutine input_grid
59 
     >                   (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
 59 
     >                   (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
60 
     >                    time,pollon,pollat,p3,ps,nz,ak,bk,stagz,
 60 
     >                    time,pollon,pollat,p3,ps,nz,ak,bk,stagz,
61 
     >                    timecheck)
 61 
     >                    timecheck)
62 
 
 62 
 
63 
c     Read grid information at <time> from file with identifier <fid>.
63 
c     Read grid information at <time> from file with identifier <fid>.
64 
c     The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>
 64 
c     The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>
65 
c     with pole position at <pollon,pollat> and grid dimension <nx,ny>.
 65 
c     with pole position at <pollon,pollat> and grid dimension <nx,ny>.
66 
c     The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either
 66 
c     The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either
67 
c     on the staggered or unstaggered grid (with <stagz> as the flag).
 67 
c     on the staggered or unstaggered grid (with <stagz> as the flag).
68 
c     The surface pressure is given in <ps(nx,ny)>. If <fid> is negative,
68 
c     The surface pressure is given in <ps(nx,ny)>. If <fid> is negative,
69 
c     only the grid dimensions and grid parameters (xmin...pollat,nz) are
69 
c     only the grid dimensions and grid parameters (xmin...pollat,nz) are
70 
c     determined and returned (this is needed for dynamical allocation of
70 
c     determined and returned (this is needed for dynamical allocation of
71 
c     memory).
 71 
c     memory).
72 
 
 72 
 
73 
      use netcdf
 73 
      use netcdf
74 
 
 74 
 
75 
      implicit none
 75 
      implicit none
76 
 
 76 
 
77 
c     Declaration of subroutine parameters
77 
c     Declaration of subroutine parameters
78 
      integer      fid                 ! File identifier
 78 
      integer      fid                 ! File identifier
79 
      real         xmin,xmax,ymin,ymax ! Domain size
 79 
      real         xmin,xmax,ymin,ymax ! Domain size
80 
      real         dx,dy               ! Horizontal resolution
 80 
      real         dx,dy               ! Horizontal resolution
81 
      integer      nx,ny,nz            ! Grid dimensions
 81 
      integer      nx,ny,nz            ! Grid dimensions
82 
      real         pollon,pollat       ! Longitude and latitude of pole
 82 
      real         pollon,pollat       ! Longitude and latitude of pole
83 
      real         p3(nx,ny,nz)        ! Staggered levels
 83 
      real         p3(nx,ny,nz)        ! Staggered levels
84 
      real         ps(nx,ny)           ! Surface pressure
 84 
      real         ps(nx,ny)           ! Surface pressure
85 
      real         time                ! Time of the grid information
 85 
      real         time                ! Time of the grid information
86 
      real         ak(nz),bk(nz)       ! Ak and Bk for layers or levels
 86 
      real         ak(nz),bk(nz)       ! Ak and Bk for layers or levels
87 
      real         stagz               ! Vertical staggering (0 or -0.5)
 87 
      real         stagz               ! Vertical staggering (0 or -0.5)
88 
      character*80 fieldname           ! Variable from which to take grid info
 88 
      character*80 fieldname           ! Variable from which to take grid info
89 
      character*80 timecheck           ! Either 'yes' or 'no'
 89 
      character*80 timecheck           ! Either 'yes' or 'no'
90
90
91 
c     Numerical and physical parameters
 91 
c     Numerical and physical parameters
92 
      real          eps                 ! Numerical epsilon
 92 
      real          eps                 ! Numerical epsilon
93 
      parameter    (eps=0.001)
 93 
      parameter    (eps=0.001)
94 
 
 94 
 
95 
c     Netcdf variables
 95 
c     Netcdf variables
96 
      integer      vardim(4)
 96 
      integer      vardim(4)
97 
      real         varmin(4),varmax(4)
 97 
      real         varmin(4),varmax(4)
98 
      real         mdv
 98 
      real         mdv
99 
      real         stag(4)
 99 
      real         stag(4)
100 
      integer      ndim
 100 
      integer      ndim
101 
      character*80 cstfile
 101 
      character*80 cstfile
102 
      integer      cstid
 102 
      integer      cstid
103 
      integer      nvars
 103 
      integer      nvars
104 
      character*80 vars(100)
 104 
      character*80 vars(100)
105 
      integer        dimids (nf90_max_var_dims)
 105 
      integer        dimids (nf90_max_var_dims),dimid
106 
      character*80   dimname(nf90_max_var_dims)
 106 
      character*80   dimname(nf90_max_var_dims)
- 
 
 107 
      character*80   stdname
107 
      real,allocatable, dimension (:)     :: lon,lat,lev
 108 
      real,allocatable, dimension (:)     :: lon,lat,lev
108 
      real,allocatable, dimension (:)     :: times
 109 
      real,allocatable, dimension (:)     :: times
109 
      real,allocatable, dimension (:,:)   :: tmp2
 110 
      real,allocatable, dimension (:,:)   :: tmp2
110 
      real,allocatable, dimension (:,:,:) :: tmp3
 111 
      real,allocatable, dimension (:,:,:) :: tmp3
111 
      real,allocatable, dimension (:)     :: aktmp,bktmp
 112 
      real,allocatable, dimension (:)     :: aktmp,bktmp
112 
      character*80  units
 113 
      character*80  units
- 
 
 114 
      character*80  leveltype
- 
 
 115 
      integer       nakbktmp
- 
 
 116 
      integer       vertical_swap
113 
 
 117 
 
114 
c     Auxiliary variables
 118 
c     Auxiliary variables
115 
      integer      ierr
119 
      integer      ierr
116 
      integer      i,j,k
 120 
      integer      i,j,k
117 
      integer      isok
 121 
      integer      isok
118 
      real         tmp(200)
 122 
      real         tmp(200)
119 
      character*80 varname
 123 
      character*80 varname
120 
      real         rtime
 124 
      real         rtime
121 
      integer      varid
 125 
      integer      varid
122 
      integer      cdfid
 126 
      integer      cdfid
123 
      integer      stat
 127 
      integer      stat
124 
      real         delta
 128 
      real         delta
125 
      integer      closear
 129 
      integer      closear
126 
      real         maxps,minps
 130 
      real         maxps,minps
127 
 
 131 
 
128 
c     ------ Set file identifier --------------------------------------
 132 
c     ---- Read data from netCDF file as they are ---------------------
- 
 
 133 
 
- 
 
 134 
c     Set file identifier
129 
      if (fid.lt.0) then
 135 
      if (fid.lt.0) then
130 
        cdfid = -fid
 136 
        cdfid = -fid
131 
      else
137 
      else
132 
        cdfid = fid
 138 
        cdfid = fid
133 
      endif
 139 
      endif
134 
 
 140 
 
135 
c     Special handling if 3D pressure is
 141 
c     Special handling if 3D pressure is
136 
      if ( fieldname.eq.'P' ) then
 142 
      if ( fieldname.eq.'P' ) then
137 
         fieldname = 'U'
 143 
         fieldname = 'U'
138 
      endif
 144 
      endif
139 
 
 145 
 
140 
c     Get number of dimensions of variable -> ndim
 146 
c     Get number of dimensions of variable -> ndim
141 
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
 147 
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
142 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 148 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
143 
      ierr = nf90_inquire_variable(cdfid, varid, ndims  = ndim)
 149 
      ierr = nf90_inquire_variable(cdfid, varid, ndims  = ndim)
144 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 150 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
145 
      if ( ndim.ne.4 ) then
 151 
      if ( ndim.ne.4 ) then
146 
          print*,' ERROR: netCDF variables need to be 4D'
 152 
          print*,' ERROR: netCDF variables need to be 4D'
147 
          print*,'      ',trim(fieldname)
 153 
          print*,'      ',trim(fieldname)
148 
          stop
 154 
          stop
149 
      endif
 155 
      endif
150 
 
 156 
 
151 
c     Get dimensions -> vardim(1:ndim),dimname(1:ndim)
 157 
c     Get dimensions -> vardim(1:ndim),dimname(1:ndim)
152 
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
 158 
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
153 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 159 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
154 
      ierr = nf90_inquire_variable(cdfid, varid,
160 
      ierr = nf90_inquire_variable(cdfid, varid,
155 
     >                                   dimids = dimids(1:ndim))
 161 
     >                                   dimids = dimids(1:ndim))
156 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 162 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
157 
      do i=1,ndim
 163 
      do i=1,ndim
158 
           ierr = nf90_inquire_dimension(cdfid, dimids(i),
164 
           ierr = nf90_inquire_dimension(cdfid, dimids(i),
159 
     >                               name = dimname(i) )
 165 
     >                               name = dimname(i) )
160 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 166 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
161 
           ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))
 167 
           ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))
162 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 168 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
163 
      enddo
 169 
      enddo
164 
 
 170 
 
- 
 
 171 
c     Get dimension of AK,BK
- 
 
 172 
      varname = 'nhym'
- 
 
 173 
      ierr = NF90_INQ_DIMID(cdfid,varname,dimid)
- 
 
 174 
      ierr = nf90_inquire_dimension(cdfid, dimid,len=nakbktmp)
- 
 
 175 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 176 
 
165 
c     Check whether the list of dimensions is OK
 177 
c     Check whether the list of dimensions is OK
166 
      if ( ( dimname(1).ne.'lon'  ).or.
 178 
      if ( ( dimname(1).ne.'lon'  ).or.
167 
     >     ( dimname(2).ne.'lat'  ).or.
179 
     >     ( dimname(2).ne.'lat'  ).or.
168 
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
 180 
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
169 
     >     ( dimname(4).ne.'time' ) )
 181 
     >     ( dimname(4).ne.'time' ) )
170 
     >then
 182 
     >then
171 
        print*,' ERROR: the dimensions of the variable are not correct'
 183 
        print*,' ERROR: the dimensions of the variable are not correct'
172 
        print*,'        expected -> lon / lat / lev / time'
 184 
        print*,'        expected -> lon / lat / lev / time'
173 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
 185 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
174 
        stop
 186 
        stop
175 
      endif
 187 
      endif
176 
 
 188 
 
- 
 
 189 
c     Check about the type of vertical levels
- 
 
 190 
      varname=dimname(3)
- 
 
 191 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 192 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 193 
      ierr = nf90_get_att(cdfid, varid, "standard_name", leveltype)
- 
 
 194 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 195 
      if ( (leveltype.ne.'hybrid_sigma_pressure').and.
- 
 
 196 
     >     (leveltype.ne.'air_pressure'         ) )
- 
 
 197 
     >then
- 
 
 198 
         print*,' ERROR: input netCDF data must be on hybrid-sigma'
- 
 
 199 
         print*,'        or air pressure levels!',trim(leveltype)
- 
 
 200 
         stop
- 
 
 201 
      endif
- 
 
 202 
 
- 
 
 203 
c     Check that vardim(3)==#AK,BK for hybrid-sigmal levels
- 
 
 204 
      if ( (leveltype.eq.'hybrid_sigma_pressure').and.
- 
 
 205 
     >     (dimname(3).eq.'lev') )
- 
 
 206 
     >then
- 
 
 207 
        if ( nakbktmp.ne.vardim(3) ) then
- 
 
 208 
           print*,' ERROR: for hybrid-sigma pressure levels, #AK,BK'
- 
 
 209 
           print*,'        must agree with number of vertical levels'
- 
 
 210 
           print*,'        ',vardim(3),nakbktmp
- 
 
 211 
           stop
- 
 
 212 
        endif
- 
 
 213 
      endif
- 
 
 214 
 
177 
c     Allocate memory for reading arrays
 215 
c     Allocate memory for reading arrays
178 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
 216 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
179 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
 217 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
180 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
 218 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
181 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
 219 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
182 
      allocate(lon(vardim(1)),stat=stat)
 220 
      allocate(lon(vardim(1)),stat=stat)
183 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
221 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
184 
      allocate(lat(vardim(2)),stat=stat)
 222 
      allocate(lat(vardim(2)),stat=stat)
185 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
223 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
- 
 
 224 
      allocate(lev(vardim(3)),stat=stat)
- 
 
 225 
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
186 
      allocate(times(vardim(4)),stat=stat)
 226 
      allocate(times(vardim(4)),stat=stat)
187 
      if (stat.ne.0) print*,'*** error allocating array times   ***'
 227 
      if (stat.ne.0) print*,'*** error allocating array times   ***'
188 
      allocate(aktmp(vardim(3)),stat=stat)
 228 
      allocate(aktmp(nakbktmp),stat=stat)
189 
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
 229 
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
190 
      allocate(bktmp(vardim(3)),stat=stat)
 230 
      allocate(bktmp(nakbktmp),stat=stat)
191 
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
 231 
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
192 
 
 232 
 
193 
c     Get domain longitudes and latitudes
 233 
c     Get domain longitudes, latitudes and levels
194 
      varname = dimname(1)
 234 
      varname = dimname(1)
195 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 235 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
196 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 236 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
197 
      ierr = nf90_get_var(cdfid,varid,lon)
 237 
      ierr = nf90_get_var(cdfid,varid,lon)
198 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 238 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
199 
      varname = dimname(2)
 239 
      varname = dimname(2)
200 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 240 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
201 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 241 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
202 
      ierr = nf90_get_var(cdfid,varid,lat)
 242 
      ierr = nf90_get_var(cdfid,varid,lat)
203 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 243 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 244 
      varname = dimname(3)
- 
 
 245 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 246 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 247 
      ierr = nf90_get_var(cdfid,varid,lev)
- 
 
 248 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
204
249
205 
c     Get ak and bk
 250 
c     Get ak and bk
206 
      varname='hyam'
 251 
      varname='hyam'
207 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 252 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
208 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 253 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
209 
      ierr = nf90_get_var(cdfid,varid,aktmp)
 254 
      ierr = nf90_get_var(cdfid,varid,aktmp)
210 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 255 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
211 
      varname='hybm'
 256 
      varname='hybm'
212 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 257 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
213 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 258 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
214 
      ierr = nf90_get_var(cdfid,varid,bktmp)
 259 
      ierr = nf90_get_var(cdfid,varid,bktmp)
215 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 260 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
216 
 
 261 
 
217 
c     Check that unit of ak is in hPa - if necessary correct it
 262 
c     Check that unit of ak is in hPa - if necessary correct it
218 
      varname='hyam'
 263 
      varname='hyam'
219 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 264 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
220 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 265 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
221 
      ierr = nf90_get_att(cdfid, varid, "units", units)
 266 
      ierr = nf90_get_att(cdfid, varid, "units", units)
222 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 267 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
223 
      if ( units.eq.'Pa' ) then
 268 
      if ( units.eq.'Pa' ) then
224 
         do k=1,vardim(3)
 269 
         do k=1,nakbktmp
225 
            aktmp(k) = 0.01 * aktmp(k)
 270 
            aktmp(k) = 0.01 * aktmp(k)
226 
         enddo
 271 
         enddo
227 
      endif
 272 
      endif
228 
 
 273 
 
- 
 
 274 
c     Check that unit of lev is in hPa - if necessary correct it
- 
 
 275 
      if ( leveltype.eq.'air_pressure' ) then
- 
 
 276 
         varname='lev'
- 
 
 277 
         ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 278 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 279 
         ierr = nf90_get_att(cdfid, varid, "units", units)
- 
 
 280 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 281 
         if ( units.eq.'Pa' ) then
- 
 
 282 
            do k=1,vardim(3)
- 
 
 283 
               lev(k) = 0.01 * lev(k)
- 
 
 284 
            enddo
- 
 
 285 
         endif
- 
 
 286 
      endif
- 
 
 287 
 
- 
 
 288 
c     Decide whether to swap vertical levels
- 
 
 289 
      vertical_swap = 1
- 
 
 290 
      if ( leveltype.eq.'hybrid_sigma_pressure') then
- 
 
 291 
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
- 
 
 292 
     >       (aktmp(2) + bktmp(2) * 1000.) )
- 
 
 293 
     >  then
- 
 
 294 
          vertical_swap = 0
- 
 
 295 
        endif
- 
 
 296 
      elseif ( leveltype.eq.'air_pressure') then
- 
 
 297 
        if ( lev(1).gt.lev(2) ) then
- 
 
 298 
          vertical_swap = 0
- 
 
 299 
        endif
- 
 
 300 
      endif
- 
 
 301 
c      print*,' Vertical SWAP P -> ', vertical_swap
- 
 
 302 
 
229 
c     Get time information (check if time is correct)
 303 
c     Get time information (check if time is correct)
230 
      varname = 'time'
 304 
      varname = 'time'
231 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 305 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
232 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 306 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
233 
      ierr = nf90_get_var(cdfid,varid,times)
 307 
      ierr = nf90_get_var(cdfid,varid,times)
234 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 308 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
235 
      isok=0
 309 
      isok=0
236 
      do i=1,vardim(4)
 310 
      do i=1,vardim(4)
237 
        if (abs(time-times(i)).lt.eps) then
 311 
        if (abs(time-times(i)).lt.eps) then
238 
               isok = 1
 312 
               isok = 1
239 
               rtime = times(i)
 313 
               rtime = times(i)
240 
        elseif (timecheck.eq.'no') then
 314 
        elseif (timecheck.eq.'no') then
241 
               isok = 1
 315 
               isok = 1
242 
               rtime = times(1)
 316 
               rtime = times(1)
243 
        endif
 317 
        endif
244 
      enddo
 318 
      enddo
245 
      if ( isok.eq.0 ) then
 319 
      if ( isok.eq.0 ) then
246 
         print*,' ERROR: time ',rtime,' not found on netCDF file'
320 
         print*,' ERROR: time ',rtime,' not found on netCDF file'
247 
         stop
 321 
         stop
248 
      endif
 322 
      endif
249 
 
 323 
 
250 
c     Read surface pressure
 324 
c     Read surface pressure
251 
      varname='PS'
 325 
      varname='PS'
252 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 326 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
253 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 327 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
254 
      ierr = nf90_get_var(cdfid,varid,tmp2)
 328 
      ierr = nf90_get_var(cdfid,varid,tmp2)
255 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 329 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
256
330
257 
c     Check that surface pressure is in hPa
 331 
c     Check that surface pressure is in hPa
258 
      maxps = -1.e39
 332 
      maxps = -1.e19
259 
      minps =  1.e39
 333 
      minps =  1.e19
260 
      do i=1,vardim(1)
 334 
      do i=1,vardim(1)
261 
        do j=1,vardim(2)
 335 
        do j=1,vardim(2)
262 
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
 336 
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
263 
             if (tmp2(i,j).lt.minps) minps = tmp2(i,j)
 337 
             if (tmp2(i,j).lt.minps) minps = tmp2(i,j)
264 
        enddo
 338 
        enddo
265 
      enddo
 339 
      enddo
266 
      if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then
 340 
      if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then
267 
         print*,' ERROR: surface pressre PS must be in hPa'
 341 
         print*,' ERROR: surface pressre PS must be in hPa'
268 
         print*,'       ',maxps,minps
 342 
         print*,'       ',maxps,minps
269 
         stop
 343 
         stop
270 
      endif
 344 
      endif
271 
 
 345 
 
272 
c     Calculate layer and level pressures
 - 
 
273 
      do i=1,vardim(1)
 - 
 
274 
         do j=1,vardim(2)
 - 
 
275 
               do k=1,vardim(3)
 - 
 
276 
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
 346 
c     ---- Define output of subroutine --------------------------------
277 
               enddo
 - 
 
278 
         enddo
 - 
 
279 
      enddo
 - 
 
280 
 
 347 
 
281 
c     Set the grid dimensions and constants
 348 
c     Set the grid dimensions and constants
282 
      nx      = vardim(1)
 349 
      nx      = vardim(1)
283 
      ny      = vardim(2)
 350 
      ny      = vardim(2)
284 
      nz      = vardim(3)
 351 
      nz      = vardim(3)
285 
      xmin    = lon(1)
 352 
      xmin    = lon(1)
286 
      ymin    = lat(1)
 353 
      ymin    = lat(1)
287 
      xmax    = lon(nx)
 354 
      xmax    = lon(nx)
288 
      ymax    = lat(ny)
 355 
      ymax    = lat(ny)
289 
      dx      = (xmax-xmin)/real(nx-1)
 356 
      dx      = (xmax-xmin)/real(nx-1)
290 
      dy      = (ymax-ymin)/real(ny-1)
 357 
      dy      = (ymax-ymin)/real(ny-1)
291 
      pollon  = 0.
 358 
      pollon  = 0.
292 
      pollat  = 90.
 359 
      pollat  = 90.
293 
      stagz   = 0.
 360 
      stagz   = 0.
294 
      delta   = xmax-xmin-360.
 361 
      delta   = xmax-xmin-360.
295 
      if (abs(delta+dx).lt.eps) then
 362 
      if (abs(delta+dx).lt.eps) then
296 
          xmax    = xmax + dx
 363 
          xmax    = xmax + dx
297 
          nx      = nx + 1
 364 
          nx      = nx + 1
298 
          closear = 1
 365 
          closear = 1
299 
      else
 366 
      else
300 
          closear = 0
 367 
          closear = 0
301 
      endif
 368 
      endif
302 
 
 369 
 
303 
c     Save the output arrays (if fid>0) - close arrays on request
 370 
c     Save the output arrays (if fid>0) - close arrays on request
304 
      if ( fid.gt.0 ) then
 371 
      if ( fid.gt.0 ) then
305 
 
 372 
 
- 
 
 373 
c        Calculate layer pressures
- 
 
 374 
         if (leveltype.eq.'hybrid_sigma_pressure' ) then
- 
 
 375 
            do i=1,vardim(1)
- 
 
 376 
              do j=1,vardim(2)
- 
 
 377 
                 do k=1,vardim(3)
- 
 
 378 
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
- 
 
 379 
                 enddo
- 
 
 380 
              enddo
- 
 
 381 
           enddo
- 
 
 382 
         elseif (leveltype.eq.'air_pressure' ) then
- 
 
 383 
           do i=1,vardim(1)
- 
 
 384 
              do j=1,vardim(2)
- 
 
 385 
                 do k=1,vardim(3)
- 
 
 386 
                  tmp3(i,j,k)=lev(k)
- 
 
 387 
                 enddo
- 
 
 388 
              enddo
- 
 
 389 
           enddo
- 
 
 390 
         endif
- 
 
 391 
 
- 
 
 392 
c        Get PS - close array on demand
306 
         do j=1,vardim(2)
 393 
         do j=1,vardim(2)
307 
           do i=1,vardim(1)
 394 
           do i=1,vardim(1)
308 
             ps(i,j) = tmp2(i,j)
 395 
             ps(i,j) = tmp2(i,j)
309 
           enddo
 396 
           enddo
310 
           if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)
 397 
           if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)
311 
         enddo
 398 
         enddo
312 
 
 399 
 
- 
 
 400 
c        Get P3 - close array on demand + vertical swap
313 
         do j=1,vardim(2)
 401 
         do j=1,vardim(2)
314 
           do k=1,vardim(3)
 402 
           do k=1,vardim(3)
315 
             do i=1,vardim(1)
 403 
             do i=1,vardim(1)
- 
 
 404 
               if ( vertical_swap.eq.1 ) then
316 
               p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)
 405 
                  p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)
- 
 
 406 
               else
- 
 
 407 
                  p3(i,j,k) = tmp3(i,j,k)
- 
 
 408 
               endif
317 
             enddo
 409 
             enddo
318 
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
 410 
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
319 
           enddo
 411 
           enddo
320 
         enddo
 412 
         enddo
321 
 
 413 
 
- 
 
 414 
c        Get AK,BK - vertical swap on demand
- 
 
 415 
         if ( leveltype.eq.'hybrid_sigma_pressure' ) then
322 
         do k=1,vardim(3)
 416 
           do k=1,vardim(3)
- 
 
 417 
              if ( vertical_swap.eq.1 ) then
323 
            ak(k) = aktmp(vardim(3)-k+1)
 418 
                 ak(k) = aktmp(vardim(3)-k+1)
324 
            bk(k) = bktmp(vardim(3)-k+1)
 419 
                 bk(k) = bktmp(vardim(3)-k+1)
- 
 
 420 
              endif
- 
 
 421 
           enddo
- 
 
 422 
         elseif (leveltype.eq.'air_pressure' ) then
- 
 
 423 
           do k=1,vardim(3)
- 
 
 424 
              if ( vertical_swap.eq.1 ) then
- 
 
 425 
                 ak(k) = lev(vardim(3)-k+1)
- 
 
 426 
                 bk(k) = 0.
- 
 
 427 
              else
- 
 
 428 
                ak(k) = lev(k)
- 
 
 429 
                bk(k) = 0.
- 
 
 430 
              endif
325 
         enddo
 431 
           enddo
- 
 
 432 
         endif
326 
 
 433 
 
327 
      endif
 434 
      endif
328 
 
 435 
 
- 
 
 436 
 
329 
      return
 437 
      return
330
438
331 
      end
 439 
      end
332 
 
 440 
 
333 
c     ------------------------------------------------------------
 441 
c     ------------------------------------------------------------
334 
c     Read wind information
 442 
c     Read wind information
335 
c     ------------------------------------------------------------
 443 
c     ------------------------------------------------------------
336 
 
 444 
 
337 
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
 445 
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
338 
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
 446 
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
339 
     >                       timecheck)
 447 
     >                       timecheck)
340 
 
 448 
 
341 
c     Read the wind component <fieldname> from the file with identifier
 449 
c     Read the wind component <fieldname> from the file with identifier
342 
c     <fid> and save it in the 3d array <field>. The vertical staggering
450 
c     <fid> and save it in the 3d array <field>. The vertical staggering
343 
c     information is provided in <stagz> and gives the reference to either
 451 
c     information is provided in <stagz> and gives the reference to either
344 
c     the layer or level field from <input_grid>. A consistency check is
 452 
c     the layer or level field from <input_grid>. A consistency check is
345 
c     performed to have an agreement with the grid specified by <xmin,xmax,
 453 
c     performed to have an agreement with the grid specified by <xmin,xmax,
346 
c     ymin,ymax,dx,dy,nx,ny,nz>.
 454 
c     ymin,ymax,dx,dy,nx,ny,nz>.
347 
 
 455 
 
348 
      use netcdf
 456 
      use netcdf
349 
 
 457 
 
350 
      implicit none
 458 
      implicit none
351 
 
 459 
 
352 
c     Declaration of variables and parameters
 460 
c     Declaration of variables and parameters
353 
      integer      fid                 ! File identifier
 461 
      integer      fid                 ! File identifier
354 
      character*80 fieldname           ! Name of the wind field
 462 
      character*80 fieldname           ! Name of the wind field
355 
      integer      nx,ny,nz            ! Dimension of fields
 463 
      integer      nx,ny,nz            ! Dimension of fields
356 
      real         field(nx,ny,nz)     ! 3d wind field
 464 
      real         field(nx,ny,nz)     ! 3d wind field
357 
      real         stagz               ! Staggering in the z direction
 465 
      real         stagz               ! Staggering in the z direction
358 
      real         mdv                 ! Missing data flag
 466 
      real         mdv                 ! Missing data flag
359 
      real         xmin,xmax,ymin,ymax ! Domain size
 467 
      real         xmin,xmax,ymin,ymax ! Domain size
360 
      real         dx,dy               ! Horizontal resolution
 468 
      real         dx,dy               ! Horizontal resolution
361 
      real         time                ! Time
 469 
      real         time                ! Time
362 
      character*80 timecheck           ! Either 'yes' or 'no'
 470 
      character*80 timecheck           ! Either 'yes' or 'no'
363 
 
 471 
 
364 
c     Numerical and physical parameters
 472 
c     Numerical and physical parameters
365 
      real        eps                 ! Numerical epsilon
 473 
      real        eps                 ! Numerical epsilon
366 
      parameter  (eps=0.001)
 474 
      parameter  (eps=0.001)
367 
      real        notimecheck         ! 'Flag' for no time check
 475 
      real        notimecheck         ! 'Flag' for no time check
368 
      parameter  (notimecheck=7.26537)
 476 
      parameter  (notimecheck=7.26537)
369 
 
 477 
 
370 
c     Netcdf variables
 478 
c     Netcdf variables
371 
      integer      ierr
 479 
      integer      ierr
372 
      character*80 varname
 480 
      character*80 varname
373 
      integer      vardim(4)
 481 
      integer      vardim(4)
374 
      real         varmin(4),varmax(4)
 482 
      real         varmin(4),varmax(4)
375 
      real         stag(4)
 483 
      real         stag(4)
376 
      integer      ndim
 484 
      integer      ndim
377 
      real         times(10)
 485 
      real         times(10)
378 
      integer      ntimes
 486 
      integer      ntimes
379 
      character*80 cstfile
 487 
      character*80 cstfile
380 
      integer      cstid
 488 
      integer      cstid
381 
      real         aklay(200),bklay(200),aklev(200),bklev(200)
 489 
      real         aklay(200),bklay(200),aklev(200),bklev(200)
382 
      real         ps(nx,ny)
 490 
      real         ps(nx,ny)
383 
      integer      dimids (nf90_max_var_dims)
 491 
      integer      dimids (nf90_max_var_dims)
384 
      character*80 dimname(nf90_max_var_dims)
 492 
      character*80 dimname(nf90_max_var_dims)
385 
      integer      varid
 493 
      integer      varid
386 
      integer      cdfid
 494 
      integer      cdfid
387 
      real,allocatable, dimension (:)     :: lon,lat
 495 
      real,allocatable, dimension (:)     :: lon,lat,lev
388 
      real,allocatable, dimension (:,:)   :: tmp2
 496 
      real,allocatable, dimension (:,:)   :: tmp2
389 
      real,allocatable, dimension (:,:,:) :: tmp3
 497 
      real,allocatable, dimension (:,:,:) :: tmp3
- 
 
 498 
      real,allocatable, dimension (:)     :: aktmp,bktmp
- 
 
 499 
      character*80  leveltype
- 
 
 500 
      integer       vertical_swap
- 
 
 501 
      character*80  units
- 
 
 502 
      integer       nakbktmp
- 
 
 503 
      integer       dimid
390 
 
 504 
 
391 
c     Auxiliary variables
 505 
c     Auxiliary variables
392 
      integer      isok
 506 
      integer      isok
393 
      integer      i,j,k
 507 
      integer      i,j,k
394 
      integer      nz1
 508 
      integer      nz1
395 
      real         rtime
 509 
      real         rtime
396 
      integer      closear
 510 
      integer      closear
397 
      integer      stat
 511 
      integer      stat
398 
      real         delta
 512 
      real         delta
399 
 
 513 
 
400 
c     Get the number of dimensions -> ndim
 514 
c     Get the number of dimensions -> ndim
401 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
 515 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
402 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 516 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
403 
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
 517 
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
404 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 518 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
405 
 
 519 
 
406 
c     Get the dimensions of the arrays -> varid(1:ndim)
 520 
c     Get the dimensions of the arrays -> varid(1:ndim)
407 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
 521 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
408 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 522 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
409 
      ierr = nf90_inquire_variable(fid, varid,
523 
      ierr = nf90_inquire_variable(fid, varid,
410 
     >                                   dimids = dimids(1:ndim))
 524 
     >                                   dimids = dimids(1:ndim))
411 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 525 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
412 
      do i=1,ndim
 526 
      do i=1,ndim
413 
           ierr = nf90_inquire_dimension(fid, dimids(i),
527 
           ierr = nf90_inquire_dimension(fid, dimids(i),
414 
     >                               name = dimname(i) )
 528 
     >                               name = dimname(i) )
415 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 529 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
416 
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
 530 
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
417 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 531 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
418 
      enddo
 532 
      enddo
419 
 
 533 
 
420 
c     Check whether the list of dimensions is OK
 534 
c     Check whether the list of dimensions is OK
421 
      if ( ( dimname(1).ne.'lon'  ).or.
 535 
      if ( ( dimname(1).ne.'lon'  ).or.
422 
     >     ( dimname(2).ne.'lat'  ).or.
 536 
     >     ( dimname(2).ne.'lat'  ).or.
423 
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
 537 
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
424 
     >     ( dimname(4).ne.'time' ) )
 538 
     >     ( dimname(4).ne.'time' ) )
425 
     >then
 539 
     >then
426 
        print*,' ERROR: the dimensions of the variable are not correct'
 540 
        print*,' ERROR: the dimensions of the variable are not correct'
427 
        print*,'        expected -> lon / lat / lev / time'
 541 
        print*,'        expected -> lon / lat / lev / time'
428 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
 542 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
429 
        stop
 543 
        stop
430 
      endif
 544 
      endif
431 
 
 545 
 
- 
 
 546 
c     Get dimension of AK,BK
- 
 
 547 
      varname = 'nhym'
- 
 
 548 
      ierr = NF90_INQ_DIMID(fid,varname,dimid)
- 
 
 549 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 550 
      ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
- 
 
 551 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 552 
 
- 
 
 553 
c     Check about the type of vertical levels
- 
 
 554 
      varname=dimname(3)
- 
 
 555 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 556 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 557 
      ierr = nf90_get_att(fid, varid, "standard_name", leveltype)
- 
 
 558 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 559 
      if ( (leveltype.ne.'hybrid_sigma_pressure').and.
- 
 
 560 
     >     (leveltype.ne.'air_pressure'         ) )
- 
 
 561 
     >then
- 
 
 562 
         print*,' ERROR: input netCDF data must be on hybrid-sigma'
- 
 
 563 
         print*,'        or air pressure levels!',trim(leveltype)
- 
 
 564 
         stop
- 
 
 565 
      endif
- 
 
 566 
 
- 
 
 567 
c     Check that vardim(3)==#AK,BK for hybrid-sigmal levels
- 
 
 568 
      if ( (leveltype.eq.'hybrid_sigma_pressure').and.
- 
 
 569 
     >     (dimname(3).eq.'lev') )
- 
 
 570 
     >then
- 
 
 571 
        if ( nakbktmp.ne.vardim(3) ) then
- 
 
 572 
           print*,' ERROR: for hybrid-sigma pressure levels, #AK,BK'
- 
 
 573 
           print*,'        must agree with number of vertical levels'
- 
 
 574 
           print*,'        ',vardim(3),nakbktmp
- 
 
 575 
           stop
- 
 
 576 
        endif
- 
 
 577 
      endif
- 
 
 578 
 
432 
c     Allocate memory for reading arrays - depending on <closear>
 579 
c     Allocate memory for reading arrays - depending on <closear>
433 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
 580 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
434 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
 581 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
435 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
 582 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
436 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
 583 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
437 
      allocate(lon(vardim(1)),stat=stat)
 584 
      allocate(lon(vardim(1)),stat=stat)
438 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
 585 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
439 
      allocate(lat(vardim(2)),stat=stat)
 586 
      allocate(lat(vardim(2)),stat=stat)
440 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
 587 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
- 
 
 588 
      allocate(lev(vardim(3)),stat=stat)
- 
 
 589 
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
- 
 
 590 
      allocate(aktmp(nakbktmp),stat=stat)
- 
 
 591 
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
- 
 
 592 
      allocate(bktmp(nakbktmp),stat=stat)
- 
 
 593 
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
441 
 
 594 
 
442 
c     Get domain boundaries
 595 
c     Get domain boundaries - longitude, latitude, levels
443 
      varname = dimname(1)
 596 
      varname = dimname(1)
444 
      ierr = NF90_INQ_VARID(fid,varname,varid)
 597 
      ierr = NF90_INQ_VARID(fid,varname,varid)
445 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 598 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
446 
      ierr = nf90_get_var(fid,varid,lon)
 599 
      ierr = nf90_get_var(fid,varid,lon)
447 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 600 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
448 
      varname = dimname(2)
 601 
      varname = dimname(2)
449 
      ierr = NF90_INQ_VARID(fid,varname,varid)
 602 
      ierr = NF90_INQ_VARID(fid,varname,varid)
450 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 603 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
451 
      ierr = nf90_get_var(fid,varid,lat)
 604 
      ierr = nf90_get_var(fid,varid,lat)
452 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 605 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 606 
      varname = dimname(3)
- 
 
 607 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 608 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 609 
      ierr = nf90_get_var(fid,varid,lev)
- 
 
 610 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 611 
 
- 
 
 612 
c     Get ak and bk
- 
 
 613 
      varname='hyam'
- 
 
 614 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 615 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 616 
      ierr = nf90_get_var(fid,varid,aktmp)
- 
 
 617 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 618 
      varname='hybm'
- 
 
 619 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 620 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 621 
      ierr = nf90_get_var(fid,varid,bktmp)
- 
 
 622 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 623 
 
- 
 
 624 
c     Check that unit of ak is in hPa - if necessary correct it
- 
 
 625 
      varname='hyam'
- 
 
 626 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 627 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 628 
      ierr = nf90_get_att(fid, varid, "units", units)
- 
 
 629 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 630 
      if ( units.eq.'Pa' ) then
- 
 
 631 
         do k=1,nakbktmp
- 
 
 632 
            aktmp(k) = 0.01 * aktmp(k)
- 
 
 633 
         enddo
- 
 
 634 
      endif
- 
 
 635 
 
- 
 
 636 
c     Check that unit of lev is in hPa - if necessary correct it
- 
 
 637 
      if ( leveltype.eq.'air_pressure' ) then
- 
 
 638 
         varname='lev'
- 
 
 639 
         ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 640 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 641 
         ierr = nf90_get_att(fid, varid, "units", units)
- 
 
 642 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 643 
         if ( units.eq.'Pa' ) then
- 
 
 644 
            do k=1,vardim(3)
- 
 
 645 
               lev(k) = 0.01 * lev(k)
- 
 
 646 
            enddo
- 
 
 647 
         endif
- 
 
 648 
      endif
- 
 
 649 
 
- 
 
 650 
c     Decide whether to swap vertical levels
- 
 
 651 
      vertical_swap = 1
- 
 
 652 
      if ( leveltype.eq.'hybrid_sigma_pressure') then
- 
 
 653 
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
- 
 
 654 
     >       (aktmp(2) + bktmp(2) * 1000.) )
- 
 
 655 
     >  then
- 
 
 656 
          vertical_swap = 0
- 
 
 657 
        endif
- 
 
 658 
      elseif ( leveltype.eq.'air_pressure') then
- 
 
 659 
        if ( lev(1).gt.lev(2) ) then
- 
 
 660 
          vertical_swap = 0
- 
 
 661 
        endif
- 
 
 662 
      endif
- 
 
 663 
c      print*,' Vertical SWAP ',trim(fieldname),' -> ', vertical_swap
453 
 
 664 
 
454 
c     Read data
665 
c     Read data
455 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
 666 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
456 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 667 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
457 
      ierr = nf90_get_var(fid,varid,tmp3)
 668 
      ierr = nf90_get_var(fid,varid,tmp3)
458 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 669 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
459
670
460 
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
 671 
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
461 
      if ( vardim(3).eq.1 ) then
 672 
      if ( vardim(3).eq.1 ) then
462 
         do i=1,vardim(1)
 673 
         do i=1,vardim(1)
463 
            do j=1,vardim(2)
 674 
            do j=1,vardim(2)
464 
               do k=1,vardim(3)
 675 
               do k=1,vardim(3)
465 
                  tmp3(i,j,k) = tmp3(i,j,1)
 676 
                  tmp3(i,j,k) = tmp3(i,j,1)
466 
               enddo
 677 
               enddo
467 
            enddo
 678 
            enddo
468 
         enddo
 679 
         enddo
469 
      endif
 680 
      endif
470 
 
 681 
 
471 
c     Save the ouput array - close on request
 682 
c     Decide whether to close arrays
472 
      delta = varmax(1)-varmin(1)-360.
 683 
      delta = varmax(1)-varmin(1)-360.
473 
      if (abs(delta+dx).lt.eps) then
 684 
      if (abs(delta+dx).lt.eps) then
474 
          closear = 1
 685 
          closear = 1
475 
      else
 686 
      else
476 
          closear = 0
 687 
          closear = 0
477 
      endif
 688 
      endif
478 
 
 689 
 
- 
 
 690 
c     Save output array - close array and swap on demand
479 
      do j=1,vardim(2)
 691 
      do j=1,vardim(2)
480 
        do k=1,vardim(3)
 692 
        do k=1,vardim(3)
481 
          do i=1,vardim(1)
 693 
          do i=1,vardim(1)
- 
 
 694 
             if ( vertical_swap.eq.1 ) then
482 
             field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
 695 
                 field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
- 
 
 696 
             else
- 
 
 697 
                 field(i,j,k) = tmp3(i,j,k)
- 
 
 698 
             endif
483 
          enddo
 699 
          enddo
484 
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
 700 
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
485 
        enddo
 701 
        enddo
486 
      enddo
 702 
      enddo
487
703
488 
c     Exit point
 704 
c     Exit point
489 
      return
 705 
      return
490
706
491 
      end
 707 
      end
492 
 
 708 
 
493 
c     ------------------------------------------------------------
 709 
c     ------------------------------------------------------------
494 
c     Close input file
 710 
c     Close input file
495 
c     ------------------------------------------------------------
 711 
c     ------------------------------------------------------------
496 
 
 712 
 
497 
      subroutine input_close(fid)
 713 
      subroutine input_close(fid)
498 
 
 714 
 
499 
c     Close the input file with file identifier <fid>.
 715 
c     Close the input file with file identifier <fid>.
500 
 
 716 
 
501 
      use netcdf
 717 
      use netcdf
502 
 
 718 
 
503 
      implicit none
 719 
      implicit none
504 
 
 720 
 
505 
c     Declaration of subroutine parameters
 721 
c     Declaration of subroutine parameters
506 
      integer fid
 722 
      integer fid
507 
 
 723 
 
508 
c     Auxiliary variables
 724 
c     Auxiliary variables
509 
      integer ierr
 725 
      integer ierr
510 
 
 726 
 
511 
c     Close file
 727 
c     Close file
512 
      ierr = NF90_CLOSE(fid)
 728 
      ierr = NF90_CLOSE(fid)
513 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 729 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
514
730
515 
      end
 731 
      end
516
732
517 
c     ------------------------------------------------------------
 733 
c     ------------------------------------------------------------
518 
c     Get a list of variables on netCDF file
 734 
c     Get a list of variables on netCDF file
519 
c     ------------------------------------------------------------
 735 
c     ------------------------------------------------------------
520 
 
 736 
 
521 
      subroutine input_getvars(fid,vnam,nvars)
 737 
      subroutine input_getvars(fid,vnam,nvars)
522 
 
 738 
 
523 
c     List of variables on netCDF file
 739 
c     List of variables on netCDF file
524 
 
 740 
 
525 
      use netcdf
 741 
      use netcdf
526 
 
 742 
 
527 
      implicit none
 743 
      implicit none
528 
 
 744 
 
529 
c     Declaration of subroutine parameters
 745 
c     Declaration of subroutine parameters
530 
      integer      fid
 746 
      integer      fid
531 
      integer      nvars
 747 
      integer      nvars
532 
      character*80 vnam(200)
 748 
      character*80 vnam(200)
533 
 
 749 
 
534 
c     Auxiliary variables
 750 
c     Auxiliary variables
535 
      integer ierr
 751 
      integer ierr
536 
      integer i
 752 
      integer i
537 
      integer nDims,nGlobalAtts,unlimdimid
 753 
      integer nDims,nGlobalAtts,unlimdimid
538 
 
 754 
 
539 
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
 755 
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
540 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 756 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
541
757
542 
      do i=1,nVars
 758 
      do i=1,nVars
543 
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
 759 
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
544 
      enddo
 760 
      enddo
545
761
546 
      end
 762 
      end