Subversion Repositories lagranto.ecmwf

Rev

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

Rev 23 Rev 25
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 172... Line 184...
172 
        print*,'        expected -> lon / lat / lev / time'
 184 
        print*,'        expected -> lon / lat / lev / time'
173 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
 185 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
174 
        stop
 186 
        stop
175 
      endif
 187 
      endif
176 
 
 188 
 
- 
 
 189 
c     Check about the type of vertical levels
- 
 
 190 
      varname=dimname(3)
- 
 
 191 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 192 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 193 
      ierr = nf90_get_att(cdfid, varid, "standard_name", leveltype)
- 
 
 194 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 195 
      if ( (leveltype.ne.'hybrid_sigma_pressure').and.
- 
 
 196 
     >     (leveltype.ne.'air_pressure'         ) )
- 
 
 197 
     >then
- 
 
 198 
         print*,' ERROR: input netCDF data must be on hybrid-sigma'
- 
 
 199 
         print*,'        or air pressure levels!',trim(leveltype)
- 
 
 200 
         stop
- 
 
 201 
      endif
- 
 
 202 
 
- 
 
 203 
c     Check that vardim(3)==#AK,BK for hybrid-sigmal levels
- 
 
 204 
      if ( (leveltype.eq.'hybrid_sigma_pressure').and.
- 
 
 205 
     >     (dimname(3).eq.'lev') )
- 
 
 206 
     >then
- 
 
 207 
        if ( nakbktmp.ne.vardim(3) ) then
- 
 
 208 
           print*,' ERROR: for hybrid-sigma pressure levels, #AK,BK'
- 
 
 209 
           print*,'        must agree with number of vertical levels'
- 
 
 210 
           print*,'        ',vardim(3),nakbktmp
- 
 
 211 
           stop
- 
 
 212 
        endif
- 
 
 213 
      endif
- 
 
 214 
 
177 
c     Allocate memory for reading arrays
 215 
c     Allocate memory for reading arrays
178 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
 216 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
179 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
 217 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
180 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
 218 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
181 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
 219 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
182 
      allocate(lon(vardim(1)),stat=stat)
 220 
      allocate(lon(vardim(1)),stat=stat)
183 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
221 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
184 
      allocate(lat(vardim(2)),stat=stat)
 222 
      allocate(lat(vardim(2)),stat=stat)
185 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
223 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
- 
 
 224 
      allocate(lev(vardim(3)),stat=stat)
- 
 
 225 
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
186 
      allocate(times(vardim(4)),stat=stat)
 226 
      allocate(times(vardim(4)),stat=stat)
187 
      if (stat.ne.0) print*,'*** error allocating array times   ***'
 227 
      if (stat.ne.0) print*,'*** error allocating array times   ***'
188 
      allocate(aktmp(vardim(3)),stat=stat)
 228 
      allocate(aktmp(nakbktmp),stat=stat)
189 
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
 229 
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
190 
      allocate(bktmp(vardim(3)),stat=stat)
 230 
      allocate(bktmp(nakbktmp),stat=stat)
191 
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
 231 
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
192 
 
 232 
 
193 
c     Get domain longitudes and latitudes
 233 
c     Get domain longitudes, latitudes and levels
194 
      varname = dimname(1)
 234 
      varname = dimname(1)
195 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 235 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
196 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 236 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
197 
      ierr = nf90_get_var(cdfid,varid,lon)
 237 
      ierr = nf90_get_var(cdfid,varid,lon)
198 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 238 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
199 
      varname = dimname(2)
 239 
      varname = dimname(2)
200 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 240 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
201 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 241 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
202 
      ierr = nf90_get_var(cdfid,varid,lat)
 242 
      ierr = nf90_get_var(cdfid,varid,lat)
203 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 243 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 244 
      varname = dimname(3)
- 
 
 245 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 246 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 247 
      ierr = nf90_get_var(cdfid,varid,lev)
- 
 
 248 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
204
249
205 
c     Get ak and bk
 250 
c     Get ak and bk
206 
      varname='hyam'
 251 
      varname='hyam'
207 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 252 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
208 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 253 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
Line 219... Line 264...
219 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 264 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
220 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 265 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
221 
      ierr = nf90_get_att(cdfid, varid, "units", units)
 266 
      ierr = nf90_get_att(cdfid, varid, "units", units)
222 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 267 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
223 
      if ( units.eq.'Pa' ) then
 268 
      if ( units.eq.'Pa' ) then
