Subversion Repositories lagranto.ecmwf

Rev

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

Rev 19 Rev 21
Line 29... Line 29...
29 
      subroutine input_open (fid,filename)
 29 
      subroutine input_open (fid,filename)
30 
 
 30 
 
31 
c     Open the input file with filename <filename> and return the
 31 
c     Open the input file with filename <filename> and return the
32 
c     file identifier <fid> for further reference.
32 
c     file identifier <fid> for further reference.
33 
 
 33 
 
- 
 
 34 
      use netcdf
- 
 
 35 
 
34 
      implicit none
 36 
      implicit none
35 
 
 37 
 
36 
c     Declaration of subroutine parameters
 38 
c     Declaration of subroutine parameters
37 
      integer      fid              ! File identifier
 39 
      integer      fid              ! File identifier
38 
      character*80 filename         ! Filename
 40 
      character*80 filename         ! Filename
39 
 
 41 
 
40 
c     Declaration of auxiliary variables
 42 
c     Declaration of auxiliary variables
41 
      integer      ierr
 43 
      integer      ierr
42 
 
 44 
 
43 
c     Open IVE netcdf file
 45 
c     Open netcdf file
44 
      call cdfopn(filename,fid,ierr)
 46 
      ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
45 
      if (ierr.ne.0) goto 900
 47 
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
46
48 
 
47 
c     Exception handling
 49 
c     Exception handling
48 
      return
 50 
      return
49
- 
 
50 
 900  print*,'Cannot open input file ',trim(filename)
 - 
 
51 
      stop
 - 
 
52 
 
 51 
 
53 
      end
 52 
      end
54 
 
 - 
 
55 
 
 53
56 
c     ------------------------------------------------------------
 54 
c     ------------------------------------------------------------
57 
c     Read information about the grid
 55 
c     Read information about the grid
58 
c     ------------------------------------------------------------
 56 
c     ------------------------------------------------------------
59
57
60 
      subroutine input_grid
58 
      subroutine input_grid
Line 70... Line 68...
70 
c     The surface pressure is given in <ps(nx,ny)>. If <fid> is negative,
68 
c     The surface pressure is given in <ps(nx,ny)>. If <fid> is negative,
71 
c     only the grid dimensions and grid parameters (xmin...pollat,nz) are
69 
c     only the grid dimensions and grid parameters (xmin...pollat,nz) are
72 
c     determined and returned (this is needed for dynamical allocation of
70 
c     determined and returned (this is needed for dynamical allocation of
73 
c     memory).
 71 
c     memory).
74 
 
 72 
 
- 
 
 73 
      use netcdf
- 
 
 74 
 
75 
      implicit none
 75 
      implicit none
76 
 
 76 
 
77 
c     Declaration of subroutine parameters
77 
c     Declaration of subroutine parameters
78 
      integer      fid                 ! File identifier
 78 
      integer      fid                 ! File identifier
79 
      real         xmin,xmax,ymin,ymax ! Domain size
 79 
      real         xmin,xmax,ymin,ymax ! Domain size
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         times(10)
 - 
 
104 
      integer      ntimes
 - 
 
105 
      real         aklay(nz),bklay(nz),aklev(nz),bklev(nz)
 103 
      real         aklay(nz),bklay(nz)
106 
      integer      nvars
 104 
      integer      nvars
107 
      character*80 vars(100)
 105 
      character*80 vars(100)
- 
 
 106 
      integer        dimids (nf90_max_var_dims)
- 
 
 107 
      character*80   dimname(nf90_max_var_dims)
- 
 
 108 
      real,allocatable, dimension (:)     :: lon,lat,lev
- 
 
 109 
      real,allocatable, dimension (:)     :: times
- 
 
 110 
      real,allocatable, dimension (:,:)   :: tmp2
- 
 
 111 
      real,allocatable, dimension (:,:,:) :: tmp3
- 
 
 112 
      character*80  units
108 
 
 113 
 
109 
c     Auxiliary varaibles
 114 
c     Auxiliary variables
110 
      integer      ierr
115 
      integer      ierr
111 
      integer      i,j,k
 116 
      integer      i,j,k
112 
      integer      isok
 117 
      integer      isok
113 
      real         tmp(200)
 118 
      real         tmp(200)
114 
      character*80 varname
 119 
      character*80 varname
115 
      real         rtime
 120 
      real         rtime
- 
 
 121 
      integer      varid
116 
      integer      is2d
 122 
      integer      cdfid
117 
      integer      plev
 123 
      integer      stat
118 
 
 - 
 
119 
c     Init the flag for 2D variables - assume a 3D field
 124 
      real         delta
120 
      is2d = 0
 125 
      integer      closear
- 
 
 126 
      real         maxps,minps
121 
 
 127 
 
122 
c     Init the flag for pressure levels (PS is not needed)
 - 
 
123 
      plev = 0
 - 
 
124 
 
 - 
 
125 
c     Inquire dimensions and grid constants if <fid> is negative
 128 
