Subversion Repositories lagranto.ecmwf

Rev

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

Rev 25 Rev 27
Line 184... Line 184...
184 
        print*,'        expected -> lon / lat / lev / time'
 184 
        print*,'        expected -> lon / lat / lev / time'
185 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
 185 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
186 
        stop
 186 
        stop
187 
      endif
 187 
      endif
188 
 
 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 
 
 - 
 
215 
c     Allocate memory for reading arrays
 189 
c     Allocate memory for reading arrays
216 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
 190 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
217 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
 191 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
218 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
 192 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
219 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
 193 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
Line 269... Line 243...
269 
         do k=1,nakbktmp
 243 
         do k=1,nakbktmp
270 
            aktmp(k) = 0.01 * aktmp(k)
 244 
            aktmp(k) = 0.01 * aktmp(k)
271 
         enddo
 245 
         enddo
272 
      endif
 246 
      endif
273 
 
 247 
 
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
 248 
c     Decide whether to swap vertical levels - highest pressure at index 1
289 
      vertical_swap = 1
 249 
      vertical_swap = 1
290 
      if ( leveltype.eq.'hybrid_sigma_pressure') then
 - 
 
291 
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
 250 
      if ( (aktmp(1) + bktmp(1) * 1000.).gt.
292 
     >       (aktmp(2) + bktmp(2) * 1000.) )
 251 
     >       (aktmp(2) + bktmp(2) * 1000.) )
293 
     >  then
 252 
     >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
 253 
          vertical_swap = 0
299 
        endif
 - 
 
300 
      endif
 254 
      endif
301 
c      print*,' Vertical SWAP P -> ', vertical_swap
 - 
 
302 
 
 255 
 
303 
c     Get time information (check if time is correct)
 256 
c     Get time information (check if time is correct)
304 
      varname = 'time'
 257 
      varname = 'time'
305 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 258 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
306 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 259 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
Line 343... Line 296...
343 
         stop
 296 
         stop
344 
      endif
 297 
      endif
345 
 
 298 
 
346 
c     ---- Define output of subroutine --------------------------------
 299 
c     ---- Define output of subroutine --------------------------------
347 
 
 300 
 
- 
 
 301 
c     If not full list of vertical levels, reduce AK,BK arrays
- 
 
 302 
      if ( (leveltype.eq.'hybrid_sigma_pressure').and.
- 
 
 303 
     >     (nakbktmp.ne.vardim(3) ) )
- 
 
 304 
     >then
- 
 
 305 
         print*,' WARNING: only subset of vertical levels used...'
- 
 
 306 
         do k=1,vardim(3)
- 
 
 307 
            if ( vertical_swap.eq.1 ) then
- 
 
 308 
               aktmp(k) = aktmp( k+nakbktmp-vardim(3) )
- 
 
 309 
               bktmp(k) = bktmp( k+nakbktmp-vardim(3) )
- 
 
 310 
            endif
- 
 
 311 
         enddo
- 
 
 312 
      endif
- 
 
 313 
 
348 
c     Set the grid dimensions and constants
 314 
c     Set the grid dimensions and constants
349 
      nx      = vardim(1)
 315 
      nx      = vardim(1)
350 
      ny      = vardim(2)
 316 
      ny      = vardim(2)
351 
      nz      = vardim(3)
 317 
      nz      = vardim(3)
352 
      xmin    = lon(1)
 318 
      xmin    = lon(1)
Line 369... Line 335...
369 
 
 335 
 
370 
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
371 
      if ( fid.gt.0 ) then
 337 
      if ( fid.gt.0 ) then
372 
 
 338 
 
373 
c        Calculate layer pressures
 339 
c        Calculate layer pressures
374 
         if (leveltype.eq.'hybrid_sigma_pressure' ) then
 - 
 
375 
            do i=1,vardim(1)
 340 
         do i=1,vardim(1)
376 
              do j=1,vardim(2)
 341 
              do j=1,vardim(2)
377 
                 do k=1,vardim(3)
 342 
                 do k=1,vardim(3)
378 
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
 343 
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
379 
                 enddo
 344 
                 enddo
380 
              enddo
 345 
              enddo
381 
           enddo
 346 
         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 
 
 347 
 
392 
c        Get PS - close array on demand
 348 
c        Get PS - close array on demand
393 
         do j=1,vardim(2)
 349 
         do j=1,vardim(2)
394 
           do i=1,vardim(1)
 350 
           do i=1,vardim(1)
395 
             ps(i,j) = tmp2(i,j)
 351 
             ps(i,j) = tmp2(i,j)
Line 410... Line 366...
410 
             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)
411 
           enddo
 367 
           enddo
412 
         enddo
 368 
         enddo
413 
 
 369 
 
414 
c        Get AK,BK - vertical swap on demand
 370 
c        Get AK,BK - vertical swap on demand
415 
         if ( leveltype.eq.'hybrid_sigma_pressure' ) then
 - 
 
416 
           do k=1,vardim(3)
 371 
         do k=1,vardim(3)
417 
              if ( vertical_swap.eq.1 ) then
 372 
              if ( vertical_swap.eq.1 ) then
418 
                 ak(k) = aktmp(vardim(3)-k+1)
 373 
                 ak(k) = aktmp(vardim(3)-k+1)
419 
                 bk(k) = bktmp(vardim(3)-k+1)
 374 
                 bk(k) = bktmp(vardim(3)-k+1)
420 
              endif
 375 
              endif
421 
           enddo
 376 
         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
 - 
 
431 
           enddo
 - 
 
432 
         endif
 - 
 
433 
 
 377 
 
434 
      endif
 378 
      endif
435 
 
 379 
 
436 
 
 380 
 
437 
      return
 381 
      return
Line 548... Line 492...
548 
      ierr = NF90_INQ_DIMID(fid,varname,dimid)
 492 
      ierr = NF90_INQ_DIMID(fid,varname,dimid)
549 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 493 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
550 
      ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
 494 
      ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
551 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 495 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
552 
 
 496 
 
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 
 
 - 
 
579 
c     Allocate memory for reading arrays - depending on <closear>
 497 
c     Allocate memory for reading arrays - depending on <closear>
580 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
 498 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
581 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
 499 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
582 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
 500 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
583 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
 501 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
Line 631... Line 549...
631 
         do k=1,nakbktmp
 549 
         do k=1,nakbktmp
632 
            aktmp(k) = 0.01 * aktmp(k)
 550 
            aktmp(k) = 0.01 * aktmp(k)
633 
         enddo
 551 
         enddo
634 
      endif
 552 
      endif
635 
 
 553 
 
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
 554 
c     Decide whether to swap vertical levels
651 
      vertical_swap = 1
 555 
      vertical_swap = 1
652 
      if ( leveltype.eq.'hybrid_sigma_pressure') then
 - 
 
653 
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
 556 
      if ( (aktmp(1) + bktmp(1) * 1000.).gt.
654 
     >       (aktmp(2) + bktmp(2) * 1000.) )
 557 
     >     (aktmp(2) + bktmp(2) * 1000.) )
655 
     >  then
 558 
     >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
 559 
          vertical_swap = 0
661 
        endif
 - 
 
662 
      endif
 560 
      endif
663 
c      print*,' Vertical SWAP ',trim(fieldname),' -> ', vertical_swap
 - 
 
664 
 
 561 
 
665 
c     Read data
562 
c     Read data
666 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
 563 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
667 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 564 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
668 
      ierr = nf90_get_var(fid,varid,tmp3)
 565 
      ierr = nf90_get_var(fid,varid,tmp3)