WebSVN – lagranto.ecmwf – Diff – /trunk/lib/ioinp_cdo.f



      real         mdv
      real         stag(4)
      integer      ndim
      character*80 cstfile
      integer      cstid
 
      integer      nvars
      character*80 vars(100)
      integer        dimids (nf90_max_var_dims)
      character*80   dimname(nf90_max_var_dims)
      real,allocatable, dimension (:)     :: lon,lat,lev
      real,allocatable, dimension (:)     :: times
      real,allocatable, dimension (:,:)   :: tmp2
      real,allocatable, dimension (:,:,:) :: tmp3
      real,allocatable, dimension (:)     :: aktmp,bktmp
      character*80  units
 
c     Auxiliary variables
      integer      ierr       
      integer      i,j,k

      if (stat.ne.0) print*,'*** error allocating array lon     ***' 
      allocate(lat(vardim(2)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lat     ***' 
      allocate(times(vardim(4)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array times   ***'
      allocate(aktmp(vardim(3)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
      allocate(bktmp(vardim(3)),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
 
c     Get domain longitudes and latitudes
      varname = dimname(1)
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)

      
c     Get ak and bk
      varname='hyam'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,aktmp)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      varname='hybm'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_var(cdfid,varid,bktmp)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
 
c     Check that unit of ak is in hPa - if necessary correct it
      varname='hyam'
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      ierr = nf90_get_att(cdfid, varid, "units", units)
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      if ( units.eq.'Pa' ) then
 
         do k=1,vardim(3)
            aktmp(k) = 0.01 * aktmp(k)
         enddo
      endif
 
c     Get time information (check if time is correct)
      varname = 'time'

 
c     Calculate layer and level pressures
      do i=1,vardim(1)
         do j=1,vardim(2)
               do k=1,vardim(3)
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
               enddo
         enddo
      enddo
 
c     Set the grid dimensions and constants

         enddo
 
         do j=1,vardim(2)
           do k=1,vardim(3)
             do i=1,vardim(1)
               p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)
             enddo
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
           enddo
         enddo
 
         do k=1,vardim(3)
            ak(k) = aktmp(vardim(3)-k+1)
            bk(k) = bktmp(vardim(3)-k+1)
         enddo
 
      endif
 
      return
      
      end

            enddo
         enddo
      endif
 
c     Save the ouput array - close on request
 
      delta = varmax(1)-varmin(1)-360.
      if (abs(delta+dx).lt.eps) then
          closear = 1
      else
          closear = 0
      endif
 
      do j=1,vardim(2)
        do k=1,vardim(3)
          do i=1,vardim(1)
             field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
          enddo
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
        enddo
      enddo
         

Rev 21	Rev 23
Line 98...	Line 98...
98	`real mdv`	98	`real mdv`
99	`real stag(4)`	99	`real stag(4)`
100	`integer ndim`	100	`integer ndim`
101	`character*80 cstfile`	101	`character*80 cstfile`
102	`integer cstid`	102	`integer cstid`
103	`real aklay(nz),bklay(nz)`	-
104	`integer nvars`	103	`integer nvars`
105	`character*80 vars(100)`	104	`character*80 vars(100)`
106	`integer dimids (nf90_max_var_dims)`	105	`integer dimids (nf90_max_var_dims)`
107	`character*80 dimname(nf90_max_var_dims)`	106	`character*80 dimname(nf90_max_var_dims)`
108	`real,allocatable, dimension (:) :: lon,lat,lev`	107	`real,allocatable, dimension (:) :: lon,lat,lev`
109	`real,allocatable, dimension (:) :: times`	108	`real,allocatable, dimension (:) :: times`
110	`real,allocatable, dimension (:,:) :: tmp2`	109	`real,allocatable, dimension (:,:) :: tmp2`
111	`real,allocatable, dimension (:,:,:) :: tmp3`	110	`real,allocatable, dimension (:,:,:) :: tmp3`
-		111	`real,allocatable, dimension (:) :: aktmp,bktmp`
112	`character*80 units`	112	`character*80 units`
113		113
114	`c Auxiliary variables`	114	`c Auxiliary variables`
115	`integer ierr`	115	`integer ierr`
116	`integer i,j,k`	116	`integer i,j,k`
Line 183...	Line 183...
183	`if (stat.ne.0) print,' error allocating array lon *'`	183	`if (stat.ne.0) print,' error allocating array lon *'`
184	`allocate(lat(vardim(2)),stat=stat)`	184	`allocate(lat(vardim(2)),stat=stat)`
185	`if (stat.ne.0) print,' error allocating array lat *'`	185	`if (stat.ne.0) print,' error allocating array lat *'`
186	`allocate(times(vardim(4)),stat=stat)`	186	`allocate(times(vardim(4)),stat=stat)`
187	`if (stat.ne.0) print,' error allocating array times *'`	187	`if (stat.ne.0) print,' error allocating array times *'`
-		188	`allocate(aktmp(vardim(3)),stat=stat)`
-		189	`if (stat.ne.0) print,' error allocating array aktmp *'`
-		190	`allocate(bktmp(vardim(3)),stat=stat)`
-		191	`if (stat.ne.0) print,' error allocating array bktmp *'`
188		192
189	`c Get domain longitudes and latitudes`	193	`c Get domain longitudes and latitudes`
190	`varname = dimname(1)`	194	`varname = dimname(1)`
191	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	195	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
192	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	196	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
Line 200...	Line 204...
200		204
201	`c Get ak and bk`	205	`c Get ak and bk`
202	`varname='hyam'`	206	`varname='hyam'`
203	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	207	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
204	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	208	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
205	`ierr = nf90_get_var(cdfid,varid,ak)`	209	`ierr = nf90_get_var(cdfid,varid,aktmp)`
206	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	210	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
207	`varname='hybm'`	211	`varname='hybm'`
208	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	212	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
209	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	213	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
210	`ierr = nf90_get_var(cdfid,varid,bk)`	214	`ierr = nf90_get_var(cdfid,varid,bktmp)`
211	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	215	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
212		216
213	`c Check that unit of ak is in hPa - if necessary correct it`	217	`c Check that unit of ak is in hPa - if necessary correct it`
214	`varname='hyam'`	218	`varname='hyam'`
215	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`	219	`ierr = NF90_INQ_VARID(cdfid,varname,varid)`
216	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	220	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
217	`ierr = nf90_get_att(cdfid, varid, "units", units)`	221	`ierr = nf90_get_att(cdfid, varid, "units", units)`
218	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`	222	`IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)`
219	`if ( units.eq.'Pa' ) then`	223	`if ( units.eq.'Pa' ) then`
220	`print*,'Hallo'`	-
221	`do i=1,vardim(3)`	224	`do k=1,vardim(3)`
222	`ak(k) = 0.01 * ak(k)`	225	`aktmp(k) = 0.01 * aktmp(k)`
223	`enddo`	226	`enddo`
224	`endif`	227	`endif`
225		228
226	`c Get time information (check if time is correct)`	229	`c Get time information (check if time is correct)`
227	`varname = 'time'`	230	`varname = 'time'`
Line 268...	Line 271...
268		271
269	`c Calculate layer and level pressures`	272	`c Calculate layer and level pressures`
270	`do i=1,vardim(1)`	273	`do i=1,vardim(1)`
271	`do j=1,vardim(2)`	274	`do j=1,vardim(2)`
272	`do k=1,vardim(3)`	275	`do k=1,vardim(3)`
273	`tmp3(i,j,k)=ak(k)+bk(k)*tmp2(i,j)`	276	`tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)`
274	`enddo`	277	`enddo`
275	`enddo`	278	`enddo`
276	`enddo`	279	`enddo`
277		280
278	`c Set the grid dimensions and constants`	281	`c Set the grid dimensions and constants`
Line 308...	Line 311...
308	`enddo`	311	`enddo`
309		312
310	`do j=1,vardim(2)`	313	`do j=1,vardim(2)`
311	`do k=1,vardim(3)`	314	`do k=1,vardim(3)`
312	`do i=1,vardim(1)`	315	`do i=1,vardim(1)`
313	`p3(i,j,k) = tmp3(i,j,k)`	316	`p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)`
314	`enddo`	317	`enddo`
315	`if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)`	318	`if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)`
316	`enddo`	319	`enddo`
317	`enddo`	320	`enddo`
-		321
-		322	`do k=1,vardim(3)`
-		323	`ak(k) = aktmp(vardim(3)-k+1)`
-		324	`bk(k) = bktmp(vardim(3)-k+1)`
-		325	`enddo`
-		326
318	`endif`	327	`endif`
319		328
320	`return`	329	`return`
321		330
322	`end`	331	`end`
Line 458...	Line 467...
458	`enddo`	467	`enddo`
459	`enddo`	468	`enddo`
460	`endif`	469	`endif`
461		470
462	`c Save the ouput array - close on request`	471	`c Save the ouput array - close on request`
463	`dx = (varmax(1)-varmin(1))/real(vardim(1)-1)`	-
464	`delta = varmax(1)-varmin(1)-360.`	472	`delta = varmax(1)-varmin(1)-360.`
465	`if (abs(delta+dx).lt.eps) then`	473	`if (abs(delta+dx).lt.eps) then`
466	`closear = 1`	474	`closear = 1`
467	`else`	475	`else`
468	`closear = 0`	476	`closear = 0`
469	`endif`	477	`endif`
470		478
471	`do j=1,vardim(2)`	479	`do j=1,vardim(2)`
472	`do k=1,vardim(3)`	480	`do k=1,vardim(3)`
473	`do i=1,vardim(1)`	481	`do i=1,vardim(1)`
474	`field(i,j,k) = tmp3(i,j,k)`	482	`field(i,j,k) = tmp3(i,j,vardim(3)-k+1)`
475	`enddo`	483	`enddo`
476	`if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)`	484	`if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)`
477	`enddo`	485	`enddo`
478	`enddo`	486	`enddo`
479		487

Subversion Repositories lagranto.ecmwf

(root)/trunk/lib/ioinp_cdo.f – Rev 21 → 23