Subversion Repositories lagranto.ecmwf

Rev

Rev 25 | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 25 Rev 27
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 
 
177 
c     Allocate memory for reading arrays
 189 
c     Allocate memory for reading arrays
178 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
 190 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
179 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
 191 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
180 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
 192 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
181 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
 193 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
182 
      allocate(lon(vardim(1)),stat=stat)
 194 
      allocate(lon(vardim(1)),stat=stat)
183 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
195 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
184 
      allocate(lat(vardim(2)),stat=stat)
 196 
      allocate(lat(vardim(2)),stat=stat)
185 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
197 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
- 
 
 198 
      allocate(lev(vardim(3)),stat=stat)
- 
 
 199 
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
186 
      allocate(times(vardim(4)),stat=stat)
 200 
      allocate(times(vardim(4)),stat=stat)
187 
      if (stat.ne.0) print*,'*** error allocating array times   ***'
 201 
      if (stat.ne.0) print*,'*** error allocating array times   ***'
188 
      allocate(aktmp(vardim(3)),stat=stat)
 202 
      allocate(aktmp(nakbktmp),stat=stat)
189 
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
 203 
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
190 
      allocate(bktmp(vardim(3)),stat=stat)
 204 
      allocate(bktmp(nakbktmp),stat=stat)
191 
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
 205 
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
192 
 
 206 
 
193 
c     Get domain longitudes and latitudes
 207 
c     Get domain longitudes, latitudes and levels
194 
      varname = dimname(1)
 208 
      varname = dimname(1)
195 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 209 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
196 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 210 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
197 
      ierr = nf90_get_var(cdfid,varid,lon)
 211 
      ierr = nf90_get_var(cdfid,varid,lon)
198 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 212 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
199 
      varname = dimname(2)
 213 
      varname = dimname(2)
200 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 214 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
201 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 215 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
202 
      ierr = nf90_get_var(cdfid,varid,lat)
 216 
      ierr = nf90_get_var(cdfid,varid,lat)
203 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 217 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 218 
      varname = dimname(3)
- 
 
 219 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 220 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 221 
      ierr = nf90_get_var(cdfid,varid,lev)
- 
 
 222 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
204
223
205 
c     Get ak and bk
 224 
c     Get ak and bk
206 
      varname='hyam'
 225 
      varname='hyam'
207 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 226 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
208 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 227 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
209 
      ierr = nf90_get_var(cdfid,varid,aktmp)
 228 
      ierr = nf90_get_var(cdfid,varid,aktmp)
210 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 229 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
211 
      varname='hybm'
 230 
      varname='hybm'
212 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 231 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
213 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 232 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
214 
      ierr = nf90_get_var(cdfid,varid,bktmp)
 233 
      ierr = nf90_get_var(cdfid,varid,bktmp)
215 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 234 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
216 
 
 235 
 
217 
c     Check that unit of ak is in hPa - if necessary correct it
 236 
c     Check that unit of ak is in hPa - if necessary correct it
218 
      varname='hyam'
 237 
      varname='hyam'
219 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 238 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
220 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 239 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
221 
      ierr = nf90_get_att(cdfid, varid, "units", units)
 240 
      ierr = nf90_get_att(cdfid, varid, "units", units)
222 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 241 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
223 
      if ( units.eq.'Pa' ) then
 242 
      if ( units.eq.'Pa' ) then
224 
         do k=1,vardim(3)
 243 
         do k=1,nakbktmp
225 
            aktmp(k) = 0.01 * aktmp(k)
 244 
            aktmp(k) = 0.01 * aktmp(k)
226 
         enddo
 245 
         enddo
227 
      endif
 246 
      endif
228 
 
 247 
 
- 
 
 248 
c     Decide whether to swap vertical levels - highest pressure at index 1
- 
 
 249 
      vertical_swap = 1
- 
 
 250 
      if ( (aktmp(1) + bktmp(1) * 1000.).gt.
- 
 
 251 
     >       (aktmp(2) + bktmp(2) * 1000.) )
- 
 
 252 
     >then
- 
 
 253 
          vertical_swap = 0
- 
 
 254 
      endif
