Subversion Repositories lagranto.um

Rev

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

Rev 3 Rev 7
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 
 
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)
 52 
 900  print*,'Cannot open input file ',trim(filename)
51 
      stop
 53 
      stop
52 
 
 54 
 
53 
      end
 55 
      end
54 
 
 56 
 
55 
 
 57 
 
56 
c     ------------------------------------------------------------
 58 
c     ------------------------------------------------------------
57 
c     Read information about the grid
 59 
c     Read information about the grid
58 
c     ------------------------------------------------------------
 60 
c     ------------------------------------------------------------
59
61
60 
      subroutine input_grid
62 
      subroutine input_grid
61 
     >                   (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
 63 
     >                   (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
62 
     >                    time,pollon,pollat,p3,ps,nz,ak,bk,stagz,
 64 
     >                    time,pollon,pollat,p3,ps,nz,ak,bk,stagz,
63 
     >                    timecheck)
 65 
     >                    timecheck)
64 
 
 66 
 
65 
c     Read grid information at <time> from file with identifier <fid>.
67 
c     Read grid information at <time> from file with identifier <fid>.
66 
c     The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>
 68 
c     The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>
67 
c     with pole position at <pollon,pollat> and grid dimension <nx,ny>.
 69 
c     with pole position at <pollon,pollat> and grid dimension <nx,ny>.
68 
c     The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either
 70 
c     The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either
69 
c     on the staggered or unstaggered grid (with <stagz> as the flag).
 71 
c     on the staggered or unstaggered grid (with <stagz> as the flag).
70 
c     The surface pressure is given in <ps(nx,ny)>. If <fid> is negative,
72 
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
73 
c     only the grid dimensions and grid parameters (xmin...pollat,nz) are
72 
c     determined and returned (this is needed for dynamical allocation of
74 
c     determined and returned (this is needed for dynamical allocation of
73 
c     memory).
 75 
c     memory).
74 
 
 76 
 
- 
 
 77 
      use netcdf
- 
 
 78 
 
75 
      implicit none
 79 
      implicit none
76 
 
 80 
 
77 
c     Declaration of subroutine parameters
81 
c     Declaration of subroutine parameters
78 
      integer      fid                 ! File identifier
 82 
      integer      fid                 ! File identifier
79 
      real         xmin,xmax,ymin,ymax ! Domain size
 83 
      real         xmin,xmax,ymin,ymax ! Domain size
80 
      real         dx,dy               ! Horizontal resolution
 84 
      real         dx,dy               ! Horizontal resolution
81 
      integer      nx,ny,nz            ! Grid dimensions
 85 
      integer      nx,ny,nz            ! Grid dimensions
82 
      real         pollon,pollat       ! Longitude and latitude of pole
 86 
      real         pollon,pollat       ! Longitude and latitude of pole
83 
      real         p3(nx,ny,nz)        ! Staggered levels
 87 
      real         p3(nx,ny,nz)        ! Staggered levels
84 
      real         ps(nx,ny)           ! Surface pressure
 88 
      real         ps(nx,ny)           ! Surface pressure
85 
      real         time                ! Time of the grid information
 89 
      real         time                ! Time of the grid information
86 
      real         ak(nz),bk(nz)       ! Ak and Bk for layers or levels
 90 
      real         ak(nz),bk(nz)       ! Ak and Bk for layers or levels
87 
      real         stagz               ! Vertical staggering (0 or -0.5)
 91 
      real         stagz               ! Vertical staggering (0 or -0.5)
88 
      character*80 fieldname           ! Variable from which to take grid info
 92 
      character*80 fieldname           ! Variable from which to take grid info
89 
      character*80 timecheck           ! Either 'yes' or 'no'
 93 
      character*80 timecheck           ! Either 'yes' or 'no'
90
94
91 
c     Numerical and physical parameters
 95 
c     Numerical and physical parameters
92 
      real          eps                 ! Numerical epsilon
 96 
      real          eps                 ! Numerical epsilon
93 
      parameter    (eps=0.001)
 97 
      parameter    (eps=0.001)
94 
 
 98 
 
95 
c     Netcdf variables
 99 
c     Netcdf variables
96 
      integer      vardim(4)
 100 
      integer      vardim(4)
97 
      real         varmin(4),varmax(4)
 101 
      real         varmin(4),varmax(4)
98 
      real         mdv
 102 
      real         mdv
99 
      real         stag(4)
 103 
      real         stag(4)
100 
      integer      ndim
 104 
      integer      ndim
101 
      character*80 cstfile
 105 
      character*80 cstfile
