Subversion Repositories lagranto.ecmwf

Rev

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

Rev 36 Rev 46
Line 146... Line 146...
146 
c     Get number of dimensions of variable -> ndim
 146 
c     Get number of dimensions of variable -> ndim
147 
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
 147 
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
148 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 148 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
149 
      ierr = nf90_inquire_variable(cdfid, varid, ndims  = ndim)
 149 
      ierr = nf90_inquire_variable(cdfid, varid, ndims  = ndim)
150 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 150 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
151 
      if ( ndim.ne.4 ) then
 151 
c      if ( ndim.ne.4 ) then
152 
          print*,' ERROR: netCDF variables need to be 4D'
 152 
c          print*,' ERROR: netCDF variables need to be 4D'
153 
          print*,'      ',trim(fieldname)
 153 
c          print*,'      ',trim(fieldname)
154 
          stop
 154 
c          stop
155 
      endif
 155 
c      endif
156 
 
 156 
 
157 
c     Get dimensions -> vardim(1:ndim),dimname(1:ndim)
 157 
c     Get dimensions -> vardim(1:ndim),dimname(1:ndim)
158 
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
 158 
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
159 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 159 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
160 
      ierr = nf90_inquire_variable(cdfid, varid,
160 
      ierr = nf90_inquire_variable(cdfid, varid,
Line 166... Line 166...
166 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 166 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
167 
           ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))
 167 
           ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))
168 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 168 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
169 
      enddo
 169 
      enddo
170 
 
 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 
 
 - 
 
177 
c     Check whether the list of dimensions is OK
 171 
c     Check whether the list of dimensions is OK
- 
 
 172 
      if ( ndim.eq.4 ) then
178 
      if ( ( dimname(1).ne.'lon'  ).or.
 173 
        if ( ( dimname(1).ne.'lon'  ).or.
179 
     >     ( dimname(2).ne.'lat'  ).or.
174 
     >     ( dimname(2).ne.'lat'  ).or.
180 
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
 175 
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
181 
     >     ( dimname(4).ne.'time' ) )
 176 
     >     ( dimname(4).ne.'time' ) )
182 
     >then
 177 
     >  then
183 
        print*,' ERROR: the dimensions of the variable are not correct'
 178 
         print*,' ERROR: the dimensions of the variable are not correct'
184 
        print*,'        expected -> lon / lat / lev / time'
 179 
         print*,'        expected -> lon / lat / lev / time'