- 
 
 255 
 
229 
c     Get time information (check if time is correct)
 256 
c     Get time information (check if time is correct)
230 
      varname = 'time'
 257 
      varname = 'time'
231 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 258 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
232 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 259 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
233 
      ierr = nf90_get_var(cdfid,varid,times)
 260 
      ierr = nf90_get_var(cdfid,varid,times)
234 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 261 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
235 
      isok=0
 262 
      isok=0
236 
      do i=1,vardim(4)
 263 
      do i=1,vardim(4)
237 
        if (abs(time-times(i)).lt.eps) then
 264 
        if (abs(time-times(i)).lt.eps) then
238 
               isok = 1
 265 
               isok = 1
239 
               rtime = times(i)
 266 
               rtime = times(i)
240 
        elseif (timecheck.eq.'no') then
 267 
        elseif (timecheck.eq.'no') then
241 
               isok = 1
 268 
               isok = 1
242 
               rtime = times(1)
 269 
               rtime = times(1)
243 
        endif
 270 
        endif
244 
      enddo
 271 
      enddo
245 
      if ( isok.eq.0 ) then
 272 
      if ( isok.eq.0 ) then
246 
         print*,' ERROR: time ',rtime,' not found on netCDF file'
273 
         print*,' ERROR: time ',rtime,' not found on netCDF file'
247 
         stop
 274 
         stop
248 
      endif
 275 
      endif
249 
 
 276 
 
250 
c     Read surface pressure
 277 
c     Read surface pressure
251 
      varname='PS'
 278 
      varname='PS'
252 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 279 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
253 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 280 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
254 
      ierr = nf90_get_var(cdfid,varid,tmp2)
 281 
      ierr = nf90_get_var(cdfid,varid,tmp2)
255 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 282 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
256
283
257 
c     Check that surface pressure is in hPa
 284 
c     Check that surface pressure is in hPa
258 
      maxps = -1.e39
 285 
      maxps = -1.e19
259 
      minps =  1.e39
 286 
      minps =  1.e19
260 
      do i=1,vardim(1)
 287 
      do i=1,vardim(1)
261 
        do j=1,vardim(2)
 288 
        do j=1,vardim(2)
262 
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
 289 
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
263 
             if (tmp2(i,j).lt.minps) minps = tmp2(i,j)
 290 
             if (tmp2(i,j).lt.minps) minps = tmp2(i,j)
264 
        enddo
 291 
        enddo
265 
      enddo
 292 
      enddo
266 
      if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then
 293 
      if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then
267 
         print*,' ERROR: surface pressre PS must be in hPa'
 294 
         print*,' ERROR: surface pressre PS must be in hPa'
268 
         print*,'       ',maxps,minps
 295 
         print*,'       ',maxps,minps
269 
         stop
 296 
         stop
270 
      endif
 297 
      endif
271 
 
 298 
 
- 
 
 299 
c     ---- Define output of subroutine --------------------------------
- 
 
 300 
 
- 
 
 301 
c     If not full list of vertical levels, reduce AK,BK arrays
272 
c     Calculate layer and level pressures
 302 
      if ( (leveltype.eq.'hybrid_sigma_pressure').and.
273 
      do i=1,vardim(1)
 303 
     >     (nakbktmp.ne.vardim(3) ) )
- 
 
 304 
     >then
- 
 
 305 
         print*,' WARNING: only subset of vertical levels used...'
274 
         do j=1,vardim(2)
 306 
         do k=1,vardim(3)
- 
 
 307 
            if ( vertical_swap.eq.1 ) then
275 
               do k=1,vardim(3)
 308 
               aktmp(k) = aktmp( k+nakbktmp-vardim(3) )
276 
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
 309 
               bktmp(k) = bktmp( k+nakbktmp-vardim(3) )
277 
               enddo
 310 
            endif
278 
         enddo
 311 
         enddo
279 
      enddo
 312 
      endif
280 
 
 313 
 
281 
c     Set the grid dimensions and constants
 314 
c     Set the grid dimensions and constants
282 
      nx      = vardim(1)
 315 
      nx      = vardim(1)
