Subversion Repositories lagranto.ecmwf

Rev

Rev 46 | Show entire file | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 46 Rev 47
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
51 
 
50 
 900  print*,'Cannot open input file ',trim(filename)
 - 
 
51 
      stop
 - 
 
52 
 
 - 
 
53 
      end
 52 
      end
54 
 
 53
55 
 
 - 
 
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)
 - 
 
106 
      integer      nvars
 103 
      integer      nvars
107 
      character*80 vars(100)
 104 
      character*80 vars(100)
- 
 
 105 
      integer        dimids (nf90_max_var_dims),dimid
- 
 
 106 
      character*80   dimname(nf90_max_var_dims)
- 
 
 107 
      character*80   stdname
- 
 
 108 
      real,allocatable, dimension (:)     :: lon,lat,lev
- 
 
 109 
      real,allocatable, dimension (:)     :: times
- 
 
 110 
      real,allocatable, dimension (:,:)   :: tmp2
- 
 
 111 
      real,allocatable, dimension (:,:,:) :: tmp3
- 
 
 112 
      real,allocatable, dimension (:)     :: aktmp,bktmp
- 
 
 113 
      character*80  units
- 
 
 114 
      character*80  leveltype
- 
 
 115 
      integer       nakbktmp
- 
 
 116 
      integer       vertical_swap
108 
 
 117 
 
109 
c     Auxiliary varaibles
 118 
c     Auxiliary variables
110 
      integer      ierr
119 
      integer      ierr
111 
      integer      i,j,k
 120 
      integer      i,j,k
112 
      integer      isok
 121 
      integer      isok
113 
      real         tmp(200)
 122 
      real         tmp(200)
114 
      character*80 varname
 123 
      character*80 varname
115 
      real         rtime
 124 
      real         rtime
- 
 
 125 
      integer      varid
116 
      integer      is2d
 126 
      integer      cdfid
117 
      integer      plev
 127 
      integer      stat
118 
 
 - 
 
119 
c     Init the flag for 2D variables - assume a 3D field
 128 
      real         delta
120 
      is2d = 0
 129 
      integer      closear
- 
 
 130 
      real         maxps,minps
121 
 
 131 
 
122 
c     Init the flag for pressure levels (PS is not needed)
 132 
c     ---- Read data from netCDF file as they are ---------------------
123 
      plev = 0
 - 
 
124 
 
 133 
 
125 
c     Inquire dimensions and grid constants if <fid> is negative
 134 
c     Set file identifier
126 
      if (fid.lt.0) then
 135 
      if (fid.lt.0) then
- 
 
 136 
        cdfid = -fid
- 
 
 137 
      else
- 
 
 138 
        cdfid = fid
- 
 
 139 
      endif
127 
 
 140 
 
128 
c        Get grid info for the selected variable
 141 
c     Special handling if 3D pressure is
129 
         if ( fieldname.eq.'PLEV' ) then
 142 
      if ( fieldname.eq.'P' ) then
130 
            varname = 'PS'
 143 
         fieldname = 'U'
- 
 
 144 
      endif
- 
 
 145 
 
- 
 
 146 
c     Get number of dimensions of variable -> ndim
- 
 
 147 
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
- 
 
 148 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 149 
      ierr = nf90_inquire_variable(cdfid, varid, ndims  = ndim)
- 
 
 150 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 151 
c      if ( ndim.ne.4 ) then
- 
 
 152 
c          print*,' ERROR: netCDF variables need to be 4D'
- 
 
 153 
c          print*,'      ',trim(fieldname)
131 
            stagz   = 0.
 154 
c          stop
- 
 
 155 
c      endif
- 
 
 156 
 
- 
 
 157 
c     Get dimensions -> vardim(1:ndim),dimname(1:ndim)
- 
 
 158 
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
- 
 
 159 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 160 
      ierr = nf90_inquire_variable(cdfid, varid,
- 
 
 161 
     >                                   dimids = dimids(1:ndim))
- 
 
 162 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 163 
      do i=1,ndim
132 
            call getdef(-fid,varname,ndim,mdv,vardim,
 164 
           ierr = nf90_inquire_dimension(cdfid, dimids(i),
133 
     >                  varmin,varmax,stag,ierr)
 165 
     >                               name = dimname(i) )
- 
 
 166 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 167 
           ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))
- 
 
 168 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 169 
      enddo
- 
 
 170 
 
- 
 
 171 
c     Check whether the list of dimensions is OK
134 
            if (ierr.ne.0) goto 900
 172 
      if ( ndim.eq.4 ) then
135 
            call getcfn(-fid,cstfile,ierr)
 173 
        if ( ( dimname(1).ne.'lon'  ).or.
136 
            if (ierr.ne.0) goto 903
 174 
     >     ( dimname(2).ne.'lat'  ).or.
137 
            call cdfopn(cstfile,cstid,ierr)
 175 
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
138 
            if (ierr.ne.0) goto 903
 176 
     >     ( dimname(4).ne.'time' ) )
- 
 
 177 
     >  then
- 
 
 178 
         print*,' ERROR: the dimensions of the variable are not correct'
- 
 
 179 
         print*,'        expected -> lon / lat / lev / time'
