Subversion Repositories lagranto.ecmwf

Rev

Rev 25 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 25 Rev 27
Line 100... Line 100...
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
Line 123... Line 127...
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
Line 160... Line 166...
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' ) )
Line 181... Line 193...
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)
Line 219... Line 238...
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)
Line 253... Line 280...
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
Line 267... Line 294...
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)
Line 301... Line 334...
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     ------------------------------------------------------------
Line 382... Line 434...
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
Line 427... Line 485...
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)
Line 466... Line 574...
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