283 
      ny      = vardim(2)
 316 
      ny      = vardim(2)
284 
      nz      = vardim(3)
 317 
      nz      = vardim(3)
285 
      xmin    = lon(1)
 318 
      xmin    = lon(1)
286 
      ymin    = lat(1)
 319 
      ymin    = lat(1)
287 
      xmax    = lon(nx)
 320 
      xmax    = lon(nx)
288 
      ymax    = lat(ny)
 321 
      ymax    = lat(ny)
289 
      dx      = (xmax-xmin)/real(nx-1)
 322 
      dx      = (xmax-xmin)/real(nx-1)
290 
      dy      = (ymax-ymin)/real(ny-1)
 323 
      dy      = (ymax-ymin)/real(ny-1)
291 
      pollon  = 0.
 324 
      pollon  = 0.
292 
      pollat  = 90.
 325 
      pollat  = 90.
293 
      stagz   = 0.
 326 
      stagz   = 0.
294 
      delta   = xmax-xmin-360.
 327 
      delta   = xmax-xmin-360.
295 
      if (abs(delta+dx).lt.eps) then
 328 
      if (abs(delta+dx).lt.eps) then
296 
          xmax    = xmax + dx
 329 
          xmax    = xmax + dx
297 
          nx      = nx + 1
 330 
          nx      = nx + 1
298 
          closear = 1
 331 
          closear = 1
299 
      else
 332 
      else
300 
          closear = 0
 333 
          closear = 0
301 
      endif
 334 
      endif
302 
 
 335 
 
303 
c     Save the output arrays (if fid>0) - close arrays on request
 336 
c     Save the output arrays (if fid>0) - close arrays on request
304 
      if ( fid.gt.0 ) then
 337 
      if ( fid.gt.0 ) then
305 
 
 338 
 
- 
 
 339 
c        Calculate layer pressures
- 
 
 340 
         do i=1,vardim(1)
- 
 
 341 
              do j=1,vardim(2)
- 
 
 342 
                 do k=1,vardim(3)
- 
 
 343 
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
- 
 
 344 
                 enddo
- 
 
 345 
              enddo
- 
 
 346 
         enddo
- 
 
 347 
 
- 
 
 348 
c        Get PS - close array on demand
306 
         do j=1,vardim(2)
 349 
         do j=1,vardim(2)
307 
           do i=1,vardim(1)
 350 
           do i=1,vardim(1)
308 
             ps(i,j) = tmp2(i,j)
 351 
             ps(i,j) = tmp2(i,j)
309 
           enddo
 352 
           enddo
310 
           if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)
 353 
           if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)
311 
         enddo
 354 
         enddo
312 
 
 355 
 
- 
 
 356 
c        Get P3 - close array on demand + vertical swap
313 
         do j=1,vardim(2)
 357 
         do j=1,vardim(2)
314 
           do k=1,vardim(3)
 358 
           do k=1,vardim(3)
315 
             do i=1,vardim(1)
 359 
             do i=1,vardim(1)
- 
 
 360 
               if ( vertical_swap.eq.1 ) then
316 
               p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)
 361 
                  p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)
- 
 
 362 
               else
- 
 
 363 
                  p3(i,j,k) = tmp3(i,j,k)
- 
 
 364 
               endif
317 
             enddo
 365 
             enddo
318 
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
 366 
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
319 
           enddo
 367 
           enddo
320 
         enddo
 368 
         enddo
321 
 
 369 
 
- 
 
 370 
c        Get AK,BK - vertical swap on demand
322 
         do k=1,vardim(3)
 371 
         do k=1,vardim(3)
- 
 
 372 
              if ( vertical_swap.eq.1 ) then
323 
            ak(k) = aktmp(vardim(3)-k+1)
 373 
                 ak(k) = aktmp(vardim(3)-k+1)
324 
            bk(k) = bktmp(vardim(3)-k+1)
 374 
                 bk(k) = bktmp(vardim(3)-k+1)
- 
 
 375 
              endif
325 
         enddo
 376 
         enddo
326 
 
 377 
 
327 
      endif
 378 
      endif
328 
 
 379 
 