c     ------ Set file identifier --------------------------------------
126 
      if (fid.lt.0) then
 129 
      if (fid.lt.0) then
- 
 
 130 
        cdfid = -fid
- 
 
 131 
      else
- 
 
 132 
        cdfid = fid
- 
 
 133 
      endif
127 
 
 134 
 
128 
c        Get grid info for the selected variable
 135 
c     Special handling if 3D pressure is
129 
         if ( fieldname.eq.'PLEV' ) then
 136 
      if ( fieldname.eq.'P' ) then
130 
            varname = 'PS'
 - 
 
131 
            stagz   = 0.
 - 
 
132 
            call getdef(-fid,varname,ndim,mdv,vardim,
 - 
 
133 
     >                  varmin,varmax,stag,ierr)
 - 
 
134 
            if (ierr.ne.0) goto 900
 - 
 
135 
            call getcfn(-fid,cstfile,ierr)
 - 
 
136 
            if (ierr.ne.0) goto 903
 - 
 
137 
            call cdfopn(cstfile,cstid,ierr)
 - 
 
138 
            if (ierr.ne.0) goto 903
 - 
 
139 
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
 - 
 
140 
            if (ierr.ne.0) goto 903
 - 
 
141 
            call clscdf(cstid,ierr)
 - 
 
142 
            if (ierr.ne.0) goto 903
 - 
 
143
- 
 
144 
         elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then
- 
 
145 
            varname = 'PS'
 137 
         fieldname = 'U'
146 
            stagz   = -0.5
 - 
 
147 
            call getdef(-fid,varname,ndim,mdv,vardim,
 - 
 
148 
     >                  varmin,varmax,stag,ierr)
 - 
 
149 
            if (ierr.ne.0) goto 900
 - 
 
150 
            call getcfn(-fid,cstfile,ierr)
 - 
 
151 
            if (ierr.ne.0) goto 903
 - 
 
152 
            call cdfopn(cstfile,cstid,ierr)
 - 
 
153 
            if (ierr.ne.0) goto 903
 - 
 
154 
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
 - 
 
155 
            if (ierr.ne.0) goto 903
 - 
 
156 
            call clscdf(cstid,ierr)
 - 
 
157 
            if (ierr.ne.0) goto 903
 - 
 
158 
 
 - 
 
159 
         else
 - 
 
160 
            varname = fieldname
 - 
 
161 
            call getdef(-fid,varname,ndim,mdv,vardim,
 - 
 
162 
     >                  varmin,varmax,stag,ierr)
 - 
 
163 
            if (ierr.ne.0) goto 900
 - 
 
164
- 
 
165 
         endif
 - 
 
166 
 
 - 
 
167 
c        Set the grid dimensions and constants - vardim(3) is taken from constants file
 - 
 
168 
         nx   = vardim(1)
 - 
 
169 
         ny   = vardim(2)
 - 
 
170 
         nz   = vardim(3)
 - 
 
171 
         xmin = varmin(1)
 - 
 
172 
         ymin = varmin(2)
 - 
 
173 
         xmax = varmax(1)
 - 
 
174 
         ymax = varmax(2)
 - 
 
175 
         dx   = (xmax-xmin)/real(nx-1)
 - 
 
176 
         dy   = (ymax-ymin)/real(ny-1)
 - 
 
177 
 
 - 
 
178 
c        Get pole position - if no constants file available, set default pole
 - 
 
179 
         call getcfn(-fid,cstfile,ierr)
 - 
 
180 
         if (ierr.eq.0) then
- 
 
181 
            call cdfopn(cstfile,cstid,ierr)
 - 
 
182 
            if (ierr.ne.0) goto 903
 - 
 
183 
            call getpole(cstid,pollon,pollat,ierr)
 - 
 
184 
            if (ierr.ne.0) goto 903
 - 
 
185 
            call clscdf(cstid,ierr)
 - 
 
186 
            if (ierr.ne.0) goto 903
 - 
 
187 
         else
 - 
 
188 
            pollon = 0.
 - 
 
189 
            pollat = 90.
 - 
 
190 
         endif
 138 
      endif
191 
 
 139 
 
192 
c     Get non-constant grid parameters (surface pressure and vertical grid)
 - 
 
193 
      else
 - 
 
194
- 
 
195 
c        Special handling for fieldname 'P.ML' -> in this case the fields
 - 
 
196 
c        P and PS are available on the data file and can be read in. There
 - 
 
197 
c        is no need to reconstruct it from PS,AK and BK. This mode is
 140 
c     Get number of dimensions of variable -> ndim
198 
c        used for model-level (2D) trajectories
 141 
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
199 
         if ( fieldname.eq.'P.ML' ) then
 142 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
200 
 
 - 
 
201 
c           Get the right time to read
 - 
 
202 
            call gettimes(fid,times,ntimes,ierr)
 143 
      ierr = nf90_inquire_variable(cdfid, varid, ndims  = ndim)
