Subversion Repositories lagranto.ecmwf

Rev

Rev 36 | Rev 46 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

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