- 
 
 380 
 
329 
      return
 381 
      return
330
382
331 
      end
 383 
      end
332 
 
 384 
 
333 
c     ------------------------------------------------------------
 385 
c     ------------------------------------------------------------
334 
c     Read wind information
 386 
c     Read wind information
335 
c     ------------------------------------------------------------
 387 
c     ------------------------------------------------------------
336 
 
 388 
 
337 
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
 389 
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
338 
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
 390 
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
339 
     >                       timecheck)
 391 
     >                       timecheck)
340 
 
 392 
 
341 
c     Read the wind component <fieldname> from the file with identifier
 393 
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
394 
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
 395 
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
 396 
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,
 397 
c     performed to have an agreement with the grid specified by <xmin,xmax,
346 
c     ymin,ymax,dx,dy,nx,ny,nz>.
 398 
c     ymin,ymax,dx,dy,nx,ny,nz>.
347 
 
 399 
 
348 
      use netcdf
 400 
      use netcdf
349 
 
 401 
 
350 
      implicit none
 402 
      implicit none
351 
 
 403 
 
352 
c     Declaration of variables and parameters
 404 
c     Declaration of variables and parameters
353 
      integer      fid                 ! File identifier
 405 
      integer      fid                 ! File identifier
354 
      character*80 fieldname           ! Name of the wind field
 406 
      character*80 fieldname           ! Name of the wind field
355 
      integer      nx,ny,nz            ! Dimension of fields
 407 
      integer      nx,ny,nz            ! Dimension of fields
356 
      real         field(nx,ny,nz)     ! 3d wind field
 408 
      real         field(nx,ny,nz)     ! 3d wind field
357 
      real         stagz               ! Staggering in the z direction
 409 
      real         stagz               ! Staggering in the z direction
358 
      real         mdv                 ! Missing data flag
 410 
      real         mdv                 ! Missing data flag
359 
      real         xmin,xmax,ymin,ymax ! Domain size
 411 
      real         xmin,xmax,ymin,ymax ! Domain size
360 
      real         dx,dy               ! Horizontal resolution
 412 
      real         dx,dy               ! Horizontal resolution
361 
      real         time                ! Time
 413 
      real         time                ! Time
362 
      character*80 timecheck           ! Either 'yes' or 'no'
 414 
      character*80 timecheck           ! Either 'yes' or 'no'
363 
 
 415 
 
364 
c     Numerical and physical parameters
 416 
c     Numerical and physical parameters
365 
      real        eps                 ! Numerical epsilon
 417 
      real        eps                 ! Numerical epsilon
366 
      parameter  (eps=0.001)
 418 
      parameter  (eps=0.001)
367 
      real        notimecheck         ! 'Flag' for no time check
 419 
      real        notimecheck         ! 'Flag' for no time check
368 
      parameter  (notimecheck=7.26537)
 420 
      parameter  (notimecheck=7.26537)
369 
 
 421 
 
370 
c     Netcdf variables
 422 
c     Netcdf variables
371 
      integer      ierr
 423 
      integer      ierr
372 
      character*80 varname
 424 
      character*80 varname
373 
      integer      vardim(4)
 425 
      integer      vardim(4)
374 
      real         varmin(4),varmax(4)
 426 
      real         varmin(4),varmax(4)
375 
      real         stag(4)
 427 
      real         stag(4)
376 
      integer      ndim
 428 
      integer      ndim
377 
      real         times(10)
 429 
      real         times(10)
378 
      integer      ntimes
 430 
      integer      ntimes
379 
      character*80 cstfile
 431 
      character*80 cstfile
380 
      integer      cstid
 432 
      integer      cstid
381 
      real         aklay(200),bklay(200),aklev(200),bklev(200)
 433 
      real         aklay(200),bklay(200),aklev(200),bklev(200)
382 
      real         ps(nx,ny)
 434 
      real         ps(nx,ny)
383 
      integer      dimids (nf90_max_var_dims)
 435 
      integer      dimids (nf90_max_var_dims)
384 
      character*80 dimname(nf90_max_var_dims)
 436 
      character*80 dimname(nf90_max_var_dims)