224 
         do k=1,vardim(3)
 269 
         do k=1,nakbktmp
225 
            aktmp(k) = 0.01 * aktmp(k)
 270 
            aktmp(k) = 0.01 * aktmp(k)
226 
         enddo
 271 
         enddo
227 
      endif
 272 
      endif
228 
 
 273 
 
- 
 
 274 
c     Check that unit of lev is in hPa - if necessary correct it
- 
 
 275 
      if ( leveltype.eq.'air_pressure' ) then
- 
 
 276 
         varname='lev'
- 
 
 277 
         ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 278 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 279 
         ierr = nf90_get_att(cdfid, varid, "units", units)
- 
 
 280 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 281 
         if ( units.eq.'Pa' ) then
- 
 
 282 
            do k=1,vardim(3)
- 
 
 283 
               lev(k) = 0.01 * lev(k)
- 
 
 284 
            enddo
- 
 
 285 
         endif
- 
 
 286 
      endif
- 
 
 287 
 
- 
 
 288 
c     Decide whether to swap vertical levels
- 
 
 289 
      vertical_swap = 1
- 
 
 290 
      if ( leveltype.eq.'hybrid_sigma_pressure') then
- 
 
 291 
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
- 
 
 292 
     >       (aktmp(2) + bktmp(2) * 1000.) )
- 
 
 293 
     >  then
- 
 
 294 
          vertical_swap = 0
- 
 
 295 
        endif
- 
 
 296 
      elseif ( leveltype.eq.'air_pressure') then
- 
 
 297 
        if ( lev(1).gt.lev(2) ) then
- 
 
 298 
          vertical_swap = 0
- 
 
 299 
        endif
- 
 
 300 
      endif
- 
 
 301 
c      print*,' Vertical SWAP P -> ', vertical_swap
- 
 
 302 
 
229 
c     Get time information (check if time is correct)
 303 
c     Get time information (check if time is correct)
230 
      varname = 'time'
 304 
      varname = 'time'
231 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 305 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
232 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 306 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
233 
      ierr = nf90_get_var(cdfid,varid,times)
 307 
      ierr = nf90_get_var(cdfid,varid,times)
Line 253... Line 327...
253 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 327 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
254 
      ierr = nf90_get_var(cdfid,varid,tmp2)
 328 
      ierr = nf90_get_var(cdfid,varid,tmp2)
255 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 329 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
256
330
257 
c     Check that surface pressure is in hPa
 331 
c     Check that surface pressure is in hPa
258 
      maxps = -1.e39
 332 
      maxps = -1.e19
259 
      minps =  1.e39
 333 
      minps =  1.e19
260 
      do i=1,vardim(1)
 334 
      do i=1,vardim(1)
261 
        do j=1,vardim(2)
 335 
        do j=1,vardim(2)
262 
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
 336 
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
263 
             if (tmp2(i,j).lt.minps) minps = tmp2(i,j)
 337 
             if (tmp2(i,j).lt.minps) minps = tmp2(i,j)
264 
        enddo
 338 
        enddo
Line 267... Line 341...
267 
         print*,' ERROR: surface pressre PS must be in hPa'
 341 
         print*,' ERROR: surface pressre PS must be in hPa'
268 
         print*,'       ',maxps,minps
 342 
         print*,'       ',maxps,minps
269 
         stop
 343 
         stop
270 
      endif
 344 
      endif
271 
 
 345 
 
272 
c     Calculate layer and level pressures
 - 
 
273 
      do i=1,vardim(1)
 - 
 
274 
         do j=1,vardim(2)
 - 
 
275 
               do k=1,vardim(3)
 - 
 
276 
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
 346 
c     ---- Define output of subroutine --------------------------------
277 
               enddo
 - 
 
278 
         enddo
 - 
 
279 
      enddo
 - 
 
280 
 
 347 
 
281 
c     Set the grid dimensions and constants
 348 
c     Set the grid dimensions and constants
282 
      nx      = vardim(1)
 349 
      nx      = vardim(1)
283 
      ny      = vardim(2)
 350 
      ny      = vardim(2)
284 
      nz      = vardim(3)
 351 
      nz      = vardim(3)
Line 301... Line 368...
301 
      endif
 368 
      endif
302 
 
 369 
 
303 
c     Save the output arrays (if fid>0) - close arrays on request
 370 
c     Save the output arrays (if fid>0) - close arrays on request
304 
      if ( fid.gt.0 ) then
 371 
      if ( fid.gt.0 ) then
305 
 
 372 
 