102 
      integer      cstid
 106 
      integer      cstid
103 
      real         times(10)
 107 
      real         times
104 
      integer      ntimes
 108 
      integer      ntimes
105 
      real         aklay(nz),bklay(nz),aklev(nz),bklev(nz)
 109 
      real         aklay(nz),bklay(nz),aklev(nz),bklev(nz)
106 
      integer      nvars
 110 
      integer      nvars
107 
      character*80 vars(100)
 111 
      character*80 vars(100)
- 
 
 112 
      integer        dimids (nf90_max_var_dims)
- 
 
 113 
      character*80   dimname(nf90_max_var_dims)
108 
 
 114 
 
109 
c     Auxiliary varaibles
 115 
c     Auxiliary varaibles
110 
      integer      ierr
116 
      integer      ierr
111 
      integer      i,j,k
 117 
      integer      i,j,k
112 
      integer      isok
 118 
      integer      isok
113 
      real         tmp(200)
 119 
      real         tmp(200)
114 
      character*80 varname
 120 
      character*80 varname
115 
      real         rtime
 121 
      real         rtime
- 
 
 122 
      integer      varid
- 
 
 123 
      integer      cdfid
116 
 
 124 
 
117 
c     Inquire dimensions and grid constants if <fid> is negative
 125 
c     Set file identifier
118 
      if (fid.lt.0) then
 126 
      if (fid.lt.0) then
- 
 
 127 
        cdfid = -fid
- 
 
 128 
      else
- 
 
 129 
        cdfid = fid
- 
 
 130 
      endif
- 
 
 131
- 
 
 132 
c     Get varid
- 
 
 133 
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
- 
 
 134 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 135
- 
 
 136 
c     Get numebr of dimensions -> ndim
- 
 
 137 
      ierr = nf90_inquire_variable(cdfid, varid, ndims  = ndim)
- 
 
 138 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 139
- 
 
 140 
c     Get dim IDs
- 
 
 141 
      ierr = nf90_inquire_variable(cdfid, varid,
- 
 
 142 
     >                                   dimids = dimids(1:ndim))
- 
 
 143 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 144 
 
- 
 
 145 
c     Get dimensions -> vardim(1:ndim)
- 
 
 146 
      do i=1,ndim
- 
 
 147 
           ierr = nf90_inquire_dimension(cdfid, dimids(i),
- 
 
 148 
     >                               name = dimname(i) )
- 
 
 149 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 150 
           ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))
- 
 
 151 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 152 
      enddo
- 
 
 153
- 
 
 154 
c     Get domain min -> varmin(1:3)
- 
 
 155 
      ierr  = nf90_get_att(cdfid, varid, "xmin", varmin(1))
- 
 
 156 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 157 
      ierr  = nf90_get_att(cdfid, varid, "ymin", varmin(2))
- 
 
 158 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 159 
      ierr  = nf90_get_att(cdfid, varid, "zmin", varmin(3))
- 
 
 160 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 161
- 
 
 162 
c     Get domain max -> varmax(1:3)
- 
 
 163 
      ierr  = nf90_get_att(cdfid, varid, "xmax", varmax(1))
- 
 
 164 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 165 
      ierr  = nf90_get_att(cdfid, varid, "ymax", varmax(2))
- 
 
 166 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 167 
      ierr  = nf90_get_att(cdfid, varid, "zmax", varmax(3))
- 
 
 168 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 169
- 
 
 170 
c     Get vertical staggering -> stagz
- 
 
 171 
      ierr  = nf90_get_att(cdfid, varid, "zstag", stagz)
- 
 
 172 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 173
- 
 
 174 
c     Get missing data value -> mdv
- 
 
 175 
      ierr  = nf90_get_att(cdfid, varid, "missing_value", mdv)
- 
 
 176 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 177 
 
- 
 
 178 
c     Set the grid dimensions and constants
- 
 
 179 
      nx   = vardim(1)
- 
 
 180 
      ny   = vardim(2)
- 
 
 181 
      nz   = vardim(3)
- 
 
 182 
      xmin = varmin(1)
- 
 
 183 
      ymin = varmin(2)
- 
 
 184 
      xmax = varmax(1)
- 
 
 185 
      ymax = varmax(2)
- 
 
 186 
      dx   = (xmax-xmin)/real(nx-1)
- 
 
 187 
      dy   = (ymax-ymin)/real(ny-1)
119 
 
 188 
 
120 
c        Get grid info for the selected variable
 - 
 