385 
      integer      varid
 437 
      integer      varid
386 
      integer      cdfid
 438 
      integer      cdfid
387 
      real,allocatable, dimension (:)     :: lon,lat
 439 
      real,allocatable, dimension (:)     :: lon,lat,lev
388 
      real,allocatable, dimension (:,:)   :: tmp2
 440 
      real,allocatable, dimension (:,:)   :: tmp2
389 
      real,allocatable, dimension (:,:,:) :: tmp3
 441 
      real,allocatable, dimension (:,:,:) :: tmp3
- 
 
 442 
      real,allocatable, dimension (:)     :: aktmp,bktmp
- 
 
 443 
      character*80  leveltype
- 
 
 444 
      integer       vertical_swap
- 
 
 445 
      character*80  units
- 
 
 446 
      integer       nakbktmp
- 
 
 447 
      integer       dimid
390 
 
 448 
 
391 
c     Auxiliary variables
 449 
c     Auxiliary variables
392 
      integer      isok
 450 
      integer      isok
393 
      integer      i,j,k
 451 
      integer      i,j,k
394 
      integer      nz1
 452 
      integer      nz1
395 
      real         rtime
 453 
      real         rtime
396 
      integer      closear
 454 
      integer      closear
397 
      integer      stat
 455 
      integer      stat
398 
      real         delta
 456 
      real         delta
399 
 
 457 
 
400 
c     Get the number of dimensions -> ndim
 458 
c     Get the number of dimensions -> ndim
401 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
 459 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
402 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 460 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
403 
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
 461 
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
404 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 462 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
405 
 
 463 
 
406 
c     Get the dimensions of the arrays -> varid(1:ndim)
 464 
c     Get the dimensions of the arrays -> varid(1:ndim)
407 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
 465 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
408 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 466 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
409 
      ierr = nf90_inquire_variable(fid, varid,
467 
      ierr = nf90_inquire_variable(fid, varid,
410 
     >                                   dimids = dimids(1:ndim))
 468 
     >                                   dimids = dimids(1:ndim))
411 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 469 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
412 
      do i=1,ndim
 470 
      do i=1,ndim
413 
           ierr = nf90_inquire_dimension(fid, dimids(i),
471 
           ierr = nf90_inquire_dimension(fid, dimids(i),
414 
     >                               name = dimname(i) )
 472 
     >                               name = dimname(i) )
415 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 473 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
416 
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
 474 
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
417 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 475 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
418 
      enddo
 476 
      enddo
419 
 
 477 
 
420 
c     Check whether the list of dimensions is OK
 478 
c     Check whether the list of dimensions is OK
421 
      if ( ( dimname(1).ne.'lon'  ).or.
 479 
      if ( ( dimname(1).ne.'lon'  ).or.
422 
     >     ( dimname(2).ne.'lat'  ).or.
 480 
     >     ( dimname(2).ne.'lat'  ).or.
423 
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
 481 
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
424 
     >     ( dimname(4).ne.'time' ) )
 482 
     >     ( dimname(4).ne.'time' ) )
425 
     >then
 483 
     >then
426 
        print*,' ERROR: the dimensions of the variable are not correct'
 484 
        print*,' ERROR: the dimensions of the variable are not correct'
427 
        print*,'        expected -> lon / lat / lev / time'
 485 
        print*,'        expected -> lon / lat / lev / time'
428 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
 486 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
429 
        stop
 487 
        stop
430 
      endif
 488 
      endif
431 
 
 489 
 
- 
 
 490 
c     Get dimension of AK,BK
- 
 
 491 
      varname = 'nhym'
- 
 
 492 
      ierr = NF90_INQ_DIMID(fid,varname,dimid)
- 
 
 493 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 494 
      ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
- 
 
 495 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 496 
 
432 
c     Allocate memory for reading arrays - depending on <closear>
 497 
c     Allocate memory for reading arrays - depending on <closear>
433 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
 498 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
434 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
 499 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
435 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
 500 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
436 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
 501 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
437 
      allocate(lon(vardim(1)),stat=stat)
 502 
      allocate(lon(vardim(1)),stat=stat)