- 
 
 373 
c        Calculate layer pressures
- 
 
 374 
         if (leveltype.eq.'hybrid_sigma_pressure' ) then
- 
 
 375 
            do i=1,vardim(1)
- 
 
 376 
              do j=1,vardim(2)
- 
 
 377 
                 do k=1,vardim(3)
- 
 
 378 
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
- 
 
 379 
                 enddo
- 
 
 380 
              enddo
- 
 
 381 
           enddo
- 
 
 382 
         elseif (leveltype.eq.'air_pressure' ) then
- 
 
 383 
           do i=1,vardim(1)
- 
 
 384 
              do j=1,vardim(2)
- 
 
 385 
                 do k=1,vardim(3)
- 
 
 386 
                  tmp3(i,j,k)=lev(k)
- 
 
 387 
                 enddo
- 
 
 388 
              enddo
- 
 
 389 
           enddo
- 
 
 390 
         endif
- 
 
 391 
 
- 
 
 392 
c        Get PS - close array on demand
306 
         do j=1,vardim(2)
 393 
         do j=1,vardim(2)
307 
           do i=1,vardim(1)
 394 
           do i=1,vardim(1)
308 
             ps(i,j) = tmp2(i,j)
 395 
             ps(i,j) = tmp2(i,j)
309 
           enddo
 396 
           enddo
310 
           if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)
 397 
           if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)
311 
         enddo
 398 
         enddo
312 
 
 399 
 
- 
 
 400 
c        Get P3 - close array on demand + vertical swap
313 
         do j=1,vardim(2)
 401 
         do j=1,vardim(2)
314 
           do k=1,vardim(3)
 402 
           do k=1,vardim(3)
315 
             do i=1,vardim(1)
 403 
             do i=1,vardim(1)
- 
 
 404 
               if ( vertical_swap.eq.1 ) then
316 
               p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)
 405 
                  p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)
- 
 
 406 
               else
- 
 
 407 
                  p3(i,j,k) = tmp3(i,j,k)
- 
 
 408 
               endif
317 
             enddo
 409 
             enddo
318 
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
 410 
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
319 
           enddo
 411 
           enddo
320 
         enddo
 412 
         enddo
321 
 
 413 
 
- 
 
 414 
c        Get AK,BK - vertical swap on demand
- 
 
 415 
         if ( leveltype.eq.'hybrid_sigma_pressure' ) then
322 
         do k=1,vardim(3)
 416 
           do k=1,vardim(3)
- 
 
 417 
              if ( vertical_swap.eq.1 ) then
323 
            ak(k) = aktmp(vardim(3)-k+1)
 418 
                 ak(k) = aktmp(vardim(3)-k+1)
324 
            bk(k) = bktmp(vardim(3)-k+1)
 419 
                 bk(k) = bktmp(vardim(3)-k+1)
- 
 
 420 
              endif
- 
 
 421 
           enddo
- 
 
 422 
         elseif (leveltype.eq.'air_pressure' ) then
- 
 
 423 
           do k=1,vardim(3)
- 
 
 424 
              if ( vertical_swap.eq.1 ) then
- 
 
 425 
                 ak(k) = lev(vardim(3)-k+1)
- 
 
 426 
                 bk(k) = 0.
- 
 
 427 
              else
- 
 
 428 
                ak(k) = lev(k)
- 
 
 429 
                bk(k) = 0.
- 
 
 430 
              endif
325 
         enddo
 431 
           enddo
- 
 
 432 
         endif
326 
 
 433 
 
327 
      endif
 434 
      endif
328 
 
 435 
 
- 
 
 436 
 
329 
      return
 437 
      return
330
438
331 
      end
 439 
      end
332 
 
 440 
 
333 
c     ------------------------------------------------------------
 441 
c     ------------------------------------------------------------
Line 382... Line 490...
382 
      real         ps(nx,ny)
 490 
      real         ps(nx,ny)
383 
      integer      dimids (nf90_max_var_dims)
 491 
      integer      dimids (nf90_max_var_dims)
384 
      character*80 dimname(nf90_max_var_dims)
 492 
      character*80 dimname(nf90_max_var_dims)
385 
      integer      varid
 493 
      integer      varid
386 
      integer      cdfid
 494 
      integer      cdfid
387 
      real,allocatable, dimension (:)     :: lon,lat
 495 
      real,allocatable, dimension (:)     :: lon,lat,lev
388 
      real,allocatable, dimension (:,:)   :: tmp2
 496 
      real,allocatable, dimension (:,:)   :: tmp2