139 
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
 180 
         print*, ( trim(dimname(i))//' / ',i=1,ndim )
- 
 
 181 
         stop
- 
 
 182 
        endif
- 
 
 183 
      else
140 
            if (ierr.ne.0) goto 903
 184 
        if ( ( dimname(1).ne.'lon'  ).or.
141 
            call clscdf(cstid,ierr)
 185 
     >       ( dimname(2).ne.'lat'  ).or.
142 
            if (ierr.ne.0) goto 903
 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
143
194 
 
144 
         elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then
- 
 
145 
            varname = 'PS'
 - 
 
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)
 195 
c     Check about the type of vertical levels
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
 196 
      if ( ndim.eq.4 ) then
161 
            call getdef(-fid,varname,ndim,mdv,vardim,
 - 
 
162 
     >                  varmin,varmax,stag,ierr)
 - 
 
163 
            if (ierr.ne.0) goto 900
 - 
 
164
197 
 
- 
 
 198 
c         Get dimension of AK,BK
165 
         endif
 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)
166 
 
 203 
 
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 
 
 204 
 
- 
 
 205 
          varname=dimname(3)
178 
c        Get pole position - if no constants file available, set default pole
 206 
          ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 207 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
179 
         call getcfn(-fid,cstfile,ierr)
 208 
          ierr = nf90_get_att(cdfid, varid, "standard_name", leveltype)
180 
         if (ierr.eq.0) then
209 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
181 
            call cdfopn(cstfile,cstid,ierr)
 210 
          if ( (leveltype.ne.'hybrid_sigma_pressure').and.
182 
            if (ierr.ne.0) goto 903
 211 
     >         (leveltype.ne.'air_pressure'         ) )
- 
 
 212 
     >    then
183 
            call getpole(cstid,pollon,pollat,ierr)
 213 
             print*,' ERROR: input netCDF data must be on hybrid-sigma'
184 
            if (ierr.ne.0) goto 903
 214 
             print*,'        or air pressure levels!',trim(leveltype)
185 
            call clscdf(cstid,ierr)
 215 
             stop
186 
            if (ierr.ne.0) goto 903
 216 
          endif
187 
         else
 217 
      else
188 
            pollon = 0.
 - 
 
189 
            pollat = 90.
 218 
          leveltype='air_pressure'
190 
         endif
 219 
      endif
191 
 
 220 
 
- 
 
 221 
c     Allocate memory for reading arrays
- 
 
 222 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
- 
 
 223 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
- 
 
 224 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
- 
 
 225 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
- 
 
 226 
      allocate(lon(vardim(1)),stat=stat)
- 
 
 227 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
- 
 
 228 
      allocate(lat(vardim(2)),stat=stat)
- 
 
 229 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
- 
 
 230 
      allocate(lev(vardim(3)),stat=stat)
- 
 
 231 
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
- 
 
 232 
      allocate(times(vardim(4)),stat=stat)
- 
 
 233 
      if (stat.ne.0) print*,'*** error allocating array times   ***'
- 
 
 234 
      if ( ndim.eq.4 ) then
- 
 
 235 
        allocate(aktmp(nakbktmp),stat=stat)
- 
 
 236 
        if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
- 
 
 237 
        allocate(bktmp(nakbktmp),stat=stat)
- 
 
 238 
        if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
- 
 
 239 
      endif
- 
 
 240 
 
- 
 
 241 
c     Get domain longitudes, latitudes and levels
- 
 
 242 
      varname = dimname(1)
- 
 
 243 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 244 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 245 
      ierr = nf90_get_var(cdfid,varid,lon)
- 
 
 246 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 247 
      varname = dimname(2)
- 
 
 248 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 249 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 250 
      ierr = nf90_get_var(cdfid,varid,lat)
- 
 
 251 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 252 
      varname = dimname(3)
- 
 
 253 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 254 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 255 
      ierr = nf90_get_var(cdfid,varid,lev)
- 
 
 256 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 257
- 
 
 258 
c     Get ak and bk
- 
 
 259 
      if ( ndim.eq.4 ) then
- 
 
 260 
 
- 
 
 261 
c       Read ak, bk
- 
 
 262 
        varname='hyam'
- 
 
 263 
        ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 264 
        IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 265 
        ierr = nf90_get_var(cdfid,varid,aktmp)
- 
 
 266 
        IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 267 
        varname='hybm'
- 
 
 268 
        ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 269 
        IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 270 
        ierr = nf90_get_var(cdfid,varid,bktmp)
- 
 
 271 
        IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 272 
 
- 
 
 273 
c       Check that unit of ak is in hPa - if necessary correct it
- 
 
 274 
        varname='hyam'
- 
 
 275 
        ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 276 
        IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 277 
        ierr = nf90_get_att(cdfid, varid, "units", units)
- 
 
 278 
        IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 279 
        if ( units.eq.'Pa' ) then
- 
 
 280 
           do k=1,nakbktmp
- 
 
 281 
              aktmp(k) = 0.01 * aktmp(k)
- 
 
 282 
           enddo
- 
 
 283 
        endif
- 
 
 284 
      endif
- 
 
 285 
 
- 
 
 286 
c     Check that unit of lev is in hPa - if necessary correct it
- 
 
 287 
      if ( leveltype.eq.'air_pressure' .and. vardim(3) .gt. 1) then
- 
 
 288 
         varname='lev'
- 
 
 289 
         ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 290 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 291 
         ierr = nf90_get_att(cdfid, varid, "units", units)
- 
 
 292 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 293 
         if ( units.eq.'Pa' ) then
- 
 
 294 
            do k=1,vardim(3)
- 
 
 295 
               lev(k) = 0.01 * lev(k)
- 
 
 296 
            enddo
- 
 
 297 
         endif
- 
 
 298 
      endif
- 
 
 299 
 