203 
            if (ierr.ne.0) goto 901
 - 
 
204 
            isok=0
 - 
 
205 
            do i=1,ntimes
 - 
 
206 
               if (abs(time-times(i)).lt.eps) then
 144 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
207 
                  isok = 1
 145 
      if ( ndim.ne.4 ) then
208 
                  rtime = times(i)
 - 
 
209 
               elseif (timecheck.eq.'no') then
 146 
          print*,' ERROR: netCDF variables need to be 4D'
210 
                  isok = 1
 - 
 
211 
                  rtime = times(1)
 147 
          print*,'      ',trim(fieldname)
212 
               endif
 148 
          stop
213 
            enddo
 149 
      endif
214 
 
 150 
 
215 
c           Read surface pressure and 3D pressure
 151 
c     Get dimensions -> vardim(1:ndim),dimname(1:ndim)
216 
            varname='PS'
 - 
 
217 
            call getdat(fid,varname,rtime,0,ps,ierr)
 152 
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
218 
            if (ierr.ne.0) goto 904
 153 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
219 
            varname='P'
 - 
 
220 
            call getdat(fid,varname,rtime,0,p3,ierr)
 154 
      ierr = nf90_inquire_variable(cdfid, varid,
221 
            if (ierr.ne.0) goto 904
 155 
     >                                   dimids = dimids(1:ndim))
222 
 
 - 
 
223 
c           Set MDV to 1050. - otherwise interpolation routines don't work
 156 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
224 
            do i=1,nx
 157 
      do i=1,ndim
225 
              do j=1,ny
 158 
           ierr = nf90_inquire_dimension(cdfid, dimids(i),
226 
                do k=1,nz
 159 
     >                               name = dimname(i) )
227 
                   if ( p3(i,j,k).lt.0. ) p3(i,j,k) = 1050.
 160 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
228 
                enddo
 161 
           ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))
229 
              enddo
 162 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
230 
            enddo
 163 
      enddo
231 
 
 164 
 
- 
 
 165 
c     Check whether the list of dimensions is OK
- 
 
 166 
      if ( ( dimname(1).ne.'lon'  ).or.
- 
 
 167 
     >     ( dimname(2).ne.'lat'  ).or.
- 
 
 168 
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
- 
 
 169 
     >     ( dimname(4).ne.'time' ) )
- 
 
 170 
     >then
- 
 
 171 
        print*,' ERROR: the dimensions of the variable are not correct'
232 
c           Don't care about other stuff - finish subroutine
 172 
        print*,'        expected -> lon / lat / lev / time'
- 
 
 173 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
233 
            goto 800
 174 
        stop
- 
 
 175 
      endif
234 
 
 176 
 
- 
 
 177 
c     Allocate memory for reading arrays
- 
 
 178 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
- 
 
 179 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
- 
 
 180 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
- 
 
 181 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
- 
 
 182 
      allocate(lon(vardim(1)),stat=stat)
- 
 
 183 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
- 
 
 184 
      allocate(lat(vardim(2)),stat=stat)
- 
 
 185 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
- 
 
 186 
      allocate(times(vardim(4)),stat=stat)
- 
 
 187 
      if (stat.ne.0) print*,'*** error allocating array times   ***'
- 
 
 188 
 
- 
 
 189 
c     Get domain longitudes and latitudes
- 
 
 190 
      varname = dimname(1)
- 
 
 191 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 192 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 193 
      ierr = nf90_get_var(cdfid,varid,lon)
- 
 
 194 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 195 
      varname = dimname(2)
- 
 
 196 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 197 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 198 
      ierr = nf90_get_var(cdfid,varid,lat)
- 
 
 199 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 200
- 
 
 201 
c     Get ak and bk
- 
 
 202 
      varname='hyam'
- 
 
 203 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 204 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 205 
      ierr = nf90_get_var(cdfid,varid,ak)
- 
 
 206 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 207 
      varname='hybm'
- 
 
 208 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 209 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 210 
      ierr = nf90_get_var(cdfid,varid,bk)
- 
 
 211 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 212 
 
- 
 
 213 
c     Check that unit of ak is in hPa - if necessary correct it
- 
 
 214 
      varname='hyam'
- 
 
 215 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 216 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 217 
      ierr = nf90_get_att(cdfid, varid, "units", units)
- 
 
 218 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 219 
      if ( units.eq.'Pa' ) then
- 
 
 220 
      print*,'Hallo'
- 
 
 221 
         do i=1,vardim(3)
- 
 
 222 
            ak(k) = 0.01 * ak(k)
235 
         endif
 223 
         enddo
- 
 
 224 
      endif
236 
 
 225 
 
237 
c        Get full grid info - in particular staggering flag; set flag for 2D tracing
 - 
 
238 
         if ( fieldname.eq.'PLEV' ) then
 - 
 
239 
            varname = 'PS'
 - 
 
240 
            stagz   = 0.
 - 
 