389 
      real,allocatable, dimension (:,:,:) :: tmp3
 497 
      real,allocatable, dimension (:,:,:) :: tmp3
- 
 
 498 
      real,allocatable, dimension (:)     :: aktmp,bktmp
- 
 
 499 
      character*80  leveltype
- 
 
 500 
      integer       vertical_swap
- 
 
 501 
      character*80  units
- 
 
 502 
      integer       nakbktmp
- 
 
 503 
      integer       dimid
390 
 
 504 
 
391 
c     Auxiliary variables
 505 
c     Auxiliary variables
392 
      integer      isok
 506 
      integer      isok
393 
      integer      i,j,k
 507 
      integer      i,j,k
394 
      integer      nz1
 508 
      integer      nz1
Line 427... Line 541...
427 
        print*,'        expected -> lon / lat / lev / time'
 541 
        print*,'        expected -> lon / lat / lev / time'
428 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
 542 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
429 
        stop
 543 
        stop
430 
      endif
 544 
      endif
431 
 
 545 
 
- 
 
 546 
c     Get dimension of AK,BK
- 
 
 547 
      varname = 'nhym'
- 
 
 548 
      ierr = NF90_INQ_DIMID(fid,varname,dimid)
- 
 
 549 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 550 
      ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
- 
 
 551 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 552 
 
- 
 
 553 
c     Check about the type of vertical levels
- 
 
 554 
      varname=dimname(3)
- 
 
 555 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 556 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 557 
      ierr = nf90_get_att(fid, varid, "standard_name", leveltype)
- 
 
 558 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 559 
      if ( (leveltype.ne.'hybrid_sigma_pressure').and.
- 
 
 560 
     >     (leveltype.ne.'air_pressure'         ) )
- 
 
 561 
     >then
- 
 
 562 
         print*,' ERROR: input netCDF data must be on hybrid-sigma'
- 
 
 563 
         print*,'        or air pressure levels!',trim(leveltype)
- 
 
 564 
         stop
- 
 
 565 
      endif
- 
 
 566 
 
- 
 
 567 
c     Check that vardim(3)==#AK,BK for hybrid-sigmal levels
- 
 
 568 
      if ( (leveltype.eq.'hybrid_sigma_pressure').and.
- 
 
 569 
     >     (dimname(3).eq.'lev') )
- 
 
 570 
     >then
- 
 
 571 
        if ( nakbktmp.ne.vardim(3) ) then
- 
 
 572 
           print*,' ERROR: for hybrid-sigma pressure levels, #AK,BK'
- 
 
 573 
           print*,'        must agree with number of vertical levels'
- 
 
 574 
           print*,'        ',vardim(3),nakbktmp
- 
 
 575 
           stop
- 
 
 576 
        endif
- 
 
 577 
      endif
- 
 
 578 
 
432 
c     Allocate memory for reading arrays - depending on <closear>
 579 
c     Allocate memory for reading arrays - depending on <closear>
433 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
 580 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
434 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
 581 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
435 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
 582 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
436 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
 583 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
437 
      allocate(lon(vardim(1)),stat=stat)
 584 
      allocate(lon(vardim(1)),stat=stat)
438 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
 585 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
439 
      allocate(lat(vardim(2)),stat=stat)
 586 
      allocate(lat(vardim(2)),stat=stat)
440 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
 587 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
- 
 
 588 
      allocate(lev(vardim(3)),stat=stat)
- 
 
 589 
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
- 
 
 590 
      allocate(aktmp(nakbktmp),stat=stat)
- 
 
 591 
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
- 
 
 592 
      allocate(bktmp(nakbktmp),stat=stat)
- 
 
 593 
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
441 
 
 594 
 
442 
c     Get domain boundaries
 595 
c     Get domain boundaries - longitude, latitude, levels
443 
      varname = dimname(1)
 596 
      varname = dimname(1)
444 
      ierr = NF90_INQ_VARID(fid,varname,varid)
 597 
      ierr = NF90_INQ_VARID(fid,varname,varid)
445 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 598 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
446 
      ierr = nf90_get_var(fid,varid,lon)
 599 
      ierr = nf90_get_var(fid,varid,lon)
447 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 600 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
448 
      varname = dimname(2)
 601 
      varname = dimname(2)
449 
      ierr = NF90_INQ_VARID(fid,varname,varid)
 602 
      ierr = NF90_INQ_VARID(fid,varname,varid)