192 
c     Get non-constant grid parameters (surface pressure and vertical grid)
 300 
c     Decide whether to swap vertical levels - highest pressure at index 1
- 
 
 301 
      vertical_swap = 1
- 
 
 302 
      if ( leveltype.eq.'hybrid_sigma_pressure') then
- 
 
 303 
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
- 
 
 304 
     >       (aktmp(2) + bktmp(2) * 1000.) )
- 
 
 305 
     >  then
- 
 
 306 
          vertical_swap = 0
- 
 
 307 
        endif
- 
 
 308 
      elseif ( leveltype.eq.'air_pressure') then
- 
 
 309 
        if ( lev(1).gt.lev(2) ) then
- 
 
 310 
          vertical_swap = 0
- 
 
 311 
        endif
193 
      else
 312 
      endif
- 
 
 313 
c      print*,' Vertical SWAP P -> ', vertical_swap
194
314 
 
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
 315 
c     Get time information (check if time is correct)
198 
c        used for model-level (2D) trajectories
 316 
      varname = 'time'
199 
         if ( fieldname.eq.'P.ML' ) then
 317 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
200 
 
 - 
 
201 
c           Get the right time to read
 318 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
202 
            call gettimes(fid,times,ntimes,ierr)
 319 
      ierr = nf90_get_var(cdfid,varid,times)
203 
            if (ierr.ne.0) goto 901
 320 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
204 
            isok=0
 321 
      isok=0
205 
            do i=1,ntimes
 322 
      do i=1,vardim(4)
206 
               if (abs(time-times(i)).lt.eps) then
 323 
        if (abs(time-times(i)).lt.eps) then
207 
                  isok = 1
 324 
               isok = 1
208 
                  rtime = times(i)
 325 
               rtime = times(i)
209 
               elseif (timecheck.eq.'no') then
 326 
        elseif (timecheck.eq.'no') then
210 
                  isok = 1
 327 
               isok = 1
211 
                  rtime = times(1)
 328 
               rtime = times(1)
212 
               endif
 329 
        endif
213 
            enddo
 330 
      enddo
- 
 
 331 
      if ( isok.eq.0 ) then
- 
 
 332 
         print*,' ERROR: time ',rtime,' not found on netCDF file'
- 
 
 333 
         stop
- 
 
 334 
      endif
214 
 
 335 
 
215 
c           Read surface pressure and 3D pressure
 336 
c     Read surface pressure
216 
            varname='PS'
 337 
      varname='PS'
217 
            call getdat(fid,varname,rtime,0,ps,ierr)
 338 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
218 
            if (ierr.ne.0) goto 904
 339 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
219 
            varname='P'
 340 
      ierr = nf90_get_var(cdfid,varid,tmp2)
220 
            call getdat(fid,varname,rtime,0,p3,ierr)
 341 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 342
221 
            if (ierr.ne.0) goto 904
 343 
c     Check that surface pressure is in hPa
222 
 
 - 
 
223 
c           Set MDV to 1050. - otherwise interpolation routines don't work
 344 
      if ( leveltype.eq.'hybrid_sigma_pressure' ) then
224 
            do i=1,nx
 345 
        maxps = -1.e19
225 
              do j=1,ny
 346 
        minps =  1.e19
- 
 
 347 
        do i=1,vardim(1)
226 
                do k=1,nz
 348 
          do j=1,vardim(2)
227 
                   if ( p3(i,j,k).lt.0. ) p3(i,j,k) = 1050.
 349 
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
228 
                enddo
 350 
             if (tmp2(i,j).lt.minps) minps = tmp2(i,j)
229 
              enddo
 351 
          enddo
230 
            enddo
 352 
        enddo
231 
 
 - 
 
- 
 
 353 
        if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then
232 
c           Don't care about other stuff - finish subroutine
 354 
           print*,' ERROR: surface pressre PS must be in hPa'
- 
 
 355 
           print*,'       ',maxps,minps
233 
            goto 800
 356 
           stop
234 
 
 357 
        endif
235 
         endif
 358 
      endif
236 
 
 359 
 
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)
 - 
 
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)
 360 
c     ---- Define output of subroutine --------------------------------
249 
            if (ierr.ne.0) goto 903
 - 
 
250 
            call clscdf(cstid,ierr)
 - 
 
251 
            if (ierr.ne.0) goto 903
 - 
 
252
361 
 
253 
         elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then
- 
 
254 
            varname = 'PS'
 - 
 
255 
            stagz   = -0.5
 - 
 
256 
            call getdef(fid,varname,ndim,mdv,vardim,
 362 
c     If not full list of vertical levels, reduce AK,BK arrays
257 
     >                  varmin,varmax,stag,ierr)
 363 
      if ( (leveltype.eq.'hybrid_sigma_pressure').and.
258 
            if (ierr.ne.0) goto 900
 364 
     >     (nakbktmp.ne.vardim(3) ) )
259 
            call getcfn(fid,cstfile,ierr)
 - 
 
260 
            if (ierr.ne.0) goto 903
 365 
     >then
261 
            call cdfopn(cstfile,cstid,ierr)
 366 
c         print*,' WARNING: only subset of vertical levels used...'
262 
            if (ierr.ne.0) goto 903
 367 
         do k=1,vardim(3)
263 
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
 - 
 
264 
            if (ierr.ne.0) goto 903
 368 
            if ( vertical_swap.eq.1 ) then
265 
            call clscdf(cstid,ierr)
 369 
               aktmp(k) = aktmp( k+nakbktmp-vardim(3) )