185 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
 180 
         print*, ( trim(dimname(i))//' / ',i=1,ndim )
186 
        stop
 181 
         stop
187 
      endif
 182 
        endif
- 
 
 183 
      else
- 
 
 184 
        if ( ( dimname(1).ne.'lon'  ).or.
- 
 
 185 
     >       ( dimname(2).ne.'lat'  ).or.
- 
 
 186 
     >       ( dimname(3).ne.'time' ) )
- 
 
 187 
     >  then
- 
 
 188 
         print*,' ERROR: the dimensions of the variable are not correct'
- 
 
 189 
         print*,'        expected -> lon / lat / time'
- 
 
 190 
         print*, ( trim(dimname(i))//' / ',i=1,ndim )
- 
 
 191 
         stop
- 
 
 192 
        endif
- 
 
 193 
      endif
188 
 
 194 
 
189 
c     Check about the type of vertical levels
 195 
c     Check about the type of vertical levels
- 
 
 196 
      if ( ndim.eq.4 ) then
- 
 
 197 
 
- 
 
 198 
c         Get dimension of AK,BK
- 
 
 199 
          varname = 'nhym'
- 
 
 200 
          ierr = NF90_INQ_DIMID(cdfid,varname,dimid)
- 
 
 201 
          ierr = nf90_inquire_dimension(cdfid, dimid,len=nakbktmp)
- 
 
 202 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 203 
 
- 
 
 204 
 
190 
      varname=dimname(3)
 205 
          varname=dimname(3)
191 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 206 
          ierr = NF90_INQ_VARID(cdfid,varname,varid)
192 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 207 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
193 
      ierr = nf90_get_att(cdfid, varid, "standard_name", leveltype)
 208 
          ierr = nf90_get_att(cdfid, varid, "standard_name", leveltype)
194 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 209 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
Line 197... Line 212...
197 
     >then
 212 
     >    then
198 
         print*,' ERROR: input netCDF data must be on hybrid-sigma'
 213 
             print*,' ERROR: input netCDF data must be on hybrid-sigma'
199 
         print*,'        or air pressure levels!',trim(leveltype)
 214 
             print*,'        or air pressure levels!',trim(leveltype)
200 
         stop
 215 
             stop
201 
      endif
 216 
          endif
- 
 
 217 
      else
- 
 
 218 
          leveltype='air_pressure'
- 
 
 219 
      endif
202 
 
 220 
 
203 
c     Allocate memory for reading arrays
 221 
c     Allocate memory for reading arrays
204 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
 222 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
205 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
 223 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
206 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
 224 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
Line 211... Line 229...
211 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
229 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
212 
      allocate(lev(vardim(3)),stat=stat)
 230 
      allocate(lev(vardim(3)),stat=stat)
213 
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
 231 
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
214 
      allocate(times(vardim(4)),stat=stat)
 232 
      allocate(times(vardim(4)),stat=stat)
215 
      if (stat.ne.0) print*,'*** error allocating array times   ***'
 233 
      if (stat.ne.0) print*,'*** error allocating array times   ***'
- 
 
 234 
      if ( ndim.eq.4 ) then
216 
      allocate(aktmp(nakbktmp),stat=stat)
 235 
        allocate(aktmp(nakbktmp),stat=stat)
217 
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
 236 
        if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
218 
      allocate(bktmp(nakbktmp),stat=stat)
 237 
        allocate(bktmp(nakbktmp),stat=stat)
219 
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
 238 
        if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
- 
 
 239 
      endif
220 
 
 240 
 
221 
c     Get domain longitudes, latitudes and levels
 241 
c     Get domain longitudes, latitudes and levels
222 
      varname = dimname(1)
 242 
      varname = dimname(1)
223 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 243 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
224 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 244 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
Line 234... Line 254...
234 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 254 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
235 
      ierr = nf90_get_var(cdfid,varid,lev)
 255 
      ierr = nf90_get_var(cdfid,varid,lev)
236 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 256 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
237
257
238 
c     Get ak and bk
 258 
c     Get ak and bk
- 
 
 259 
      if ( ndim.eq.4 ) then
- 
 
 260 
 
- 
 
 261 
c       Read ak, bk
239 
      varname='hyam'
 262 
        varname='hyam'
240 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
 263 
        ierr = NF90_INQ_VARID(cdfid,varname,varid)
241 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 264 
        IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
242 
      ierr = nf90_get_var(cdfid,varid,aktmp)
 265 
        ierr = nf90_get_var(cdfid,varid,aktmp)
243 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 266 
        IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
Line 256... Line 279...
256 
      if ( units.eq.'Pa' ) then
 279 
        if ( units.eq.'Pa' ) then
257 
         do k=1,nakbktmp
 280 
           do k=1,nakbktmp
258 
            aktmp(k) = 0.01 * aktmp(k)
 281 
              aktmp(k) = 0.01 * aktmp(k)
259 
         enddo
 282 
           enddo
260 
      endif
 283 
        endif
- 
 
 284 
      endif
261 
 
 285 
 
262 
c     Check that unit of lev is in hPa - if necessary correct it
 286 
c     Check that unit of lev is in hPa - if necessary correct it
263 
      if ( leveltype.eq.'air_pressure' ) then
 287 
      if ( leveltype.eq.'air_pressure' .and. vardim(3) .gt. 1) then
264 
         varname='lev'
 288 
         varname='lev'
265 
         ierr = NF90_INQ_VARID(cdfid,varname,varid)
 289 
         ierr = NF90_INQ_VARID(cdfid,varname,varid)
266 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 290 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
267 
         ierr = nf90_get_att(cdfid, varid, "units", units)
 291 
         ierr = nf90_get_att(cdfid, varid, "units", units)
268 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 292 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
Line 315... Line 339...
315 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 339 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
316 
      ierr = nf90_get_var(cdfid,varid,tmp2)
 340 
      ierr = nf90_get_var(cdfid,varid,tmp2)
317 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 341 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
318
342
319 
c     Check that surface pressure is in hPa
 343 
c     Check that surface pressure is in hPa
- 
 
 344 
      if ( leveltype.eq.'hybrid_sigma_pressure' ) then
320 
      maxps = -1.e19
 345 
        maxps = -1.e19
321 
      minps =  1.e19
 346 
        minps =  1.e19
322 
      do i=1,vardim(1)
 347 
        do i=1,vardim(1)
323 
        do j=1,vardim(2)
 348 
          do j=1,vardim(2)
324 
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
 349 
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
Line 328... Line 353...
328 
      if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then
 353 
        if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then
329 
         print*,' ERROR: surface pressre PS must be in hPa'
 354 
           print*,' ERROR: surface pressre PS must be in hPa'
330 
         print*,'       ',maxps,minps
 355 
           print*,'       ',maxps,minps
331 
         stop
 356 
           stop
332 
      endif
 357 
        endif
- 
 
 358 
      endif
333 
 
 359 
 
334 
c     ---- Define output of subroutine --------------------------------
 360 
c     ---- Define output of subroutine --------------------------------
335 
 
 361 
 
336 
c     If not full list of vertical levels, reduce AK,BK arrays
 362 
c     If not full list of vertical levels, reduce AK,BK arrays
337 
      if ( (leveltype.eq.'hybrid_sigma_pressure').and.
 363 
      if ( (leveltype.eq.'hybrid_sigma_pressure').and.
Line 545... Line 571...
545 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 571 
        IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
546 
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
 572 
        ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
547 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 573 
        IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
548 
      enddo
 574 
      enddo
549 
 
 575
- 
 
 576 
 
550 
c     Check whether the list of dimensions is OK
 577 
c     Check whether the list of dimensions is OK
- 
 
 578 
      if (ndim .eq. 4) then
551 
      if ( ( dimname(1).ne.'lon'  ).or.
 579 
          if ( ( dimname(1).ne.'lon'  ).or.
552 
     >     ( dimname(2).ne.'lat'  ).or.
 580 
     >       ( dimname(2).ne.'lat'  ).or.
553 
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
 581 
     >       ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
554 
     >     ( dimname(4).ne.'time' ) )
 582 
     >         ( dimname(4).ne.'time' ) )
555 
     >then
 583 
     >    then
556 
        print*,' ERROR: the dimensions of the variable are not correct'
 584 
        print*,' ERROR: the dimensions of the variable are not correct'
557 
        print*,'        expected -> lon / lat / lev / time'
 585 
        print*,'        expected -> lon / lat / lev / time'
558 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
 586 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
559 
        stop
 587 
            stop
560 
      endif
 588 
          endif
- 
 
 589 
      else
- 
 
 590 
          if ( ( dimname(1).ne.'lon'  ).or.
- 
 
 591 
     >         ( dimname(2).ne.'lat'  ).or.
- 
 
 592 
     >         ( dimname(3).ne.'time' ) )
- 
 
 593 
     >    then
- 
 
 594 
        print*,' ERROR: the dimensions of the variable are not correct'
- 
 
 595 
        print*,'        expected -> lon / lat / time'
- 
 
 596 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
- 
 
 597 
            stop
- 
 
 598 
          endif
- 
 
 599 
      end if
- 
 
 600
- 
 
 601 
c     Check about the type of vertical levels
- 
 
 602 
      if (ndim .eq. 4) then
561 
 
 603 
 
562 
c     Get dimension of AK,BK
 604 
c         Get dimension of AK,BK
563 
      varname = 'nhym'
 605 
          varname = 'nhym'
564 
      ierr = NF90_INQ_DIMID(fid,varname,dimid)
 606 
          ierr = NF90_INQ_DIMID(fid,varname,dimid)
565 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 607 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
566 
      ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
 608 
          ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
567 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 609 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
568 
 
 610 
 
569 
c     Check about the type of vertical levels
 - 
 
570 
      varname=dimname(3)
 611 
          varname=dimname(3)
571 
      ierr = NF90_INQ_VARID(fid,varname,varid)
 612 
          ierr = NF90_INQ_VARID(fid,varname,varid)
572 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 613 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
573 
      ierr = nf90_get_att(fid, varid, "standard_name", leveltype)
 614 
          ierr = nf90_get_att(fid, varid, "standard_name", leveltype)
574 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 615 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
575 
      if ( (leveltype.ne.'hybrid_sigma_pressure').and.
 616 
          if ( (leveltype.ne.'hybrid_sigma_pressure').and.
576 
     >     (leveltype.ne.'air_pressure'         ) )
 617 
     >         (leveltype.ne.'air_pressure'         ) ) then
577 
     >then
 - 
 
578 
         print*,' ERROR: input netCDF data must be on hybrid-sigma'
 618 
             print*,' ERROR: input netCDF data must be on hybrid-sigma'
579 
         print*,'        or air pressure levels!',trim(leveltype)
 619 
             print*,'        or air pressure levels!',trim(leveltype)
580 
         stop
 620 
             stop
581 
      endif
 621 
          endif
- 
 
 622 
      else
- 
 
 623 
         leveltype='air_pressure'
- 
 
 624 
      endif
582 
 
 625 
 
583 
c     Allocate memory for reading arrays - depending on <closear>
 626 
c     Allocate memory for reading arrays - depending on <closear>
584 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
 627 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
585 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
 628 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
586 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
 629 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
Line 589... Line 632...
589 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
 632 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
590 
      allocate(lat(vardim(2)),stat=stat)
 633 
      allocate(lat(vardim(2)),stat=stat)
591 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
 634 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
592 
      allocate(lev(vardim(3)),stat=stat)
 635 
      allocate(lev(vardim(3)),stat=stat)
593 
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
 636 
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
- 
 
 637 
      if (ndim .eq. 4) then
594 
      allocate(aktmp(nakbktmp),stat=stat)
 638 
          allocate(aktmp(nakbktmp),stat=stat)
595 
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
 639 
          if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
596 
      allocate(bktmp(nakbktmp),stat=stat)
 640 
          allocate(bktmp(nakbktmp),stat=stat)
597 
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
 641 
          if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
- 
 
 642 
      end if
598 
 
 643 
 
599 
c     Get domain boundaries - longitude, latitude, levels
 644 
c     Get domain boundaries - longitude, latitude, levels
600 
      varname = dimname(1)
 645 
      varname = dimname(1)
601 
      ierr = NF90_INQ_VARID(fid,varname,varid)
 646 
      ierr = NF90_INQ_VARID(fid,varname,varid)
602 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 647 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
Line 616... Line 661...
616 
      varmax(1) = lon( vardim(1) )
 661 
      varmax(1) = lon( vardim(1) )
617 
      varmin(2) = lat(1)
 662 
      varmin(2) = lat(1)
618 
      varmax(2) = lat( vardim(2) )
 663 
      varmax(2) = lat( vardim(2) )
619 
 
 664 
 
620 
c     Get ak and bk
 665 
c     Get ak and bk
- 
 
 666 
      if (ndim .eq. 4) then
621 
      varname='hyam'
 667 
          varname='hyam'
622 
      ierr = NF90_INQ_VARID(fid,varname,varid)
 668 
          ierr = NF90_INQ_VARID(fid,varname,varid)
623 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 669 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
624 
      ierr = nf90_get_var(fid,varid,aktmp)
 670 
          ierr = nf90_get_var(fid,varid,aktmp)
625 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 671 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
Line 638... Line 684...
638 
      if ( units.eq.'Pa' ) then
 684 
          if ( units.eq.'Pa' ) then
639 
         do k=1,nakbktmp
 685 
             do k=1,nakbktmp
640 
            aktmp(k) = 0.01 * aktmp(k)
 686 
                aktmp(k) = 0.01 * aktmp(k)
641 
         enddo
 687 
             enddo
642 
      endif
 688 
          endif
- 
 
 689 
      end if
643 
 
 690 
 
644 
c     Check that unit of lev is in hPa - if necessary correct it
 691 
c     Check that unit of lev is in hPa - if necessary correct it
645 
      if ( leveltype.eq.'air_pressure' ) then
 692 
      if ( leveltype.eq.'air_pressure' .and. vardim(3) .gt. 1) then
646 
         varname='lev'
 693 
         varname='lev'
647 
         ierr = NF90_INQ_VARID(fid,varname,varid)
 694 
         ierr = NF90_INQ_VARID(fid,varname,varid)
648 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 695 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
649 
         ierr = nf90_get_att(fid, varid, "units", units)
 696 
         ierr = nf90_get_att(fid, varid, "units", units)
650 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 697 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
Line 694... Line 741...
694 
           enddo
 741 
           enddo
695 
         endif
 742 
         endif
696 
      endif
 743 
      endif
697
744
698 
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
 745 
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
699 
      if ( vardim(3).eq.1 ) then
 746 
c      if ( vardim(3).eq.1 ) then
700 
         do i=1,vardim(1)
 747 
c         do i=1,vardim(1)
701 
            do j=1,vardim(2)
 748 
c            do j=1,vardim(2)
702 
               do k=1,vardim(3)
 749 
c               do k=1,vardim(3)
703 
                  tmp3(i,j,k) = tmp3(i,j,1)
 750 
c                  tmp3(i,j,k) = tmp3(i,j,1)
704 
               enddo
 751 
c               enddo
705 
            enddo
 752 
c            enddo
706 
         enddo
 753 
c         enddo
707 
      endif
 754 
c      endif
708 
 
 755 
 
709 
c     Decide whether to close arrays
 756 
c     Decide whether to close arrays
710 
      delta = varmax(1)-varmin(1)-360.
 757 
      delta = varmax(1)-varmin(1)-360.
711 
      if (abs(delta+dx).lt.eps) then
 758 
      if (abs(delta+dx).lt.eps) then
712 
          closear = 1
 759 
          closear = 1
713 
      else
 760 
      else
714 
          closear = 0
 761 
          closear = 0
715 
      endif
 762 
      endif
716 
 
 763 
 
717 
c     Save output array - close array and swap on demand
 764 
c     Save output array - close array and swap on demand
- 
 
 765 
 
718 
      do j=1,vardim(2)
 766 
      do j=1,vardim(2)
719 
        do k=1,vardim(3)
 767 
        do k=1,nz
720 
          do i=1,vardim(1)
 768 
          do i=1,vardim(1)
- 
 
 769 
             if (ndim .eq. 4) then
721 
             if ( vertical_swap.eq.1 ) then
 770 
                 if ( vertical_swap.eq.1 ) then
722 
                 field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
 771 
                     field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
723 
             else
 772 
                 else
724 
                 field(i,j,k) = tmp3(i,j,k)
 773 
                     field(i,j,k) = tmp3(i,j,k)
725 
             endif
 774 
                 endif
- 
 
 775 
             else
- 
 
 776 
                field(i,j,k) = tmp3(i,j,1)
- 
 
 777 
             endif
726 
          enddo
 778 
          enddo
727 
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
 779 
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
728 
        enddo
 780 
        enddo
729 
      enddo
 781 
      enddo
730
782
731 
c     Exit point
 783 
c     Exit point
732 
      return
 784 
      return
733
785
734 
      end
 786 
      end
735 
 
 787
- 
 
 788 
 
736 
c     ------------------------------------------------------------
 789 
c     ------------------------------------------------------------
737 
c     Close input file
 790 
c     Close input file
738 
c     ------------------------------------------------------------
 791 
c     ------------------------------------------------------------
739 
 
 792 
 
740 
      subroutine input_close(fid)
 793 
      subroutine input_close(fid)