121 
         if ( fieldname.eq.'PLEV' ) then
 - 
 
122 
            varname = 'PS'
 - 
 
123 
            stagz   = 0.
 - 
 
124 
            call getdef(-fid,varname,ndim,mdv,vardim,
 - 
 
125 
     >                  varmin,varmax,stag,ierr)
 - 
 
126 
            if (ierr.ne.0) goto 900
 - 
 
127 
            call getcfn(-fid,cstfile,ierr)
 - 
 
128 
            if (ierr.ne.0) goto 903
 - 
 
129 
            call cdfopn(cstfile,cstid,ierr)
 - 
 
130 
            if (ierr.ne.0) goto 903
 - 
 
131 
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
 - 
 
132 
            if (ierr.ne.0) goto 903
 - 
 
133 
            call clscdf(cstid,ierr)
 - 
 
134 
            if (ierr.ne.0) goto 903
 - 
 
135
- 
 
136 
         elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then
- 
 
137 
            varname = 'PS'
 - 
 
138 
            stagz   = -0.5
 - 
 
139 
            call getdef(-fid,varname,ndim,mdv,vardim,
 - 
 
140 
     >                  varmin,varmax,stag,ierr)
 - 
 
141 
            if (ierr.ne.0) goto 900
 - 
 
142 
            call getcfn(-fid,cstfile,ierr)
 - 
 
143 
            if (ierr.ne.0) goto 903
 - 
 
144 
            call cdfopn(cstfile,cstid,ierr)
 - 
 
145 
            if (ierr.ne.0) goto 903
 - 
 
146 
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
 - 
 
147 
            if (ierr.ne.0) goto 903
 - 
 
148 
            call clscdf(cstid,ierr)
 - 
 
149 
            if (ierr.ne.0) goto 903
 - 
 
150 
 
 - 
 
151 
         else
 - 
 
152 
            varname = fieldname
 - 
 
153 
            call getdef(-fid,varname,ndim,mdv,vardim,
 - 
 
154 
     >                  varmin,varmax,stag,ierr)
 - 
 
155 
            if (ierr.ne.0) goto 900
 - 
 
156
- 
 
157 
         endif
 - 
 
158 
 
 - 
 
159 
c        Set the grid dimensions and constants
 - 
 
160 
         nx   = vardim(1)
 - 
 
161 
         ny   = vardim(2)
 - 
 
162 
         nz   = vardim(3)
 - 
 
163 
         xmin = varmin(1)
 - 
 
164 
         ymin = varmin(2)
 - 
 
165 
         xmax = varmax(1)
 - 
 
166 
         ymax = varmax(2)
 - 
 
167 
         dx   = (xmax-xmin)/real(nx-1)
 - 
 
168 
         dy   = (ymax-ymin)/real(ny-1)
 - 
 
169 
 
 - 
 
170 
c        We want the longitudes within -180 ... + 180
 189 
c     We want the longitudes within -180 ... + 180
171 
         if ( xmin.lt.-180.) then
 190 
      if ( xmin.lt.-180.) then
172 
            xmin = xmin + 360.
 191 
            xmin = xmin + 360.
173 
            xmax = xmax + 360.
 192 
            xmax = xmax + 360.
174 
         else if ( xmax.gt.360. ) then
 193 
      else if ( xmax.gt.360. ) then
175 
            xmin = xmin - 360.
 194 
            xmin = xmin - 360.
176 
            xmax = xmax - 360.
 195 
            xmax = xmax - 360.
177 
         endif
 196 
      endif
178 
 
 - 
 
179 
c        Get pole position and number of vertical levels
 - 
 
180 
         call getcfn(-fid,cstfile,ierr)
 - 
 
181 
         if (ierr.ne.0) goto 903
 - 
 
182 
         call cdfopn(cstfile,cstid,ierr)
 - 
 
183 
         if (ierr.ne.0) goto 903
 - 
 
184 
         call getpole(cstid,pollon,pollat,ierr)
 - 
 
185 
         if (ierr.ne.0) goto 903
 - 
 
186 
         call clscdf(cstid,ierr)
 - 
 
187 
         if (ierr.ne.0) goto 903
 - 
 
188 
 
 197 
 
189 
c     Get non-constant grid parameters (surface pressure and vertical grid)
 198 
c     Get name of constants file -> cstfile
190 
      else
 199 
      ierr  = nf90_get_att(cdfid,nf90_global,
- 
 
 200 
     >                    "constants_file_name ", cstfile )
191
201 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 202 
 
192 
c        Get full grid info - in particular staggering flag
 203 
