Subversion Repositories lagranto.ecmwf

Rev

Rev 27 | Rev 31 | Go to most recent revision | Show entire file | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 27 Rev 29
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 
 
189 
c     Allocate memory for reading arrays
 203 
c     Allocate memory for reading arrays
190 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
 204 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
191 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
 205 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
192 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
 206 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
193 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
 207 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
Line 243... Line 257...
243 
         do k=1,nakbktmp
 257 
         do k=1,nakbktmp
244 
            aktmp(k) = 0.01 * aktmp(k)
 258 
            aktmp(k) = 0.01 * aktmp(k)
245 
         enddo
 259 
         enddo
246 
      endif
 260 
      endif
247 
 
 261 
 
- 
 
 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
- 
 
 274 
      endif
- 
 
 275 
 
248 
c     Decide whether to swap vertical levels - highest pressure at index 1
 276 
c     Decide whether to swap vertical levels - highest pressure at index 1
249 
      vertical_swap = 1
 277 
      vertical_swap = 1
- 
 
 278 
      if ( leveltype.eq.'hybrid_sigma_pressure') then
250 
      if ( (aktmp(1) + bktmp(1) * 1000.).gt.
 279 
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
251 
     >       (aktmp(2) + bktmp(2) * 1000.) )
 280 
     >       (aktmp(2) + bktmp(2) * 1000.) )
252 
     >then
 281 
     >  then
253 
          vertical_swap = 0
 282 
          vertical_swap = 0
254 
      endif
 283 
        endif
- 
 
 284 
      elseif ( leveltype.eq.'air_pressure') then
- 
 
 285 
        if ( lev(1).gt.lev(2) ) then
- 
 
 286 
          vertical_swap = 0
- 
 
 287 
        endif
- 
 
 288 
      endif
- 
 
 289 
c      print*,' Vertical SWAP P -> ', vertical_swap
255 
 
 290 
 
256 
c     Get time information (check if time is correct)
 291 
c     Get time information (check if time is correct)
257 
      varname = 'time'
 292 
      varname = 'time'
258 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 293 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
259 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 294 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
Line 335... Line 370...
335 
 
 370 
 
336 
c     Save the output arrays (if fid>0) - close arrays on request
 371 
c     Save the output arrays (if fid>0) - close arrays on request
337 
      if ( fid.gt.0 ) then
 372 
      if ( fid.gt.0 ) then
338 
 
 373 
 
339 
c        Calculate layer pressures
 374 
c        Calculate layer pressures
- 
 
 375 
         if (leveltype.eq.'hybrid_sigma_pressure' ) then
340 
         do i=1,vardim(1)
 376 
            do i=1,vardim(1)
341 
              do j=1,vardim(2)
 377 
              do j=1,vardim(2)
342 
                 do k=1,vardim(3)
 378 
                 do k=1,vardim(3)
343 
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
 379 
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
344 
                 enddo
 380 
                 enddo
345 
              enddo
 381 
              enddo
346 
         enddo
 382 
           enddo
- 
 
 383 
         elseif (leveltype.eq.'air_pressure' ) then
- 
 
 384 
           do i=1,vardim(1)
- 
 
 385 
              do j=1,vardim(2)
- 
 
 386 
                 do k=1,vardim(3)
- 
 
 387 
                  tmp3(i,j,k)=lev(k)
- 
 
 388 
                 enddo
- 
 
 389 
              enddo
- 
 
 390 
           enddo
- 
 
 391 
         endif
347 
 
 392 
 
348 
c        Get PS - close array on demand
 393 
c        Get PS - close array on demand
349 
         do j=1,vardim(2)
 394 
         do j=1,vardim(2)
350 
           do i=1,vardim(1)
 395 
           do i=1,vardim(1)
351 
             ps(i,j) = tmp2(i,j)
 396 
             ps(i,j) = tmp2(i,j)
Line 366... Line 411...
366 
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
 411 
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
367 
           enddo
 412 
           enddo
368 
         enddo
 413 
         enddo
369 
 
 414 
 
370 
c        Get AK,BK - vertical swap on demand
 415 
c        Get AK,BK - vertical swap on demand
- 
 
 416 
         if ( leveltype.eq.'hybrid_sigma_pressure' ) then
371 
         do k=1,vardim(3)
 417 
           do k=1,vardim(3)
372 
              if ( vertical_swap.eq.1 ) then
 418 
              if ( vertical_swap.eq.1 ) then
373 
                 ak(k) = aktmp(vardim(3)-k+1)
 419 
                 ak(k) = aktmp(vardim(3)-k+1)
374 
                 bk(k) = bktmp(vardim(3)-k+1)
 420 
                 bk(k) = bktmp(vardim(3)-k+1)
375 
              endif
 421 
              endif