241 
            call getdef(fid,varname,ndim,mdv,vardim,
 - 
 
242 
     >                  varmin,varmax,stag,ierr)
 - 
 
243 
            if (ierr.ne.0) goto 900
 - 
 
244 
            call getcfn(fid,cstfile,ierr)
 226 
c     Get time information (check if time is correct)
245 
            if (ierr.ne.0) goto 903
 - 
 
246 
            call cdfopn(cstfile,cstid,ierr)
 - 
 
247 
            if (ierr.ne.0) goto 903
 - 
 
248 
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
 - 
 
249 
            if (ierr.ne.0) goto 903
 - 
 
250 
            call clscdf(cstid,ierr)
 - 
 
251 
            if (ierr.ne.0) goto 903
 - 
 
252
- 
 
253 
         elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then
- 
 
254 
            varname = 'PS'
 227 
      varname = 'time'
255 
            stagz   = -0.5
 - 
 
256 
            call getdef(fid,varname,ndim,mdv,vardim,
 228 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
257 
     >                  varmin,varmax,stag,ierr)
 229 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
258 
            if (ierr.ne.0) goto 900
 - 
 
259 
            call getcfn(fid,cstfile,ierr)
 230 
      ierr = nf90_get_var(cdfid,varid,times)
260 
            if (ierr.ne.0) goto 903
 - 
 
261 
            call cdfopn(cstfile,cstid,ierr)
 - 
 
262 
            if (ierr.ne.0) goto 903
 - 
 
263 
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
 - 
 
264 
            if (ierr.ne.0) goto 903
 - 
 
265 
            call clscdf(cstid,ierr)
 - 
 
266 
            if (ierr.ne.0) goto 903
 - 
 
267 
 
 - 
 
268 
         else
 - 
 
269 
            varname=fieldname
 - 
 
270 
            call getdef(fid,varname,ndim,mdv,vardim,
 - 
 
271 
     >                  varmin,varmax,stag,ierr)
 231 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
272 
            if (ierr.ne.0) goto 900
 - 
 
273 
            if (vardim(3).eq.1) is2d = 1
 - 
 
274 
         endif
 - 
 
275 
 
 - 
 
276 
c        Get time information (check if time is correct)
 - 
 
277 
         call gettimes(fid,times,ntimes,ierr)
 - 
 
278 
         if (ierr.ne.0) goto 901
 - 
 
279 
         isok=0
 232 
      isok=0
280 
         do i=1,ntimes
 233 
      do i=1,vardim(4)
281 
            if (abs(time-times(i)).lt.eps) then
 234 
        if (abs(time-times(i)).lt.eps) then
282 
               isok = 1
 235 
               isok = 1
283 
               rtime = times(i)
 236 
               rtime = times(i)
284 
            elseif (timecheck.eq.'no') then
 237 
        elseif (timecheck.eq.'no') then
285 
               isok = 1
 238 
               isok = 1
286 
               rtime = times(1)
 239 
               rtime = times(1)
287 
            endif
 240 
        endif
288 
         enddo
 241 
      enddo
289 
         if ( isok.eq.0) goto 905
 242 
      if ( isok.eq.0 ) then
- 
 
 243 
         print*,' ERROR: time ',rtime,' not found on netCDF file'
- 
 
 244 
         stop
- 
 
 245 
      endif
290 
 
 246 
 
- 
 
 247 
c     Read surface pressure
- 
 
 248 
      varname='PS'
- 
 
 249 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 250 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 251 
      ierr = nf90_get_var(cdfid,varid,tmp2)
- 
 
 252 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 253
- 
 
 254 
c     Check that surface pressure is in hPa
- 
 
 255 
      maxps = -1.e39
- 
 
 256 
      minps =  1.e39
- 
 
 257 
      do i=1,vardim(1)
- 
 
 258 
        do j=1,vardim(2)
- 
 
 259 
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
- 
 
 260 
             if (tmp2(i,j).lt.minps) minps = tmp2(i,j)
- 
 
 261 
        enddo
- 
 
 262 
      enddo
- 
 
 263 
      if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then
291 
c        If 2D tracing requested: take dummay values for PS, AKLEV,BKLEV,AKLAY,BKLAY
 264 
         print*,' ERROR: surface pressre PS must be in hPa'
292 
         if ( is2d.eq.1 ) then
 265 
         print*,'       ',maxps,minps
293
266 
         stop
- 
 
 267 
      endif
- 
 
 268 
 
- 
 
 269 
c     Calculate layer and level pressures
- 
 
 270 
      do i=1,vardim(1)
294 
            do i=1,nx
 271 
         do j=1,vardim(2)
295 
               do j=1,ny
 272 
               do k=1,vardim(3)
296 
                  ps(i,j) = 1050.
 273 
                  tmp3(i,j,k)=ak(k)+bk(k)*tmp2(i,j)
297 
               enddo
 274 
               enddo
298 
            enddo
 275 
         enddo