c     Get pole position from constants file
193 
         if ( fieldname.eq.'PLEV' ) then
 204 
      ierr = NF90_OPEN(TRIM(cstfile),nf90_nowrite, cstid)
194 
            varname = 'PS'
 205 
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
195 
            stagz   = 0.
 206 
      varname='pollon'
196 
            call getdef(fid,varname,ndim,mdv,vardim,
 207 
      ierr = NF90_INQ_VARID(cstid,varname,varid)
197 
     >                  varmin,varmax,stag,ierr)
 208 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
198 
            if (ierr.ne.0) goto 900
 209 
      ierr = nf90_get_var(cstid,varid,pollon)
199 
            call getcfn(fid,cstfile,ierr)
 210 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
200 
            if (ierr.ne.0) goto 903
 211 
      varname='pollat'
201 
            call cdfopn(cstfile,cstid,ierr)
 212 
      ierr = NF90_INQ_VARID(cstid,varname,varid)
202 
            if (ierr.ne.0) goto 903
 213 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
203 
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
 214 
      ierr = nf90_get_var(cstid,varid,pollat)
- 
 
 215 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
204 
            if (ierr.ne.0) goto 903
 216 
      ierr = NF90_CLOSE(cstid)
- 
 
 217 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 218
205 
            call clscdf(cstid,ierr)
 219 
c     Skip the rest if fid < 0:
206 
            if (ierr.ne.0) goto 903
 220 
      if ( fid.lt.0 ) goto 900
207
221
- 
 
 222 
c     Get ak and bk from constants file
208 
         elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then
223 
      ierr = NF90_OPEN(TRIM(cstfile),nf90_nowrite, cstid)
- 
 
 224 
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
209 
            varname = 'PS'
 225 
      varname='aklev'
210 
            stagz   = -0.5
 226 
      ierr = NF90_INQ_VARID(cstid,varname,varid)
- 
 
 227 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
211 
            call getdef(fid,varname,ndim,mdv,vardim,
 228 
      ierr = nf90_get_var(cstid,varid,aklev)
212 
     >                  varmin,varmax,stag,ierr)
 229 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
213 
            if (ierr.ne.0) goto 900
 230 
      varname='bklev'
214 
            call getcfn(fid,cstfile,ierr)
 231 
      ierr = NF90_INQ_VARID(cstid,varname,varid)
215 
            if (ierr.ne.0) goto 903
 232 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
216 
            call cdfopn(cstfile,cstid,ierr)
 233 
      ierr = nf90_get_var(cstid,varid,bklev)
- 
 
 234 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
217 
            if (ierr.ne.0) goto 903
 235 
      varname='aklay'
218 
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
 236 
      ierr = NF90_INQ_VARID(cstid,varname,varid)
219 
            if (ierr.ne.0) goto 903
 237 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
220 
            call clscdf(cstid,ierr)
 238 
      ierr = nf90_get_var(cstid,varid,aklay)
221 
            if (ierr.ne.0) goto 903
 239 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
222 
 
 - 
 
223 
         else
 240 
      varname='bklay'
224 
            varname=fieldname
 241 
      ierr = NF90_INQ_VARID(cstid,varname,varid)
- 
 
 242 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
225 
            call getdef(fid,varname,ndim,mdv,vardim,
 243 
      ierr = nf90_get_var(cstid,varid,bklay)
226 
     >                  varmin,varmax,stag,ierr)
 244 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
227 
            if (ierr.ne.0) goto 900
 245 
      ierr = NF90_CLOSE(cstid)
- 
 
 246 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
228
247
- 
 
 248 
c     print*,'AKLEV : '
- 
 
 249 
c     print*, (aklev(i),i=1,20)
- 
 
 250 
c     print*,'BKLEV : '
- 
 
 251 
c     print*,(bklev(i),i=1,20)
- 
 
 252 
c     print*,'AKLAY : '
- 
 
 253 
c     print*,(aklay(i),i=1,20)
229 
         endif
 254 
c     print*,'BKLAY : '
- 
 
 255 
c     print*,(bklay(i),i=1,20)
230 
 
 256
231 
c        Get time information (check if time is correct)
 257 
c     Get time information (check if time is correct)
- 
 
 258 
      varname = 'time'
232 
         call gettimes(fid,times,ntimes,ierr)
 259 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 260 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
233 
         if (ierr.ne.0) goto 901
 261 
      ierr = nf90_get_var(cdfid,varid,times)
- 
 
 262 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
234 
         isok=0
 263 
      isok=0