438 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
 503 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
439 
      allocate(lat(vardim(2)),stat=stat)
 504 
      allocate(lat(vardim(2)),stat=stat)
440 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
 505 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
- 
 
 506 
      allocate(lev(vardim(3)),stat=stat)
- 
 
 507 
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
- 
 
 508 
      allocate(aktmp(nakbktmp),stat=stat)
- 
 
 509 
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
- 
 
 510 
      allocate(bktmp(nakbktmp),stat=stat)
- 
 
 511 
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
441 
 
 512 
 
442 
c     Get domain boundaries
 513 
c     Get domain boundaries - longitude, latitude, levels
443 
      varname = dimname(1)
 514 
      varname = dimname(1)
444 
      ierr = NF90_INQ_VARID(fid,varname,varid)
 515 
      ierr = NF90_INQ_VARID(fid,varname,varid)
445 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 516 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
446 
      ierr = nf90_get_var(fid,varid,lon)
 517 
      ierr = nf90_get_var(fid,varid,lon)
447 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 518 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
448 
      varname = dimname(2)
 519 
      varname = dimname(2)
449 
      ierr = NF90_INQ_VARID(fid,varname,varid)
 520 
      ierr = NF90_INQ_VARID(fid,varname,varid)
450 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 521 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
451 
      ierr = nf90_get_var(fid,varid,lat)
 522 
      ierr = nf90_get_var(fid,varid,lat)
452 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 523 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 524 
      varname = dimname(3)
- 
 
 525 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 526 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 527 
      ierr = nf90_get_var(fid,varid,lev)
- 
 
 528 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 529 
 
- 
 
 530 
c     Get ak and bk
- 
 
 531 
      varname='hyam'
- 
 
 532 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 533 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 534 
      ierr = nf90_get_var(fid,varid,aktmp)
- 
 
 535 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 536 
      varname='hybm'
- 
 
 537 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 538 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 539 
      ierr = nf90_get_var(fid,varid,bktmp)
- 
 
 540 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 541 
 
- 
 
 542 
c     Check that unit of ak is in hPa - if necessary correct it
- 
 
 543 
      varname='hyam'
- 
 
 544 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 545 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 546 
      ierr = nf90_get_att(fid, varid, "units", units)
- 
 
 547 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 548 
      if ( units.eq.'Pa' ) then
- 
 
 549 
         do k=1,nakbktmp
- 
 
 550 
            aktmp(k) = 0.01 * aktmp(k)
- 
 
 551 
         enddo
- 
 
 552 
      endif
- 
 
 553 
 
- 
 
 554 
c     Decide whether to swap vertical levels
- 
 
 555 
      vertical_swap = 1
- 
 
 556 
      if ( (aktmp(1) + bktmp(1) * 1000.).gt.
- 
 
 557 
     >     (aktmp(2) + bktmp(2) * 1000.) )
- 
 
 558 
     >then
- 
 
 559 
          vertical_swap = 0
- 
 
 560 
      endif
453 
 
 561 
 
454 
c     Read data
562 
c     Read data
455 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
 563 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
456 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 564 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
457 
      ierr = nf90_get_var(fid,varid,tmp3)
 565 
      ierr = nf90_get_var(fid,varid,tmp3)
458 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 566 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
459
567
460 
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
 568 
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
461 
      if ( vardim(3).eq.1 ) then
 569 
      if ( vardim(3).eq.1 ) then
462 
         do i=1,vardim(1)
 570 
         do i=1,vardim(1)
463 
            do j=1,vardim(2)
 571 
            do j=1,vardim(2)
464 
               do k=1,vardim(3)
 572 
               do k=1,vardim(3)
465 
                  tmp3(i,j,k) = tmp3(i,j,1)
 573 
                  tmp3(i,j,k) = tmp3(i,j,1)
466 
               enddo
 574 
               enddo
467 
            enddo
 575 
            enddo
468 
         enddo
 576 
         enddo
469 
      endif
 577 
      endif
470 
 
 578 
 
471 
c     Save the ouput array - close on request
 579 
c     Decide whether to close arrays
472 
      delta = varmax(1)-varmin(1)-360.
 580 
      delta = varmax(1)-varmin(1)-360.