299
- 
 
300 
            do k=1,nz
 - 
 
301 
               aklev(k) = 0.
 - 
 
302 
               bklev(k) = real(nz-k)/real(nz-1)
- 
 
303 
               aklay(k) = 0.
 - 
 
304 
               bklay(k) = real(nz-k)/real(nz-1)
- 
 
305 
            enddo
 276 
      enddo
306 
 
 277 
 
307 
c        3D tracing - read PS, AKLEV,BKLEV,AKLAY;BKLAY
 278 
c     Set the grid dimensions and constants
- 
 
 279 
      nx      = vardim(1)
- 
 
 280 
      ny      = vardim(2)
- 
 
 281 
      nz      = vardim(3)
- 
 
 282 
      xmin    = lon(1)
- 
 
 283 
      ymin    = lat(1)
- 
 
 284 
      xmax    = lon(nx)
- 
 
 285 
      ymax    = lat(ny)
- 
 
 286 
      dx      = (xmax-xmin)/real(nx-1)
- 
 
 287 
      dy      = (ymax-ymin)/real(ny-1)
- 
 
 288 
      pollon  = 0.
- 
 
 289 
      pollat  = 90.
- 
 
 290 
      stagz   = 0.
- 
 
 291 
      delta   = xmax-xmin-360.
- 
 
 292 
      if (abs(delta+dx).lt.eps) then
- 
 
 293 
          xmax    = xmax + dx
- 
 
 294 
          nx      = nx + 1
- 
 
 295 
          closear = 1
308 
         else
 296 
      else
- 
 
 297 
          closear = 0
- 
 
 298 
      endif
309 
 
 299 
 
310 
c           Read the level coefficients from the constants file
 - 
 
311 
            call getcfn(fid,cstfile,ierr)
 - 
 
312 
            if (ierr.ne.0) goto 903
 - 
 
313 
            call cdfopn(cstfile,cstid,ierr)
 - 
 
314 
            if (ierr.ne.0) goto 903
 - 
 
315 
            call getlevs(cstid,vardim(3),aklev,bklev,aklay,bklay,ierr)
 - 
 
316 
            if (ierr.ne.0) goto 903
 - 
 
317 
            call clscdf(cstid,ierr)
 - 
 
318 
            if (ierr.ne.0) goto 903
 - 
 
319 
 
 - 
 
320 
c           Check whether PS is needed to get the 3d pressure field
 300 
c     Save the output arrays (if fid>0) - close arrays on request
321 
            plev = 1
 - 
 
322 
            do i=1,nz
 301 
      if ( fid.gt.0 ) then
323 
              if ( (abs(stagz).lt.eps).and.(abs(bklev(i)).gt.eps) ) then
 - 
 
324 
                plev = 0
 - 
 
325 
              endif
 - 
 
326 
              if ( (abs(stagz).gt.eps).and.(abs(bklay(i)).gt.eps) ) then
 - 
 
327 
                plev = 0
 - 
 
328 
              endif
 - 
 
329 
            enddo
 - 
 
330 
 
 302 
 
331 
c           Read surface pressure if needed
 - 
 
332 
            if ( plev.ne.1 ) then
 - 
 
333 
              varname='PS'
 303 
         do j=1,vardim(2)
334 
              call getdat(fid,varname,rtime,0,ps,ierr)
 - 
 
335 
              if (ierr.ne.0) goto 904
 - 
 
336 
            else
 - 
 
337 
              do i=1,nx
 304 
           do i=1,vardim(1)
338 
                do j=1,ny
 - 
 
339 
                    ps(i,j) = 1000.
 305 
             ps(i,j) = tmp2(i,j)
340 
                enddo
 - 
 
341 
              enddo
 306 
           enddo
342 
            endif
 - 
 
343 
 
 - 
 
344 
         endif
 - 
 
345 
 
 - 
 
346 
c        Calculate layer and level pressures
 - 
 
347 
         do i=1,nx
 - 
 
348 
            do j=1,ny
 - 
 
349 
               do k=1,nz
 - 
 
350 
                  if ( abs(stagz).lt.eps ) then
 - 
 