235 
         do i=1,ntimes
 - 
 
236 
            if (abs(time-times(i)).lt.eps) then
 264 
      if (abs(time-times).lt.eps) then
237 
               isok = 1
 265 
               isok = 1
238 
               rtime = times(i)
 266 
               rtime = times
239 
            elseif (timecheck.eq.'no') then
 267 
      elseif (timecheck.eq.'no') then
240 
               isok = 1
 268 
               isok = 1
241 
               rtime = times(1)
 269 
               rtime = times
242 
            endif
 - 
 
243 
         enddo
 270 
      endif
244 
         if ( isok.eq.0) goto 905
 271 
      if ( isok.eq.0 ) then
245 
 
 - 
 
246 
c        Read the level coefficients from the constants file
 272 
         print*,' ERROR: time ',rtime,' not found on netCDF file'
247 
         call getcfn(fid,cstfile,ierr)
 - 
 
248 
         if (ierr.ne.0) goto 903
 - 
 
249 
         call cdfopn(cstfile,cstid,ierr)
 - 
 
250 
         if (ierr.ne.0) goto 903
 - 
 
251 
         call getlevs(cstid,vardim(3),aklev,bklev,aklay,bklay,ierr)
 - 
 
252 
         if (ierr.ne.0) goto 903
 - 
 
253 
         call clscdf(cstid,ierr)
 - 
 
254 
         if (ierr.ne.0) goto 903
 - 
 
255 
 
 - 
 
256 
c        Set the grid dimensions and constants
 - 
 
257 
         nx   = vardim(1)
 - 
 
258 
         ny   = vardim(2)
 - 
 
259 
         nz   = vardim(3)
 - 
 
260 
         xmin = varmin(1)
 - 
 
261 
         ymin = varmin(2)
 - 
 
262 
         xmax = varmax(1)
 - 
 
263 
         ymax = varmax(2)
 273 
         stop
264 
         dx   = (xmax-xmin)/real(nx-1)
 - 
 
265 
         dy   = (ymax-ymin)/real(ny-1)
 - 
 
266 
 
 - 
 
267 
c        We want the longitudes within -180 ... + 180
 - 
 
268 
         if ( xmin.lt.-180.) then
 - 
 
269 
            xmin = xmin + 360.
 - 
 
270 
            xmax = xmax + 360.
 - 
 
271 
         else if ( xmax.gt.360. ) then
 - 
 
272 
            xmin = xmin - 360.
 - 
 
273 
            xmax = xmax - 360.
 - 
 
274 
         endif
 274 
      endif
275 
 
 - 
 
276 
c        Read surface pressure
 - 
 
277 
         varname='PS'
 275
278 
         call getdat(fid,varname,rtime,0,ps,ierr)
 - 
 
279 
         if (ierr.ne.0) goto 904
 276 
c     print*,'TIME : ',rtime
280 
 
 277 
 
- 
 
 278 
c     Read surface pressure
- 
 
 279 
      varname='PS'
- 
 
 280 
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
- 
 
 281 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 282 
      ierr = nf90_get_var(cdfid,varid,ps)
- 
 
 283 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 284
281 
c        Calculate layer and level pressures
 285 
c     Calculate layer and level pressures
282 
         do i=1,nx
 286 
      do i=1,nx
283 
            do j=1,ny
 287 
         do j=1,ny
284 
               do k=1,nz
 288 
               do k=1,nz
285 
                  if ( abs(stagz).lt.eps ) then
 289 
                  if ( abs(stagz).lt.eps ) then