266 
            if (ierr.ne.0) goto 903
 370 
               bktmp(k) = bktmp( k+nakbktmp-vardim(3) )
267 
 
 - 
 
268 
         else
371 
            endif
269 
            varname=fieldname
 372 
         enddo
270 
            call getdef(fid,varname,ndim,mdv,vardim,
 - 
 
271 
     >                  varmin,varmax,stag,ierr)
 - 
 
272 
            if (ierr.ne.0) goto 900
 - 
 
273 
            if (vardim(3).eq.1) is2d = 1
 - 
 
274 
         endif
 373 
      endif
275 
 
 374 
 
276 
c        Get time information (check if time is correct)
 375 
c     Set the grid dimensions and constants
277 
         call gettimes(fid,times,ntimes,ierr)
 376 
      nx      = vardim(1)
- 
 
 377 
      ny      = vardim(2)
278 
         if (ierr.ne.0) goto 901
 378 
      nz      = vardim(3)
279 
         isok=0
 379 
      xmin    = lon(1)
- 
 
 380 
      ymin    = lat(1)
- 
 
 381 
      xmax    = lon(nx)
280 
         do i=1,ntimes
 382 
      ymax    = lat(ny)
- 
 
 383 
      dx      = (xmax-xmin)/real(nx-1)
281 
            if (abs(time-times(i)).lt.eps) then
 384 
      dy      = (ymax-ymin)/real(ny-1)
- 
 
 385 
      pollon  = 0.
- 
 
 386 
      pollat  = 90.
282 
               isok = 1
 387 
      stagz   = 0.
283 
               rtime = times(i)
 388 
      delta   = xmax-xmin-360.
284 
            elseif (timecheck.eq.'no') then
 389 
      if (abs(delta+dx).lt.eps) then
285 
               isok = 1
 390 
          xmax    = xmax + dx
286 
               rtime = times(1)
 391 
          nx      = nx + 1
- 
 
 392 
          closear = 1
- 
 
 393 
      else
- 
 
 394 
          closear = 0
287 
            endif
 395 
      endif
288 
         enddo
 - 
 
289 
         if ( isok.eq.0) goto 905
 - 
 
290 
 
 396 
 
291 
c        If 2D tracing requested: take dummay values for PS, AKLEV,BKLEV,AKLAY,BKLAY
 397 
c     Save the output arrays (if fid>0) - close arrays on request
292 
         if ( is2d.eq.1 ) then
 398 
      if ( fid.gt.0 ) then
293
399 
 
- 
 
 400 
c        Calculate layer pressures
- 
 
 401 
         if (leveltype.eq.'hybrid_sigma_pressure' ) then
294 
            do i=1,nx
 402 
            do i=1,vardim(1)
295 
               do j=1,ny
 403 
              do j=1,vardim(2)
296 
                  ps(i,j) = 1050.
 404 
                 do k=1,vardim(3)
- 
 
 405 
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
297 
               enddo
 406 
                 enddo
298 
            enddo
 407 
              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
 408 
           enddo
306 
 
 - 
 
307 
c        3D tracing - read PS, AKLEV,BKLEV,AKLAY;BKLAY
 - 
 
308 
         else
 - 
 
309 
 
 - 
 
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
 409 
         elseif (leveltype.eq.'air_pressure' ) then
321 
            plev = 1
 410 
           do i=1,vardim(1)
322 
            do i=1,nz
 411 
              do j=1,vardim(2)
323 
              if ( (abs(stagz).lt.eps).and.(abs(bklev(i)).gt.eps) ) then
 - 
 
324 
                plev = 0
 412 
                 do k=1,vardim(3)
325 
              endif
 - 
 
326 
              if ( (abs(stagz).gt.eps).and.(abs(bklay(i)).gt.eps) ) then
 - 
 
327 
                plev = 0
 413 
                  tmp3(i,j,k)=lev(k)
328 
              endif
 - 
 
329 
            enddo
 414 
                 enddo
330 
 
 - 
 
331 
c           Read surface pressure if needed
 - 
 
332 
            if ( plev.ne.1 ) then
 - 
 
333 
              varname='PS'
 - 
 
334 
              call getdat(fid,varname,rtime,0,ps,ierr)
 - 
 
335 
              if (ierr.ne.0) goto 904
 - 
 
336 
            else
 - 
 
337 
              do i=1,nx
 - 
 
338 
                do j=1,ny
 - 
 
339 
                    ps(i,j) = 1000.
 - 
 
340 
                enddo
 415 
              enddo
341 
              enddo
 416 
           enddo
342 
            endif
 417 
         endif
343 
 
 418 
 
- 
 
 419 
c        Get PS - close array on demand
- 
 
 420 
         do j=1,vardim(2)
- 
 
 421 
           do i=1,vardim(1)
- 
 
 422 
             ps(i,j) = tmp2(i,j)
- 
 
 423 
           enddo
- 
 
 424 
           if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)
344 
         endif
 425 
         enddo
345 
 
 426 
 
346 
c        Calculate layer and level pressures
 427 
c        Get P3 - close array on demand + vertical swap
347 
         do i=1,nx
 428 
         do j=1,vardim(2)
348 
            do j=1,ny
 429 
           do k=1,vardim(3)
349 
               do k=1,nz
 430 
             do i=1,vardim(1)
350 
                  if ( abs(stagz).lt.eps ) then
 431 
               if ( vertical_swap.eq.1 ) then
351 
                     p3(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
 432 
                  p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)
352 
                  else
 433 
               else