473 
      if (abs(delta+dx).lt.eps) then
 581 
      if (abs(delta+dx).lt.eps) then
474 
          closear = 1
 582 
          closear = 1
475 
      else
 583 
      else
476 
          closear = 0
 584 
          closear = 0
477 
      endif
 585 
      endif
478 
 
 586 
 
- 
 
 587 
c     Save output array - close array and swap on demand
479 
      do j=1,vardim(2)
 588 
      do j=1,vardim(2)
480 
        do k=1,vardim(3)
 589 
        do k=1,vardim(3)
481 
          do i=1,vardim(1)
 590 
          do i=1,vardim(1)
- 
 
 591 
             if ( vertical_swap.eq.1 ) then
482 
             field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
 592 
                 field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
- 
 
 593 
             else
- 
 
 594 
                 field(i,j,k) = tmp3(i,j,k)
- 
 
 595 
             endif
483 
          enddo
 596 
          enddo
484 
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
 597 
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
485 
        enddo
 598 
        enddo
486 
      enddo
 599 
      enddo
487
600
488 
c     Exit point
 601 
c     Exit point
489 
      return
 602 
      return
490
603
491 
      end
 604 
      end
492 
 
 605 
 
493 
c     ------------------------------------------------------------
 606 
c     ------------------------------------------------------------
494 
c     Close input file
 607 
c     Close input file
495 
c     ------------------------------------------------------------
 608 
c     ------------------------------------------------------------
496 
 
 609 
 
497 
      subroutine input_close(fid)
 610 
      subroutine input_close(fid)
498 
 
 611 
 
499 
c     Close the input file with file identifier <fid>.
 612 
c     Close the input file with file identifier <fid>.
500 
 
 613 
 
501 
      use netcdf
 614 
      use netcdf
502 
 
 615 
 
503 
      implicit none
 616 
      implicit none
504 
 
 617 
 
505 
c     Declaration of subroutine parameters
 618 
c     Declaration of subroutine parameters
506 
      integer fid
 619 
      integer fid
507 
 
 620 
 
508 
c     Auxiliary variables
 621 
c     Auxiliary variables
509 
      integer ierr
 622 
      integer ierr
510 
 
 623 
 
511 
c     Close file
 624 
c     Close file
512 
      ierr = NF90_CLOSE(fid)
 625 
      ierr = NF90_CLOSE(fid)
513 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 626 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
514
627
515 
      end
 628 
      end
516
629
517 
c     ------------------------------------------------------------
 630 
c     ------------------------------------------------------------
518 
c     Get a list of variables on netCDF file
 631 
c     Get a list of variables on netCDF file
519 
c     ------------------------------------------------------------
 632 
c     ------------------------------------------------------------
520 
 
 633 
 
521 
      subroutine input_getvars(fid,vnam,nvars)
 634 
      subroutine input_getvars(fid,vnam,nvars)
522 
 
 635 
 
523 
c     List of variables on netCDF file
 636 
c     List of variables on netCDF file
524 
 
 637 
 
525 
      use netcdf
 638 
      use netcdf
526 
 
 639 
 
527 
      implicit none
 640 
      implicit none
528 
 
 641 
 
529 
c     Declaration of subroutine parameters
 642 
c     Declaration of subroutine parameters
530 
      integer      fid
 643 
      integer      fid
531 
      integer      nvars
 644 
      integer      nvars
532 
      character*80 vnam(200)
 645 
      character*80 vnam(200)
533 
 
 646 
 
534 
c     Auxiliary variables
 647 
c     Auxiliary variables
535 
      integer ierr
 648 
      integer ierr
536 
      integer i
 649 
      integer i
537 
      integer nDims,nGlobalAtts,unlimdimid
 650 
      integer nDims,nGlobalAtts,unlimdimid
538 
 
 651 
 
539 
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
 652 
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
540 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 653 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
541
654
542 
      do i=1,nVars
 655 
      do i=1,nVars
543 
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
 656 
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
544 
      enddo
 657 
      enddo
545
658
546 
      end
 659 
      end