286 
                     p3(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
 290 
                     p3(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
287 
                  else
 291 
                  else
288 
                     p3(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
 292 
                     p3(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
289 
                  endif
 293 
                  endif
290 
               enddo
 294 
               enddo
291 
            enddo
 - 
 
292 
         enddo
 295 
         enddo
- 
 
 296 
      enddo
293 
 
 297 
 
294 
c        Set the ak and bk for the vertical grid
 298 
c     Set the ak and bk for the vertical grid
295 
         do k=1,nz
 299 
      do k=1,nz
296 
            if ( abs(stagz).lt.eps ) then
 300 
            if ( abs(stagz).lt.eps ) then
297 
               ak(k)=aklev(k)
 301 
               ak(k)=aklev(k)
298 
               bk(k)=bklev(k)
 302 
               bk(k)=bklev(k)
299 
            else
 303 
            else
300 
               ak(k)=aklay(k)
 304 
               ak(k)=aklay(k)
301 
               bk(k)=bklay(k)
 305 
               bk(k)=bklay(k)
302 
            endif
 306 
            endif
303 
         enddo
 - 
 
304 
 
 - 
 
305 
      endif
 307 
      enddo
306
308
307 
c     Exception handling
 309 
c     Exit point
- 
 
 310 
 900  continue
- 
 
 311
308 
      return
 312 
      return
309
313
310 
 900  print*,'Cannot retrieve grid dimension from ',fid
 - 
 
311 
      stop
 - 
 
312 
 901  print*,'Cannot retrieve grid parameters from ',fid
 - 
 
313 
      stop
 - 
 
314 
 902  print*,'Grid inconsistency detected for ',fid
 - 
 
315 
      stop
 - 
 
316 
 903  print*,'Problem with level coefficients from ',trim(cstfile)
 - 
 
317 
      stop
 - 
 
318 
 904  print*,'Cannot read surface pressure from ',trim(cstfile)
 - 
 
319 
      stop
 - 
 
320 
 905  print*,'Cannot find time ',time,' on ',fid
 - 
 
321 
      stop
 - 
 
322 
 906  print*,'Unable to get grid info ',fid
 - 
 
323 
      stop
 - 
 
324
- 
 
325 
      end
 314 
      end
326 
 
 315 
 
327 
c     ------------------------------------------------------------
 316 
c     ------------------------------------------------------------
328 
c     Read wind information
 317 
c     Read wind information
329 
c     ------------------------------------------------------------
 318 
c     ------------------------------------------------------------
330 
 
 319 
 
331 
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
 320 
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
332 
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
 321 
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
333 
     >                       timecheck)
 322 
     >                       timecheck)
334 
 
 323 
 
335 
c     Read the wind component <fieldname> from the file with identifier
 324 
c     Read the wind component <fieldname> from the file with identifier
336 
c     <fid> and save it in the 3d array <field>. The vertical staggering
325 
c     <fid> and save it in the 3d array <field>. The vertical staggering
337 
c     information is provided in <stagz> and gives the reference to either
 326 
c     information is provided in <stagz> and gives the reference to either
338 
c     the layer or level field from <input_grid>. A consistency check is
 327 
c     the layer or level field from <input_grid>. A consistency check is
339 
c     performed to have an agreement with the grid specified by <xmin,xmax,
 328 
c     performed to have an agreement with the grid specified by <xmin,xmax,
340 
c     ymin,ymax,dx,dy,nx,ny,nz>.
 329 
c     ymin,ymax,dx,dy,nx,ny,nz>.
341 
 
 330 
 
- 
 
 331 
      use netcdf
- 
 
 332 
 
342 
      implicit none
 333 
      implicit none
343 
 
 334 
 
344 
c     Declaration of variables and parameters
 335 
c     Declaration of variables and parameters
345 
      integer      fid                 ! File identifier
 336 
      integer      fid                 ! File identifier
346 
      character*80 fieldname           ! Name of the wind field
 337 
      character*80 fieldname           ! Name of the wind field
347 
      integer      nx,ny,nz            ! Dimension of fields
 338 
      integer      nx,ny,nz            ! Dimension of fields
348 
      real         field(nx,ny,nz)     ! 3d wind field
 339 
      real         field(nx,ny,nz)     ! 3d wind field
349 
      real         stagz               ! Staggering in the z direction
 340 
      real         stagz               ! Staggering in the z direction
350 
      real         mdv                 ! Missing data flag
 341 
      real         mdv                 ! Missing data flag
351 
      real         xmin,xmax,ymin,ymax ! Domain size
 342 
      real         xmin,xmax,ymin,ymax ! Domain size
352 
      real         dx,dy               ! Horizontal resolution
 343 
      real         dx,dy               ! Horizontal resolution
353 
      real         time                ! Time
 344 
      real         time                ! Time
354 
      character*80 timecheck           ! Either 'yes' or 'no'
 345 
      character*80 timecheck           ! Either 'yes' or 'no'
355 
 
 346 
 
356 
c     Numerical and physical parameters
 347 
c     Numerical and physical parameters
357 
      real        eps                 ! Numerical epsilon
 348 
      real        eps                 ! Numerical epsilon
358 
      parameter  (eps=0.001)
 349 
      parameter  (eps=0.001)
359 
      real        notimecheck         ! 'Flag' for no time check
 350 
      real        notimecheck         ! 'Flag' for no time check
360 
      parameter  (notimecheck=7.26537)
 351 
      parameter  (notimecheck=7.26537)
361 
 
 352 
 
362 
c     Netcdf variables
 353 
c     Netcdf variables
363 
      integer      ierr
 354 
      integer      ierr
364 
      character*80 varname
 355 
      character*80 varname
365 
      integer      vardim(4)
 356 
      integer      vardim(4)
366 
      real         varmin(4),varmax(4)
 357 
      real         varmin(4),varmax(4)
367 
      real         stag(4)
 358 
      real         stag(4)
368 
      integer      ndim
 359 
      integer      ndim
369 
      real         times(10)
 360 
      real         times(10)
370 
      integer      ntimes
 361 
      integer      ntimes
371 
      character*80 cstfile
 362 
      character*80 cstfile
372 
      integer      cstid
 363 
      integer      cstid
373 
      real         aklay(200),bklay(200),aklev(200),bklev(200)
 364 
      real         aklay(200),bklay(200),aklev(200),bklev(200)
374 
      real         ps(nx,ny)
 365 
      real         ps(nx,ny)
- 
 
 366 
      integer      dimids (nf90_max_var_dims)
- 
 
 367 
      character*80 dimname(nf90_max_var_dims)
- 
 
 368 
      integer      varid
- 
 
 369 
      integer      cdfid
375 
 
 370 
 
376 
c     Auxiliary variables
 371 
c     Auxiliary variables
377 
      integer      isok
 372 
      integer      isok
378 
      integer      i,j,k
 373 
      integer      i,j,k
379 
      integer      nz1
 374 
      integer      nz1
380 
      real         rtime
 375 
      real         rtime
381 
 
 376 
 
382 
c     Read variable definition - for P, PLEV and PLAY: load also ak,bk
 - 
 
383 
      if ( ( fieldname.eq.'PLEV' ).or.
 - 
 
384 
     >     ( fieldname.eq.'PLAY' ).or.
 - 
 
385 
     >     ( fieldname.eq.'P'    ) )
 - 
 
386 
     >then
 - 
 
387 
         call getcfn(fid,cstfile,ierr)
 - 
 
388 
         if (ierr.ne.0) goto 905
 - 
 
389 
         call cdfopn(cstfile,cstid,ierr)
 - 
 
390 
         if (ierr.ne.0) goto 905
 - 
 
391 
         call getlevs(cstid,nz1,aklev,bklev,aklay,bklay,ierr)
 - 
 
392 
         if (ierr.ne.0) goto 905
 - 
 
393 
         call clscdf(cstid,ierr)
 - 
 
394 
         if (ierr.ne.0) goto 905
 - 
 
395 
         varname = 'PS'
 - 
 
396 
         call getdef(fid,varname,ndim,mdv,vardim,
 - 
 
397 
     >               varmin,varmax,stag,ierr)
 - 
 
398 
         vardim(3) = nz1
 - 
 
399 
         if (ierr.ne.0) goto 900
 - 
 
400 
      else
 - 
 
401 
         varname = fieldname
 - 
 
402 
         call getdef(fid,varname,ndim,mdv,vardim,
 - 
 
403 
     >               varmin,varmax,stag,ierr)
 - 
 
404 
         if (ierr.ne.0) goto 900
 - 
 
405 
         stagz=stag(3)
 - 
 
406 
      endif
 - 
 
407 
 
 - 
 
408 
c     Get time information (set time to first one in the file)
 - 
 
409 
      call gettimes(fid,times,ntimes,ierr)
 - 
 
410 
      if (ierr.ne.0) goto 902
 - 
 
411 
      isok=0
 - 
 
412 
      do i=1,ntimes
 - 
 
413 
         if (abs(time-times(i)).lt.eps) then
 - 
 
414 
            isok = 1
 - 
 
415 
            rtime = times(i)
 - 
 
416 
         elseif (timecheck.eq.'no') then
 - 
 
417 
            isok = 1
 - 
 
418 
            rtime = times(1)
 - 
 
419 
         endif
 - 
 
420 
      enddo
 - 
 
421 
      if ( isok.eq.0 ) goto 904
 - 
 
422 
 
 - 
 
423 
c     Read  field
 377 
c     Get varid
424 
      if ( ( fieldname.eq.'P' ).or.(fieldname.eq.'PLAY') )  then       ! P, PLAY
 - 
 
425 
         stagz   = -0.5
 - 
 
426 
         varname = 'PS'
 - 
 
427 
         call getdat(fid,varname,rtime,0,ps,ierr)
 378 
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
428 
         if (ierr.ne.0) goto 903
 - 
 
429 
         do i=1,nx
 - 
 
430 
            do j=1,ny
 - 
 
431 
               do k=1,nz
 - 
 
432 
                  field(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
 379 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
433 
               enddo
 - 
 
434 
            enddo
 - 
 
435 
         enddo
 - 
 
436
380
437 
      elseif ( fieldname.eq.'PLEV' )  then                             ! PLEV
 381 
c     Get number of vertical levels -> vardim(3)
438 
         stagz   = 0.
 - 
 
439 
         varname = 'PS'
 382 
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
440 
         call getdat(fid,varname,rtime,0,ps,ierr)
 383 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
441 
         if (ierr.ne.0) goto 903
 384 
      ierr = nf90_inquire_variable(fid, varid,
442 
         do i=1,nx
 - 
 
443 
            do j=1,ny
 - 
 
444 
               do k=1,nz
 385 
     >                                   dimids = dimids(1:ndim))
445 
                  field(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
 386 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
446 
               enddo
 387 
      do i=1,ndim
447 
            enddo
 388 
           ierr = nf90_inquire_dimension(fid, dimids(i),
448 
         enddo
 - 
 
449 
 
 - 
 
450 
      else                                                             ! Other fields
 389 
     >                               name = dimname(i) )
451 
         varname=fieldname
 390 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
452 
         call getdat(fid,varname,rtime,0,field,ierr)
 391 
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
453 
         if (ierr.ne.0) goto 903
 392 
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
454 
 
 - 
 
455 
      endif
 393 
      enddo
456 
 
 394 
 
- 
 
 395 
c     Read data
- 
 
 396 
      varname=fieldname
- 
 
 397 
      ierr = nf90_get_var(fid,varid,field)
- 
 
 398 
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
- 
 
 399
457 
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
 400 
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
458 
      if ( vardim(3).eq.1 ) then
 401 
      if ( vardim(3).eq.1 ) then
459 
         do i=1,nx
 402 
         do i=1,nx
460 
            do j=1,ny
 403 
            do j=1,ny
461 
               do k=1,nz
 404 
               do k=1,nz
462 
                  field(i,j,k) = field(i,j,1)
 405 
                  field(i,j,k) = field(i,j,1)
463 
               enddo
 406 
               enddo
464 
            enddo
 407 
            enddo
465 
         enddo
 408 
         enddo
466 
      endif
 409 
      endif
467
410
468 
 
 - 
 
469 
 
 - 
 
470 
c     Exception handling
 411 
c     Exit point
471 
      return
 412 
      return
472
413
473 
 900  print*,'Cannot retrieve definition for ',trim(varname),'  ',fid
 - 
 
474 
      stop
 - 
 
475 
 901  print*,'Grid inconsistency detected for ',trim(varname),'  ',fid
 - 
 
476 
      stop
 - 
 
477 
 902  print*,'Cannot retrieve time for ',trim(varname),'  ',fid
 - 
 
478 
      stop
 - 
 
479 
 903  print*,'Cannot load wind component ',trim(varname),'  ',fid
 - 
 
480 
      stop
 - 
 
481 
 904  print*,'Cannot load time ',time,' for ',trim(varname),'  ',fid
 - 
 
482 
      stop
 - 
 
483 
 905  print*,'Cannot load time vertical grid AK, BK from file  ',fid
 - 
 
484 
      stop
 - 
 
485
- 
 
486 
      end
 414 
      end
487 
 
 415 
 
488 
c     ------------------------------------------------------------
 416 
c     ------------------------------------------------------------
489 
c     Close input file
 417 
c     Close input file
490 
c     ------------------------------------------------------------
 418 
c     ------------------------------------------------------------
491 
 
 419 
 
492 
      subroutine input_close(fid)
 420 
      subroutine input_close(fid)
493 
 
 421 
 
494 
c     Close the input file with file identifier <fid>.
 422 
c     Close the input file with file identifier <fid>.
495 
 
 423 
 
- 
 
 424 
      use netcdf
- 
 
 425 
 
496 
      implicit none
 426 
      implicit none
497 
 
 427 
 
498 
c     Declaration of subroutine parameters
 428 
c     Declaration of subroutine parameters
499 
      integer fid
 429 
      integer fid
500 
 
 430 
 
501 
c     Auxiliary variables
 431 
c     Auxiliary variables
502 
      integer ierr
 432 
      integer ierr
503 
 
 433 
 
504 
c     Close file
 434 
c     Close file
505 
      call clscdf(fid,ierr)
 435 
      ierr = NF90_CLOSE(fid)
- 
 
 436 
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
506
437
507 
      end
 438 
      end