353 
                     p3(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
 434 
                  p3(i,j,k) = tmp3(i,j,k)
354 
                  endif
 435 
               endif
355 
               enddo
 436 
             enddo
- 
 
 437 
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
356 
            enddo
 438 
           enddo
357 
         enddo
 439 
         enddo
358 
 
 440 
 
359 
c        Set the ak and bk for the vertical grid
 441 
c        Get AK,BK - vertical swap on demand
- 
 
 442 
         if ( leveltype.eq.'hybrid_sigma_pressure' ) then
360 
         do k=1,nz
 443 
           do k=1,vardim(3)
361 
            if ( abs(stagz).lt.eps ) then
 444 
              if ( vertical_swap.eq.1 ) then
362 
               ak(k)=aklev(k)
 445 
                 ak(k) = aktmp(vardim(3)-k+1)
363 
               bk(k)=bklev(k)
 446 
                 bk(k) = bktmp(vardim(3)-k+1)
364 
            else
 447 
              else
365 
               ak(k)=aklay(k)
 448 
                 ak(k) = aktmp(k)
366 
               bk(k)=bklay(k)
 449 
                 bk(k) = bktmp(k)
367 
            endif
 450 
              endif
368 
         enddo
 451 
           enddo
- 
 
 452 
         elseif (leveltype.eq.'air_pressure' ) then
- 
 
 453 
           do k=1,vardim(3)
- 
 
 454 
              if ( vertical_swap.eq.1 ) then
- 
 
 455 
                 ak(k) = lev(vardim(3)-k+1)
- 
 
 456 
                 bk(k) = 0.
- 
 
 457 
              else
- 
 
 458 
                ak(k) = lev(k)
- 
 
 459 
                bk(k) = 0.
- 
 
 460 
              endif
- 
 
 461 
           enddo
- 
 
 462 
         endif
369 
 
 463 
 
370 
      endif
 464 
      endif
371
465 
 
372 
c     Exit point for subroutine
 - 
 
373 
 800  continue
 - 
 
374 
      return
 - 
 
375
466 
 
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
 467 
      return
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
468
392 
      end
 469 
      end
393 
 
 470 
 
394 
c     ------------------------------------------------------------
 471 
c     ------------------------------------------------------------
395 
c     Read wind information
 472 
c     Read wind information
Line 404... Line 481...
404 
c     information is provided in <stagz> and gives the reference to either
 481 
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
 482 
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,
 483 
c     performed to have an agreement with the grid specified by <xmin,xmax,
407 
c     ymin,ymax,dx,dy,nx,ny,nz>.
 484 
c     ymin,ymax,dx,dy,nx,ny,nz>.
408 
 
 485 
 
- 
 
 486 
      use netcdf
- 
 
 487 
 
409 
      implicit none
 488 
      implicit none
410 
 
 489 
 
411 
c     Declaration of variables and parameters
 490 
c     Declaration of variables and parameters
412 
      integer      fid                 ! File identifier
 491 
      integer      fid                 ! File identifier
413 
      character*80 fieldname           ! Name of the wind field
 492 
      character*80 fieldname           ! Name of the wind field
Line 437... Line 516...
437 
      integer      ntimes
 516 
      integer      ntimes
438 
      character*80 cstfile
 517 
      character*80 cstfile
439 
      integer      cstid
 518 
      integer      cstid
440 
      real         aklay(200),bklay(200),aklev(200),bklev(200)
 519 
      real         aklay(200),bklay(200),aklev(200),bklev(200)
441 
      real         ps(nx,ny)
 520 
      real         ps(nx,ny)
- 
 
 521 
      integer      dimids (nf90_max_var_dims)
- 
 
 522 
      character*80 dimname(nf90_max_var_dims)
- 
 
 523 
      integer      varid
- 
 
 524 
      integer      cdfid
- 
 
 525 
      real,allocatable, dimension (:)     :: lon,lat,lev
- 
 
 526 
      real,allocatable, dimension (:,:)   :: tmp2
- 
 
 527 
      real,allocatable, dimension (:,:,:) :: tmp3
- 
 
 528 
      real,allocatable, dimension (:)     :: aktmp,bktmp
- 
 
 529 
      character*80  leveltype
- 
 
 530 
      integer       vertical_swap
- 
 
 531 
      character*80  units
- 
 
 532 
      integer       nakbktmp
- 
 
 533 
      integer       dimid
442 
 
 534 
 
443 
c     Auxiliary variables
 535 
c     Auxiliary variables
444 
      integer      isok
 536 
      integer      isok
445 
      integer      i,j,k
 537 
      integer      i,j,k
446 
      integer      nz1
 538 
      integer      nz1
447 
      real         rtime
 539 
      real         rtime
- 
 
 540 
      integer      closear
- 
 
 541 
      integer      stat
- 
 
 542 
      real         delta
- 
 
 543 
      integer      pressure
- 
 
 544 
 
- 
 
 545 
c     Init mdv
- 
 
 546 
      mdv = -999.
- 
 
 547
- 
 
 548 
c     Special handling if 3D pressure is
- 
 
 549 
      if ( fieldname.eq.'P' ) then
- 
 
 550 
         fieldname = 'U'
- 
 
 551 
         pressure  = 1
- 
 
 552 
      else
- 
 
 553 
         pressure = 0
- 
 
 554 
      endif
448 
 
 555 
 
- 
 
 556 
c     Get the number of dimensions -> ndim
- 
 
 557 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
- 
 
 558 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 559 
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
- 
 
 560 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 561 
 
449 
c     Read variable definition - for P, PLEV and PLAY: load also ak,bk
 562 
c     Get the dimensions of the arrays -> varid(1:ndim)
- 
 
 563 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
- 
 
 564 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 565 
      ierr = nf90_inquire_variable(fid, varid,
- 
 
 566 
     >                                dimids = dimids(1:ndim))
- 
 
 567 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 568 
      do i=1,ndim
- 
 
 569 
        ierr = nf90_inquire_dimension(fid, dimids(i),
- 
 
 570 
     >                               name = dimname(i) )
- 
 
 571 
        IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 572 
        ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
- 
 
 573 
        IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 574 
      enddo
- 
 
 575
- 
 
 576 
 
- 
 
 577 
c     Check whether the list of dimensions is OK
- 
 
 578 
      if (ndim .eq. 4) then
450 
      if ( ( fieldname.eq.'PLEV' ).or.
 579 
          if ( ( dimname(1).ne.'lon'  ).or.
451 
     >     ( fieldname.eq.'PLAY' ).or.
 580 
     >       ( dimname(2).ne.'lat'  ).or.
- 
 
 581 
     >       ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
452 
     >     ( fieldname.eq.'P'    ) )
 582 
     >         ( dimname(4).ne.'time' ) )
453 
     >then
 583 
     >    then
454 
         call getcfn(fid,cstfile,ierr)
 - 
 
455 
         if (ierr.ne.0) goto 905
 - 
 
456 
         call cdfopn(cstfile,cstid,ierr)
 584 
        print*,' ERROR: the dimensions of the variable are not correct'
457 
         if (ierr.ne.0) goto 905
 - 
 
458 
         call getlevs(cstid,nz1,aklev,bklev,aklay,bklay,ierr)
 585 
        print*,'        expected -> lon / lat / lev / time'
459 
         if (ierr.ne.0) goto 905
 - 
 
460 
         call clscdf(cstid,ierr)
 - 
 
461 
         if (ierr.ne.0) goto 905
 586 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
462 
         varname = 'PS'
 587 
            stop
463 
         call getdef(fid,varname,ndim,mdv,vardim,
 - 
 
464 
     >               varmin,varmax,stag,ierr)
 - 
 
465 
         vardim(3) = nz1
 588 
          endif
466 
         if (ierr.ne.0) goto 900
 - 
 
467 
 
 - 
 
468 
      else
 589 
      else
469 
         varname = fieldname
 590 
          if ( ( dimname(1).ne.'lon'  ).or.
470 
         call getdef(fid,varname,ndim,mdv,vardim,
 - 
 
471 
     >               varmin,varmax,stag,ierr)
 591 
     >         ( dimname(2).ne.'lat'  ).or.
472 
         if (ierr.ne.0) goto 900
 592 
     >         ( dimname(3).ne.'time' ) )
473 
         stagz=stag(3)
 - 
 
474 
      endif
 593 
     >    then
475 
 
 - 
 
476 
c     Get time information (set time to first one in the file)
 594 
        print*,' ERROR: the dimensions of the variable are not correct'
477 
      call gettimes(fid,times,ntimes,ierr)
 595 
        print*,'        expected -> lon / lat / time'
478 
      if (ierr.ne.0) goto 902
 - 
 
479 
      isok=0
 - 
 
480 
      do i=1,ntimes
 - 
 
481 
         if (abs(time-times(i)).lt.eps) then
 596 
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
482 
            isok = 1
 597 
            stop
483 
            rtime = times(i)
 - 
 
484 
         elseif (timecheck.eq.'no') then
 - 
 
485 
            isok = 1
 598 
          endif
486 
            rtime = times(1)
 - 
 
487 
         endif
 599 
      end if
488 
      enddo
 - 
 
489 
      if ( isok.eq.0 ) goto 904
 - 
 
490 
 
 600
- 
 
 601 
c     Check about the type of vertical levels
- 
 
 602 
      if (ndim .eq. 4) then
- 
 
 603 
 
- 
 
 604 
c         Get dimension of AK,BK
- 
 
 605 
          varname = 'nhym'
- 
 
 606 
          ierr = NF90_INQ_DIMID(fid,varname,dimid)
- 
 
 607 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 608 
          ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
- 
 
 609 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 610 
 
- 
 
 611 
          varname=dimname(3)
- 
 
 612 
          ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 613 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 614 
          ierr = nf90_get_att(fid, varid, "standard_name", leveltype)
- 
 
 615 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 616 
          if ( (leveltype.ne.'hybrid_sigma_pressure').and.
- 
 
 617 
     >         (leveltype.ne.'air_pressure'         ) ) then
- 
 
 618 
             print*,' ERROR: input netCDF data must be on hybrid-sigma'
- 
 
 619 
             print*,'        or air pressure levels!',trim(leveltype)
491 
c     Read  field
 620 
             stop
- 
 
 621 
          endif
- 
 
 622 
      else
- 
 
 623 
         leveltype='air_pressure'
- 
 
 624 
      endif
- 
 
 625 
 
- 
 
 626 
c     Allocate memory for reading arrays - depending on <closear>
- 
 
 627 
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
- 
 
 628 
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
- 
 
 629 
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
- 
 
 630 
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
- 
 
 631 
      allocate(lon(vardim(1)),stat=stat)
- 
 
 632 
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
- 
 
 633 
      allocate(lat(vardim(2)),stat=stat)
- 
 
 634 
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
- 
 
 635 
      allocate(lev(vardim(3)),stat=stat)
492 
      if ( ( fieldname.eq.'P' ).or.(fieldname.eq.'PLAY') )  then       ! P, PLAY
 636 
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
- 
 
 637 
      if (ndim .eq. 4) then
- 
 
 638 
          allocate(aktmp(nakbktmp),stat=stat)
- 
 
 639 
          if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
- 
 
 640 
          allocate(bktmp(nakbktmp),stat=stat)
- 
 
 641 
          if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
- 
 
 642 
      end if
- 
 
 643 
 
- 
 
 644 
c     Get domain boundaries - longitude, latitude, levels
- 
 
 645 
      varname = dimname(1)
- 
 
 646 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 647 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 648 
      ierr = nf90_get_var(fid,varid,lon)
- 
 
 649 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 650 
      varname = dimname(2)
- 
 
 651 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 652 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 653 
      ierr = nf90_get_var(fid,varid,lat)
- 
 
 654 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 655 
      varname = dimname(3)
- 
 
 656 
      ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 657 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 658 
      ierr = nf90_get_var(fid,varid,lev)
- 
 
 659 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 660 
      varmin(1) = lon(1)
- 
 
 661 
      varmax(1) = lon( vardim(1) )
- 
 
 662 
      varmin(2) = lat(1)
- 
 
 663 
      varmax(2) = lat( vardim(2) )
- 
 
 664 
 
- 
 
 665 
c     Get ak and bk
- 
 
 666 
      if (ndim .eq. 4) then
- 
 
 667 
          varname='hyam'
- 
 
 668 
          ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 669 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 670 
          ierr = nf90_get_var(fid,varid,aktmp)
- 
 
 671 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 672 
          varname='hybm'
- 
 
 673 
          ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 674 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 675 
          ierr = nf90_get_var(fid,varid,bktmp)
- 
 
 676 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 677 
 
- 
 
 678 
c         Check that unit of ak is in hPa - if necessary correct it
- 
 
 679 
          varname='hyam'
- 
 
 680 
          ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 681 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 682 
          ierr = nf90_get_att(fid, varid, "units", units)
- 
 
 683 
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 684 
          if ( units.eq.'Pa' ) then
- 
 
 685 
             do k=1,nakbktmp
- 
 
 686 
                aktmp(k) = 0.01 * aktmp(k)
493 
         stagz   = -0.5
 687 
             enddo
- 
 
 688 
          endif
- 
 
 689 
      end if
- 
 
 690 
 
- 
 
 691 
c     Check that unit of lev is in hPa - if necessary correct it
- 
 
 692 
      if ( leveltype.eq.'air_pressure' .and. vardim(3) .gt. 1) then
494 
         varname = 'PS'
 693 
         varname='lev'
- 
 
 694 
         ierr = NF90_INQ_VARID(fid,varname,varid)
- 
 
 695 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
495 
         call getdat(fid,varname,rtime,0,ps,ierr)
 696 
         ierr = nf90_get_att(fid, varid, "units", units)
- 
 
 697 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
496 
         if (ierr.ne.0) goto 903
 698 
         if ( units.eq.'Pa' ) then
- 
 
 699 
            do k=1,vardim(3)
- 
 
 700 
               lev(k) = 0.01 * lev(k)
497 
         do i=1,nx
 701 
            enddo
- 
 
 702 
         endif
- 
 
 703 
      endif
- 
 
 704 
 
- 
 
 705 
c     Decide whether to swap vertical levels
- 
 
 706 
      vertical_swap = 1
- 
 
 707 
      if ( leveltype.eq.'hybrid_sigma_pressure') then
- 
 
 708 
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
- 
 
 709 
     >       (aktmp(2) + bktmp(2) * 1000.) )