351 
                     p3(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
 - 
 
352 
                  else
 - 
 
353 
                     p3(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
 307 
           if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)
354 
                  endif
 - 
 
355 
               enddo
 - 
 
356 
            enddo
 - 
 
357 
         enddo
 308 
         enddo
358 
 
 309 
 
359 
c        Set the ak and bk for the vertical grid
 - 
 
360 
         do k=1,nz
 310 
         do j=1,vardim(2)
361 
            if ( abs(stagz).lt.eps ) then
 311 
           do k=1,vardim(3)
362 
               ak(k)=aklev(k)
 312 
             do i=1,vardim(1)
363 
               bk(k)=bklev(k)
 313 
               p3(i,j,k) = tmp3(i,j,k)
364 
            else
 314 
             enddo
365 
               ak(k)=aklay(k)
 315 
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
366 
               bk(k)=bklay(k)
 - 
 
367 
            endif
 316 
           enddo
368 
         enddo
 317 
         enddo
369 
 
 - 
 
370 
      endif
 318 
      endif
371
319 
 
372 
c     Exit point for subroutine
 - 
 
373 
 800  continue
 - 
 
374 
      return
 320 
      return
375
321
376 
c     Exception handling
 - 
 
377 
 900  print*,'Cannot retrieve grid dimension from ',fid
 - 
 
378 
      stop
 - 
 
379 
 901  print*,'Cannot retrieve grid parameters from ',fid
 - 
 
380 
      stop
 - 
 
381 
 902  print*,'Grid inconsistency detected for ',fid
 - 
 
382 
      stop
 - 
 
383 
 903  print*,'Problem with level coefficients from ',trim(cstfile)
 - 
 
384 
      stop
 - 
 
385 
 904  print*,'Cannot read surface pressure from ',trim(cstfile)
 - 
 
386 
      stop
 - 
 
387 
 905  print*,'Cannot find time ',time,' on ',fid
 - 
 
388 
      stop
 - 
 
389 
 906  print*,'Unable to get grid info ',fid
 - 
 
390 
      stop
 - 
 
391
- 
 
392 
      end
 322 
      end
393 
 
 323 
 
394 
c     ------------------------------------------------------------
 324 
c     ------------------------------------------------------------
395 
c     Read wind information
 325 
c     Read wind information
396 
c     ------------------------------------------------------------
 326 
c     ------------------------------------------------------------
Line 404... Line 334...
404 
c     information is provided in <stagz> and gives the reference to either
 334 
c     information is provided in <stagz> and gives the reference to either
405 
c     the layer or level field from <input_grid>. A consistency check is
 335 
c     the layer or level field from <input_grid>. A consistency check is
406 
c     performed to have an agreement with the grid specified by <xmin,xmax,
 336 
c     performed to have an agreement with the grid specified by <xmin,xmax,
407 
c     ymin,ymax,dx,dy,nx,ny,nz>.
 337 
c     ymin,ymax,dx,dy,nx,ny,nz>.
408 
 
 338 
 
- 
 
 339 
      use netcdf
- 
 
 340 
 
409 
      implicit none
 341 
      implicit none
410 
 
 342 
 
411 
c     Declaration of variables and parameters
 343 
c     Declaration of variables and parameters
412 
      integer      fid                 ! File identifier
 344 
      integer      fid                 ! File identifier
413 
      character*80 fieldname           ! Name of the wind field
 345 
      character*80 fieldname           ! Name of the wind field
Line 437... Line 369...
437 
      integer      ntimes
 369 
      integer      ntimes
438 
      character*80 cstfile
 370 
      character*80 cstfile
439 
      integer      cstid
 371 
      integer      cstid
440 
      real         aklay(200),bklay(200),aklev(200),bklev(200)
 372 
      real         aklay(200),bklay(200),aklev(200),bklev(200)
441 
      real         ps(nx,ny)
 373 
      real         ps(nx,ny)
- 
 
 374 
      integer      dimids (nf90_max_var_dims)
- 
 
 375 
      character*80 dimname(nf90_max_var_dims)
- 
 
 376 
      integer      varid
- 
 
 377 
      integer      cdfid
- 
 
 378 
      real,allocatable, dimension (:)     :: lon,lat
- 
 
 379 
      real,allocatable, dimension (:,:)   :: tmp2
- 
 
 380 
      real,allocatable, dimension (:,:,:) :: tmp3
442 
 
 381 
 
443 
c     Auxiliary variables
 382 
c     Auxiliary variables
444 
      integer      isok
 383 
      integer      isok
445 
      integer      i,j,k
 384 
      integer      i,j,k
446 
      integer      nz1
 385 
      integer      nz1
447 
      real         rtime
 386 
      real         rtime
448 
 
 - 
 
449 
c     Read variable definition - for P, PLEV and PLAY: load also ak,bk
 - 
 
450 
      if ( ( fieldname.eq.'PLEV' ).or.
 387 
      integer      closear
451 
     >     ( fieldname.eq.'PLAY' ).or.
 388 
      integer      stat
452 
     >     ( fieldname.eq.'P'    ) )
 389 
      real         delta
453 
     >then
 390 
 
454 
         call getcfn(fid,cstfile,ierr)
 391 
c     Get the number of dimensions -> ndim
455 
         if (ierr.ne.0) goto 905
 - 
 
456 
         call cdfopn(cstfile,cstid,ierr)
 392 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
457 
         if (ierr.ne.0) goto 905
 - 
 
458 
         call getlevs(cstid,nz1,aklev,bklev,aklay,bklay,ierr)
 393 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
459 
         if (ierr.ne.0) goto 905
 - 
 
460 
         call clscdf(cstid,ierr)
 - 
 
461 
         if (ierr.ne.0) goto 905
 - 
 
462 
         varname = 'PS'
 - 
 
463 
         call getdef(fid,varname,ndim,mdv,vardim,
 394 
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
464 
     >               varmin,varmax,stag,ierr)
 395 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
465 
         vardim(3) = nz1
 - 
 
466 
         if (ierr.ne.0) goto 900
 - 
 
467 
 
 396 
 
468 
      else
 - 
 
469 
         varname = fieldname
 397 
c     Get the dimensions of the arrays -> varid(1:ndim)
470 
         call getdef(fid,varname,ndim,mdv,vardim,
 398 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
471 
     >               varmin,varmax,stag,ierr)
 399 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
472 
         if (ierr.ne.0) goto 900
 400 
      ierr = nf90_inquire_variable(fid, varid,
473 
         stagz=stag(3)
 - 
 
474 
      endif
 - 
 
475 
 
 - 
 
476 
c     Get time information (set time to first one in the file)
 401 
     >                                   dimids = dimids(1:ndim))
477 
      call gettimes(fid,times,ntimes,ierr)
 402 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
478 
      if (ierr.ne.0) goto 902
 - 
 
479 
      isok=0
 - 
 
480 
      do i=1,ntimes
 403 
      do i=1,ndim
481 
         if (abs(time-times(i)).lt.eps) then
 404 
           ierr = nf90_inquire_dimension(fid, dimids(i),
482 
            isok = 1
 - 
 
483 
            rtime = times(i)
 405 
     >                               name = dimname(i) )
484 
         elseif (timecheck.eq.'no') then
 406 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
485 
            isok = 1
 407 
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
486 
            rtime = times(1)
 408 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
487 
         endif
 - 
 
488 
      enddo
 409 
      enddo
489 
      if ( isok.eq.0 ) goto 904
 - 
 
490 
 
 - 
 
491 
c     Read  field
 - 
 
492 
      if ( ( fieldname.eq.'P' ).or.(fieldname.eq.'PLAY') )  then       ! P, PLAY
 - 
 
493 
         stagz   = -0.5
 - 
 
494 
         varname = 'PS'
 - 
 
495 
         call getdat(fid,varname,rtime,0,ps,ierr)
 - 
 
496 
         if (ierr.ne.0) goto 903
 - 
 
497 
         do i=1,nx
 - 
 
498 
            do j=1,ny
 - 
 
499 
               do k=1,nz
 - 
 
500 
                  field(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
 - 
 
501 
               enddo
 - 
 
502 
            enddo
 - 
 
503 
         enddo
 - 
 
504
- 
 
505 
      elseif ( fieldname.eq.'PLEV' )  then                             ! PLEV
 - 
 
506 
         stagz   = 0.
 - 
 
507 
         varname = 'PS'
 - 
 
508 
         call getdat(fid,varname,rtime,0,ps,ierr)
 - 
 
509 
         if (ierr.ne.0) goto 903
 - 
 
510 
         do i=1,nx
 - 
 
511 
            do j=1,ny
 - 
 
512 
               do k=1,nz
 - 
 
513 
                  field(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
 - 
 
514 
               enddo
 - 
 
515 
            enddo
 - 
 
516 
         enddo
 - 
 
517 
 
 - 
 
518 
      else                                                             ! Other fields
 - 
 
519 
         varname=fieldname
 - 
 
520 
         call getdat(fid,varname,rtime,0,field,ierr)
 - 
 
521 
         if (ierr.ne.0) goto 903
 - 
 
522 
 
 410 
 
- 
 
 411 
c     Check whether the list of dimensions is OK
- 
 
 412 
      if ( ( dimname(1).ne.'lon'  ).or.
- 
 
 413 
     >     ( dimname(2).ne.'lat'  ).or.
- 
 
 414 
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
- 
 
 415 
     >     ( dimname(4).ne.'time' ) )
- 
 
 416 
     >then
- 
 
 417 
        print*,' ERROR: the dimensions of the variable are not correct'
- 
 
 418 
        print*,'        expected -> lon / lat / lev / time'
- 
 
 419 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
- 
 
 420 
        stop
523 
      endif
 421 
      endif
524 
 
 422 
 
- 
 
 423 
c     Allocate memory for reading arrays - depending on <closear>
- 
 
 424 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
- 
 
 425 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
- 
 
 426 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
- 
 
 427 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
- 
 
 428 
      allocate(lon(vardim(1)),stat=stat)
- 
 
 429 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
- 
 
 430 
      allocate(lat(vardim(2)),stat=stat)
- 
 
 431 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
- 
 
 432 
 
- 
 
 433 
c     Get domain boundaries
- 
 
 434 
      varname = dimname(1)
- 
 
 435 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 436 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 437 
      ierr = nf90_get_var(fid,varid,lon)
- 
 
 438 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 439 
      varname = dimname(2)
- 
 
 440 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 441 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 442 
      ierr = nf90_get_var(fid,varid,lat)
- 
 
 443 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 444 
 
- 
 
 445 
c     Read data
- 
 
 446 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
- 
 
 447 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 448 
      ierr = nf90_get_var(fid,varid,tmp3)
- 
 
 449 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 450
525 
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
 451 
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
526 
      if ( vardim(3).eq.1 ) then
 452 
      if ( vardim(3).eq.1 ) then
527 
         do i=1,nx
 453 
         do i=1,vardim(1)
528 
            do j=1,ny
 454 
            do j=1,vardim(2)
529 
               do k=1,nz
 455 
               do k=1,vardim(3)
530 
                  field(i,j,k) = field(i,j,1)
 456 
                  tmp3(i,j,k) = tmp3(i,j,1)
531 
               enddo
 457 
               enddo
532 
            enddo
 458 
            enddo
533 
         enddo
 459 
         enddo
534 
      endif
 460 
      endif
535
- 
 
536 
 
 461 
 
- 
 
 462 
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.
- 
 
 465 
      if (abs(delta+dx).lt.eps) then
- 
 
 466 
          closear = 1
- 
 
 467 
      else
- 
 
 468 
          closear = 0
- 
 
 469 
      endif
537 
 
 470 
 
- 
 
 471 
      do j=1,vardim(2)
- 
 
 472 
        do k=1,vardim(3)
- 
 
 473 
          do i=1,vardim(1)
- 
 
 474 
             field(i,j,k) = tmp3(i,j,k)
- 
 
 475 
          enddo
- 
 
 476 
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
- 
 
 477 
        enddo
- 
 
 478 
      enddo
- 
 
 479
538 
c     Exception handling
 480 
c     Exit point
539 
      return
 481 
      return
540
482
541 
 900  print*,'Cannot retrieve definition for ',trim(varname),'  ',fid
 - 
 
542 
      stop
 - 
 
543 
 901  print*,'Grid inconsistency detected for ',trim(varname),'  ',fid
 - 
 
544 
      stop
 - 
 
545 
 902  print*,'Cannot retrieve time for ',trim(varname),'  ',fid
 - 
 
546 
      stop
 - 
 
547 
 903  print*,'Cannot load wind component ',trim(varname),'  ',fid
 - 
 
548 
      stop
 - 
 
549 
 904  print*,'Cannot load time ',time,' for ',trim(varname),'  ',fid
 - 
 
550 
      stop
 - 
 
551 
 905  print*,'Cannot load time vertical grid AK, BK from file  ',fid
 - 
 
552 
      stop
 - 
 
553
- 
 
554 
      end
 483 
      end
555 
 
 484 
 
556 
c     ------------------------------------------------------------
 485 
c     ------------------------------------------------------------
557 
c     Close input file
 486 
c     Close input file
558 
c     ------------------------------------------------------------
 487 
c     ------------------------------------------------------------
559 
 
 488 
 
560 
      subroutine input_close(fid)
 489 
      subroutine input_close(fid)
561 
 
 490 
 
562 
c     Close the input file with file identifier <fid>.
 491 
c     Close the input file with file identifier <fid>.
563 
 
 492 
 
- 
 
 493 
      use netcdf
- 
 
 494 
 
564 
      implicit none
 495 
      implicit none
565 
 
 496 
 
566 
c     Declaration of subroutine parameters
 497 
c     Declaration of subroutine parameters
567 
      integer fid
 498 
      integer fid
568 
 
 499 
 
569 
c     Auxiliary variables
 500 
c     Auxiliary variables
570 
      integer ierr
 501 
      integer ierr
571 
 
 502 
 
572 
c     Close file
 503 
c     Close file
- 
 
 504 
      ierr = NF90_CLOSE(fid)
- 
 
 505 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 506
- 
 
 507 
      end
- 
 
 508
- 
 
 509 
c     ------------------------------------------------------------
- 
 
 510 
c     Get a list of variables on netCDF file
- 
 
 511 
c     ------------------------------------------------------------
- 
 
 512 
 
- 
 
 513 
      subroutine input_getvars(fid,vnam,nvars)
- 
 
 514 
 
- 
 
 515 
c     List of variables on netCDF file
- 
 
 516 
 
- 
 
 517 
      use netcdf
- 
 
 518 
 
- 
 
 519 
      implicit none
- 
 
 520 
 
- 
 
 521 
c     Declaration of subroutine parameters
- 
 
 522 
      integer      fid
- 
 
 523 
      integer      nvars
573 
      call clscdf(fid,ierr)
 524 
      character*80 vnam(200)
- 
 
 525 
 
- 
 
 526 
c     Auxiliary variables
- 
 
 527 
      integer ierr
- 
 
 528 
      integer i
- 
 
 529 
      integer nDims,nGlobalAtts,unlimdimid
- 
 
 530 
 
- 
 
 531 
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
- 
 
 532 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 533
- 
 
 534 
      do i=1,nVars
- 
 
 535 
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
- 
 
 536 
      enddo
574
537
575 
      end
 538 
      end