450 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 603 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
451 
      ierr = nf90_get_var(fid,varid,lat)
 604 
      ierr = nf90_get_var(fid,varid,lat)
452 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 605 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 606 
      varname = dimname(3)
- 
 
 607 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 608 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 609 
      ierr = nf90_get_var(fid,varid,lev)
- 
 
 610 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 611 
 
- 
 
 612 
c     Get ak and bk
- 
 
 613 
      varname='hyam'
- 
 
 614 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 615 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 616 
      ierr = nf90_get_var(fid,varid,aktmp)
- 
 
 617 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 618 
      varname='hybm'
- 
 
 619 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 620 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 621 
      ierr = nf90_get_var(fid,varid,bktmp)
- 
 
 622 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 623 
 
- 
 
 624 
c     Check that unit of ak is in hPa - if necessary correct it
- 
 
 625 
      varname='hyam'
- 
 
 626 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 627 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 628 
      ierr = nf90_get_att(fid, varid, "units", units)
- 
 
 629 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 630 
      if ( units.eq.'Pa' ) then
- 
 
 631 
         do k=1,nakbktmp
- 
 
 632 
            aktmp(k) = 0.01 * aktmp(k)
- 
 
 633 
         enddo
- 
 
 634 
      endif
- 
 
 635 
 
- 
 
 636 
c     Check that unit of lev is in hPa - if necessary correct it
- 
 
 637 
      if ( leveltype.eq.'air_pressure' ) then
- 
 
 638 
         varname='lev'
- 
 
 639 
         ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 640 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 641 
         ierr = nf90_get_att(fid, varid, "units", units)
- 
 
 642 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 643 
         if ( units.eq.'Pa' ) then
- 
 
 644 
            do k=1,vardim(3)
- 
 
 645 
               lev(k) = 0.01 * lev(k)
- 
 
 646 
            enddo
- 
 
 647 
         endif
- 
 
 648 
      endif
- 
 
 649 
 
- 
 
 650 
c     Decide whether to swap vertical levels
- 
 
 651 
      vertical_swap = 1
- 
 
 652 
      if ( leveltype.eq.'hybrid_sigma_pressure') then
- 
 
 653 
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
- 
 
 654 
     >       (aktmp(2) + bktmp(2) * 1000.) )
- 
 
 655 
     >  then
- 
 
 656 
          vertical_swap = 0
- 
 
 657 
        endif
- 
 
 658 
      elseif ( leveltype.eq.'air_pressure') then
- 
 
 659 
        if ( lev(1).gt.lev(2) ) then
- 
 
 660 
          vertical_swap = 0
- 
 
 661 
        endif
- 
 
 662 
      endif
- 
 
 663 
c      print*,' Vertical SWAP ',trim(fieldname),' -> ', vertical_swap
453 
 
 664 
 
454 
c     Read data
665 
c     Read data
455 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
 666 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
456 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 667 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
457 
      ierr = nf90_get_var(fid,varid,tmp3)
 668 
      ierr = nf90_get_var(fid,varid,tmp3)
Line 466... Line 677...
466 
               enddo
 677 
               enddo
467 
            enddo
 678 
            enddo
468 
         enddo
 679 
         enddo
469 
      endif
 680 
      endif
470 
 
 681 
 
471 
c     Save the ouput array - close on request
 682 
c     Decide whether to close arrays
472 
      delta = varmax(1)-varmin(1)-360.
 683 
      delta = varmax(1)-varmin(1)-360.
473 
      if (abs(delta+dx).lt.eps) then
 684 
      if (abs(delta+dx).lt.eps) then
474 
          closear = 1
 685 
          closear = 1
475 
      else
 686 
      else
476 
          closear = 0
 687 
          closear = 0
477 
      endif
 688 
      endif
478 
 
 689 
 
- 
 
 690 
c     Save output array - close array and swap on demand
479 
      do j=1,vardim(2)
 691 
      do j=1,vardim(2)
480 
        do k=1,vardim(3)
 692 
        do k=1,vardim(3)
481 
          do i=1,vardim(1)
 693 
          do i=1,vardim(1)
- 
 
 694 
             if ( vertical_swap.eq.1 ) then
482 
             field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
 695 
                 field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
- 
 
 696 
             else
- 
 
 697 
                 field(i,j,k) = tmp3(i,j,k)
- 
 
 698 
             endif
483 
          enddo
 699 
          enddo
484 
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
 700 
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
485 
        enddo
 701 
        enddo
486 
      enddo
 702 
      enddo
487
703