- 
 
 710 
     >  then
- 
 
 711 
          vertical_swap = 0
- 
 
 712 
        endif
- 
 
 713 
      elseif ( leveltype.eq.'air_pressure') then
- 
 
 714 
        if ( lev(1).gt.lev(2) ) then
- 
 
 715 
          vertical_swap = 0
- 
 
 716 
        endif
- 
 
 717 
      endif
- 
 
 718 
 
- 
 
 719 
c     Read data /or set 3D pressure field)
- 
 
 720 
      if ( pressure.eq.0 ) then
- 
 
 721 
         ierr = NF90_INQ_VARID(fid,fieldname,varid)
- 
 
 722 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 723 
         ierr = nf90_get_var(fid,varid,tmp3)
- 
 
 724 
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 725 
      else
- 
 
 726 
         if (leveltype.eq.'hybrid_sigma_pressure' ) then
498 
            do j=1,ny
 727 
            do i=1,vardim(1)
- 
 
 728 
              do j=1,vardim(2)
499 
               do k=1,nz
 729 
                 do k=1,vardim(3)
500 
                  field(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
 730 
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
501 
               enddo
 731 
                 enddo
502 
            enddo
 732 
              enddo
503 
         enddo
 733 
           enddo
504
- 
 
505 
      elseif ( fieldname.eq.'PLEV' )  then                             ! PLEV
 - 
 
506 
         stagz   = 0.
 - 
 
507 
         varname = 'PS'
 - 
 
508 
         call getdat(fid,varname,rtime,0,ps,ierr)
 734 
         elseif (leveltype.eq.'air_pressure' ) then
509 
         if (ierr.ne.0) goto 903
 - 
 
510 
         do i=1,nx
 735 
           do i=1,vardim(1)
511 
            do j=1,ny
 736 
              do j=1,vardim(2)
512 
               do k=1,nz
 737 
                 do k=1,vardim(3)
513 
                  field(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
 738 
                  tmp3(i,j,k)=lev(k)
514 
               enddo
 739 
                 enddo
515 
            enddo
 740 
              enddo
516 
         enddo
 741 
           enddo
517 
 
 - 
 
518 
      else                                                             ! Other fields
 - 
 
519 
         varname=fieldname
 742 
         endif
520 
         call getdat(fid,varname,rtime,0,field,ierr)
 - 
 
521 
         if (ierr.ne.0) goto 903
 - 
 
522 
 
 - 
 
523 
      endif
 743 
      endif
524 
 
 744
525 
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
526 
      if ( vardim(3).eq.1 ) then
 746 
c      if ( vardim(3).eq.1 ) then
527 
         do i=1,nx
 747 
c         do i=1,vardim(1)
528 
            do j=1,ny
 748 
c            do j=1,vardim(2)
- 
 
 749 
c               do k=1,vardim(3)
- 
 
 750 
c                  tmp3(i,j,k) = tmp3(i,j,1)
- 
 
 751 
c               enddo
- 
 
 752 
c            enddo
- 
 
 753 
c         enddo
- 
 
 754 
c      endif
- 
 
 755 
 
- 
 
 756 
c     Decide whether to close arrays
- 
 
 757 
      delta = varmax(1)-varmin(1)-360.
- 
 
 758 
      if (abs(delta+dx).lt.eps) then
- 
 
 759 
          closear = 1
- 
 
 760 
      else
- 
 
 761 
          closear = 0
- 
 
 762 
      endif
- 
 
 763 
 
- 
 
 764 
c     Save output array - close array and swap on demand
- 
 
 765 
 
- 
 
 766 
      do j=1,vardim(2)
529 
               do k=1,nz
 767 
        do k=1,nz
- 
 
 768 
          do i=1,vardim(1)
- 
 
 769 
             if (ndim .eq. 4) then
- 
 
 770 
                 if ( vertical_swap.eq.1 ) then
- 
 
 771 
                     field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
- 
 
 772 
                 else
- 
 
 773 
                     field(i,j,k) = tmp3(i,j,k)
- 
 
 774 
                 endif
- 
 
 775 
             else
530 
                  field(i,j,k) = field(i,j,1)
 776 
                field(i,j,k) = tmp3(i,j,1)
- 
 
 777 
             endif
531 
               enddo
 778 
          enddo
- 
 
 779 
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
532 
            enddo
 780 
        enddo
533 
         enddo
 781 
      enddo
534 
      endif
 - 
 
535
782
536 
 
 - 
 
537 
 
 - 
 
538 
c     Exception handling
 783 
c     Exit point
539 
      return
 784 
      return
540
785
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
 786 
      end
555 
 
 787
- 
 
 788 
 
556 
c     ------------------------------------------------------------
 789 
c     ------------------------------------------------------------
557 
c     Close input file
 790 
c     Close input file
558 
c     ------------------------------------------------------------
 791 
c     ------------------------------------------------------------
559 
 
 792 
 
560 
      subroutine input_close(fid)
 793 
      subroutine input_close(fid)
561 
 
 794 
 
562 
c     Close the input file with file identifier <fid>.
 795 
c     Close the input file with file identifier <fid>.
563 
 
 796 
 
- 
 
 797 
      use netcdf
- 
 
 798 
 
564 
      implicit none
 799 
      implicit none
565 
 
 800 
 
566 
c     Declaration of subroutine parameters
 801 
c     Declaration of subroutine parameters
567 
      integer fid
 802 
      integer fid
568 
 
 803 
 
569 
c     Auxiliary variables
 804 
c     Auxiliary variables
570 
      integer ierr
 805 
      integer ierr
571 
 
 806 
 
572 
c     Close file
 807 
c     Close file
573 
      call clscdf(fid,ierr)
 808 
      ierr = NF90_CLOSE(fid)
- 
 
 809 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
574
810
575 
      end
 811 
      end
576
812
577 
c     ------------------------------------------------------------
 813 
c     ------------------------------------------------------------
578 
c     Get a list of variables on netCDF file
 814 
c     Get a list of variables on netCDF file
Line 580... Line 816...
580 
 
 816 
 
581 
      subroutine input_getvars(fid,vnam,nvars)
 817 
      subroutine input_getvars(fid,vnam,nvars)
582 
 
 818 
 
583 
c     List of variables on netCDF file
 819 
c     List of variables on netCDF file
584 
 
 820 
 
- 
 
 821 
      use netcdf
- 
 
 822 
 
585 
      implicit none
 823 
      implicit none
586 
 
 824 
 
587 
c     Declaration of subroutine parameters
 825 
c     Declaration of subroutine parameters
588 
      integer      fid
 826 
      integer      fid
589 
      integer      nvars
 827 
      integer      nvars
590 
      character*80 vnam(200)
 828 
      character*80 vnam(200)
591 
 
 829 
 
592 
c     Auxiliary variables
 830 
c     Auxiliary variables
593 
      integer ierr
 831 
      integer ierr
- 
 
 832 
      integer i
- 
 
 833 
      integer nDims,nGlobalAtts,unlimdimid
594
834 
 
- 
 
 835 
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
- 
 
 836 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 837
595 
c     Get list and save
838 
      do i=1,nVars
596 
      call getvars(fid,nvars,vnam,ierr)
 839 
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
- 
 
 840 
      enddo
597
841
598 
      end
 842 
      end