Subversion Repositories lagranto.ecmwf

Rev

Rev 36 | Rev 46 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

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