376 
         enddo
 422 
           enddo
- 
 
 423 
         elseif (leveltype.eq.'air_pressure' ) then
- 
 
 424 
           do k=1,vardim(3)
- 
 
 425 
              if ( vertical_swap.eq.1 ) then
- 
 
 426 
                 ak(k) = lev(vardim(3)-k+1)
- 
 
 427 
                 bk(k) = 0.
- 
 
 428 
              else
- 
 
 429 
                ak(k) = lev(k)
- 
 
 430 
                bk(k) = 0.
- 
 
 431 
              endif
- 
 
 432 
           enddo
- 
 
 433 
         endif
377 
 
 434 
 
378 
      endif
 435 
      endif
379 
 
 436 
 
380 
 
 437 
 
381 
      return
 438 
      return
Line 453... Line 510...
453 
      real         rtime
 510 
      real         rtime
454 
      integer      closear
 511 
      integer      closear
455 
      integer      stat
 512 
      integer      stat
456 
      real         delta
 513 
      real         delta
457 
 
 514 
 
- 
 
 515 
c     Init mdv
- 
 
 516 
      mdv = -999.
- 
 
 517 
 
458 
c     Get the number of dimensions -> ndim
 518 
c     Get the number of dimensions -> ndim
459 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
 519 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
460 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 520 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
461 
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
 521 
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
462 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 522 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
Line 492... Line 552...
492 
      ierr = NF90_INQ_DIMID(fid,varname,dimid)
 552 
      ierr = NF90_INQ_DIMID(fid,varname,dimid)
493 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 553 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
494 
      ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
 554 
      ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
495 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 555 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
496 
 
 556 
 
- 
 
 557 
c     Check about the type of vertical levels
- 
 
 558 
      varname=dimname(3)
- 
 
 559 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 560 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 561 
      ierr = nf90_get_att(fid, varid, "standard_name", leveltype)
- 
 
 562 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 563 
      if ( (leveltype.ne.'hybrid_sigma_pressure').and.
- 
 
 564 
     >     (leveltype.ne.'air_pressure'         ) )
- 
 
 565 
     >then
- 
 
 566 
         print*,' ERROR: input netCDF data must be on hybrid-sigma'
- 
 
 567 
         print*,'        or air pressure levels!',trim(leveltype)
- 
 
 568 
         stop
- 
 
 569 
      endif
- 
 
 570 
 
497 
c     Allocate memory for reading arrays - depending on <closear>
 571 
c     Allocate memory for reading arrays - depending on <closear>
498 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
 572 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
499 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
 573 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
500 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
 574 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
501 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
 575 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
Line 549... Line 623...
549 
         do k=1,nakbktmp
 623 
         do k=1,nakbktmp
550 
            aktmp(k) = 0.01 * aktmp(k)
 624 
            aktmp(k) = 0.01 * aktmp(k)
551 
         enddo
 625 
         enddo
552 
      endif
 626 
      endif
553 
 
 627 
 
- 
 
 628 
c     Check that unit of lev is in hPa - if necessary correct it
- 
 
 629 
      if ( leveltype.eq.'air_pressure' ) then
- 
 
 630 
         varname='lev'
- 
 
 631 
         ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 632 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 633 
         ierr = nf90_get_att(fid, varid, "units", units)
- 
 
 634 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 635 
         if ( units.eq.'Pa' ) then
- 
 
 636 
            do k=1,vardim(3)
- 
 
 637 
               lev(k) = 0.01 * lev(k)
- 
 
 638 
            enddo
- 
 
 639 
         endif
- 
 
 640 
      endif
- 
 
 641 
 
554 
c     Decide whether to swap vertical levels
 642 
c     Decide whether to swap vertical levels
555 
      vertical_swap = 1
 643 
      vertical_swap = 1
- 
 
 644 
      if ( leveltype.eq.'hybrid_sigma_pressure') then
556 
      if ( (aktmp(1) + bktmp(1) * 1000.).gt.
 645 
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
557 
     >     (aktmp(2) + bktmp(2) * 1000.) )
 646 
     >       (aktmp(2) + bktmp(2) * 1000.) )
558 
     >then
 647 
     >  then
559 
          vertical_swap = 0
 648 
          vertical_swap = 0
560 
      endif
 649 
        endif
- 
 
 650 
      elseif ( leveltype.eq.'air_pressure') then
- 
 
 651 
        if ( lev(1).gt.lev(2) ) then
- 
 
 652 
          vertical_swap = 0
- 
 
 653 
        endif
- 
 
 654 
      endif
561 
 
 655 
 
562 
c     Read data
656 
c     Read data
563 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
 657 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
564 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 658 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
565 
      ierr = nf90_get_var(fid,varid,tmp3)
 659 
      ierr = nf90_get_var(fid,varid,tmp3)