Subversion Repositories lagranto.um

Rev

Rev 3 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
c     ************************************************************
2
c     * This package provides input routines to read the wind    *
3
c     * and other fields from IVE necdf files. The routines are  *
4
c     *                                                          *
5
c     * 1) input_open  : to open a data file                     *
6
c     * 2) input_grid  : to read the grid information, including *
7
c     *                  the vertical levels                     *
8
c     * 3) input_wind  : to read the wind components             *
9
c     * 4) input_close : to close an input file                  *
10
c     *                                                          *
11
c     * The file is characterised by an filename <filename> and  *
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  *
14
c     * rotated grid is given by <pollon,pollat>. The vertical   *
15
c     * grid is characterised by the surface pressure <ps> and   *
16
c     * the pressure at staggered <slev> and unstaggered <ulev>  *
17
c     * levels. The number of levels is given by <nz>. Finally,  *
18
c     * the retrieval of the wind <field> with name <fieldname>  *
19
c     * is characterised by a <time> and a missing data value    *
20
c     * <mdv>.                                                   *
21
c     *                                                          *
22
c     * Author: Michael Sprenger, Autumn 2008                    *
23
c     ************************************************************
24
 
25
c     ------------------------------------------------------------
26
c     Open input file
27
c     ------------------------------------------------------------
28
 
29
      subroutine input_open (fid,filename)
30
 
31
c     Open the input file with filename <filename> and return the
32
c     file identifier <fid> for further reference. 
33
 
7 michaesp 34
      use netcdf
35
 
3 michaesp 36
      implicit none
37
 
38
c     Declaration of subroutine parameters
39
      integer      fid              ! File identifier
40
      character*80 filename         ! Filename
41
 
42
c     Declaration of auxiliary variables
43
      integer      ierr
44
 
7 michaesp 45
c     Open netcdf file
46
      ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
47
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
48
 
3 michaesp 49
c     Exception handling
50
      return
51
 
52
 900  print*,'Cannot open input file ',trim(filename)
53
      stop
54
 
55
      end
56
 
57
 
58
c     ------------------------------------------------------------
59
c     Read information about the grid
60
c     ------------------------------------------------------------
61
 
62
      subroutine input_grid 
63
     >                   (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
64
     >                    time,pollon,pollat,p3,ps,nz,ak,bk,stagz,
65
     >                    timecheck)
66
 
67
c     Read grid information at <time> from file with identifier <fid>. 
68
c     The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>
69
c     with pole position at <pollon,pollat> and grid dimension <nx,ny>.
70
c     The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either
71
c     on the staggered or unstaggered grid (with <stagz> as the flag).
72
c     The surface pressure is given in <ps(nx,ny)>. If <fid> is negative, 
73
c     only the grid dimensions and grid parameters (xmin...pollat,nz) are 
74
c     determined and returned (this is needed for dynamical allocation of 
75
c     memory).
76
 
7 michaesp 77
      use netcdf
78
 
3 michaesp 79
      implicit none
80
 
81
c     Declaration of subroutine parameters 
82
      integer      fid                 ! File identifier
83
      real         xmin,xmax,ymin,ymax ! Domain size
84
      real         dx,dy               ! Horizontal resolution
85
      integer      nx,ny,nz            ! Grid dimensions
86
      real         pollon,pollat       ! Longitude and latitude of pole
87
      real         p3(nx,ny,nz)        ! Staggered levels
88
      real         ps(nx,ny)           ! Surface pressure
89
      real         time                ! Time of the grid information
90
      real         ak(nz),bk(nz)       ! Ak and Bk for layers or levels
91
      real         stagz               ! Vertical staggering (0 or -0.5)
92
      character*80 fieldname           ! Variable from which to take grid info
93
      character*80 timecheck           ! Either 'yes' or 'no'
94
 
95
c     Numerical and physical parameters
96
      real          eps                 ! Numerical epsilon
97
      parameter    (eps=0.001)
98
 
99
c     Netcdf variables
100
      integer      vardim(4)
101
      real         varmin(4),varmax(4)
102
      real         mdv
103
      real         stag(4)
104
      integer      ndim
105
      character*80 cstfile
106
      integer      cstid
7 michaesp 107
      real         times
3 michaesp 108
      integer      ntimes
109
      real         aklay(nz),bklay(nz),aklev(nz),bklev(nz)
110
      integer      nvars
111
      character*80 vars(100)
7 michaesp 112
      integer        dimids (nf90_max_var_dims)
113
      character*80   dimname(nf90_max_var_dims)
3 michaesp 114
 
115
c     Auxiliary varaibles
116
      integer      ierr       
117
      integer      i,j,k
118
      integer      isok
119
      real         tmp(200)
120
      character*80 varname
121
      real         rtime
7 michaesp 122
      integer      varid
123
      integer      cdfid
3 michaesp 124
 
7 michaesp 125
c     Set file identifier
3 michaesp 126
      if (fid.lt.0) then
7 michaesp 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)
3 michaesp 144
 
7 michaesp 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)
3 michaesp 177
 
7 michaesp 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)
3 michaesp 188
 
7 michaesp 189
c     We want the longitudes within -180 ... + 180
190
      if ( xmin.lt.-180.) then
3 michaesp 191
            xmin = xmin + 360.
192
            xmax = xmax + 360.
7 michaesp 193
      else if ( xmax.gt.360. ) then
3 michaesp 194
            xmin = xmin - 360.
195
            xmax = xmax - 360.
7 michaesp 196
      endif
3 michaesp 197
 
7 michaesp 198
c     Get name of constants file -> cstfile
199
      ierr  = nf90_get_att(cdfid,nf90_global,
200
     >                    "constants_file_name ", cstfile )
201
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
3 michaesp 202
 
7 michaesp 203
c     Get pole position from constants file
204
      ierr = NF90_OPEN(TRIM(cstfile),nf90_nowrite, cstid)
205
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
206
      varname='pollon'
207
      ierr = NF90_INQ_VARID(cstid,varname,varid)
208
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
209
      ierr = nf90_get_var(cstid,varid,pollon)
210
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
211
      varname='pollat'
212
      ierr = NF90_INQ_VARID(cstid,varname,varid)
213
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
214
      ierr = nf90_get_var(cstid,varid,pollat)
215
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
216
      ierr = NF90_CLOSE(cstid)
217
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
218
 
219
c     Skip the rest if fid < 0: 
220
      if ( fid.lt.0 ) goto 900      
221
 
222
c     Get ak and bk from constants file
223
      ierr = NF90_OPEN(TRIM(cstfile),nf90_nowrite, cstid)
224
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
225
      varname='aklev'
226
      ierr = NF90_INQ_VARID(cstid,varname,varid)
227
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
228
      ierr = nf90_get_var(cstid,varid,aklev)
229
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
230
      varname='bklev'
231
      ierr = NF90_INQ_VARID(cstid,varname,varid)
232
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
233
      ierr = nf90_get_var(cstid,varid,bklev)
234
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
235
      varname='aklay'
236
      ierr = NF90_INQ_VARID(cstid,varname,varid)
237
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
238
      ierr = nf90_get_var(cstid,varid,aklay)
239
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
240
      varname='bklay'
241
      ierr = NF90_INQ_VARID(cstid,varname,varid)
242
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
243
      ierr = nf90_get_var(cstid,varid,bklay)
244
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
245
      ierr = NF90_CLOSE(cstid)
246
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
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)
254
c     print*,'BKLAY : '
255
c     print*,(bklay(i),i=1,20)
256
 
257
c     Get time information (check if time is correct)
258
      varname = 'time'
259
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
260
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
261
      ierr = nf90_get_var(cdfid,varid,times)
262
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
263
      isok=0
264
      if (abs(time-times).lt.eps) then
3 michaesp 265
               isok = 1
7 michaesp 266
               rtime = times
267
      elseif (timecheck.eq.'no') then
3 michaesp 268
               isok = 1
7 michaesp 269
               rtime = times
270
      endif
271
      if ( isok.eq.0 ) then
272
         print*,' ERROR: time ',rtime,' not found on netCDF file' 
273
         stop
274
      endif
275
 
276
c     print*,'TIME : ',rtime
3 michaesp 277
 
7 michaesp 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
 
285
c     Calculate layer and level pressures
286
      do i=1,nx
287
         do j=1,ny
3 michaesp 288
               do k=1,nz
289
                  if ( abs(stagz).lt.eps ) then
290
                     p3(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
291
                  else
292
                     p3(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
293
                  endif
294
               enddo
295
         enddo
7 michaesp 296
      enddo
3 michaesp 297
 
7 michaesp 298
c     Set the ak and bk for the vertical grid
299
      do k=1,nz
3 michaesp 300
            if ( abs(stagz).lt.eps ) then
301
               ak(k)=aklev(k)
302
               bk(k)=bklev(k)
303
            else
304
               ak(k)=aklay(k)
305
               bk(k)=bklay(k)
306
            endif
7 michaesp 307
      enddo
3 michaesp 308
 
7 michaesp 309
c     Exit point
310
 900  continue    
311
 
3 michaesp 312
      return
313
 
314
      end
315
 
316
c     ------------------------------------------------------------
317
c     Read wind information
318
c     ------------------------------------------------------------
319
 
320
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
321
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
322
     >                       timecheck)
323
 
324
c     Read the wind component <fieldname> from the file with identifier
325
c     <fid> and save it in the 3d array <field>. The vertical staggering 
326
c     information is provided in <stagz> and gives the reference to either
327
c     the layer or level field from <input_grid>. A consistency check is
328
c     performed to have an agreement with the grid specified by <xmin,xmax,
329
c     ymin,ymax,dx,dy,nx,ny,nz>.
330
 
7 michaesp 331
      use netcdf
332
 
3 michaesp 333
      implicit none
334
 
335
c     Declaration of variables and parameters
336
      integer      fid                 ! File identifier
337
      character*80 fieldname           ! Name of the wind field
338
      integer      nx,ny,nz            ! Dimension of fields
339
      real         field(nx,ny,nz)     ! 3d wind field
340
      real         stagz               ! Staggering in the z direction
341
      real         mdv                 ! Missing data flag
342
      real         xmin,xmax,ymin,ymax ! Domain size
343
      real         dx,dy               ! Horizontal resolution
344
      real         time                ! Time
345
      character*80 timecheck           ! Either 'yes' or 'no'
346
 
347
c     Numerical and physical parameters
348
      real        eps                 ! Numerical epsilon
349
      parameter  (eps=0.001)
350
      real        notimecheck         ! 'Flag' for no time check
351
      parameter  (notimecheck=7.26537)
352
 
353
c     Netcdf variables
354
      integer      ierr
355
      character*80 varname
356
      integer      vardim(4)
357
      real         varmin(4),varmax(4)
358
      real         stag(4)
359
      integer      ndim
360
      real         times(10)
361
      integer      ntimes
362
      character*80 cstfile
363
      integer      cstid
364
      real         aklay(200),bklay(200),aklev(200),bklev(200)
365
      real         ps(nx,ny)
7 michaesp 366
      integer      dimids (nf90_max_var_dims)
367
      character*80 dimname(nf90_max_var_dims)
368
      integer      varid
369
      integer      cdfid
3 michaesp 370
 
371
c     Auxiliary variables
372
      integer      isok
373
      integer      i,j,k
374
      integer      nz1
375
      real         rtime
376
 
7 michaesp 377
c     Get varid
378
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
379
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
380
 
381
c     Get number of vertical levels -> vardim(3)    
382
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
383
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
384
      ierr = nf90_inquire_variable(fid, varid, 
385
     >                                   dimids = dimids(1:ndim))
386
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
387
      do i=1,ndim
388
           ierr = nf90_inquire_dimension(fid, dimids(i), 
389
     >                               name = dimname(i) )
390
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
391
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
392
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
3 michaesp 393
      enddo
394
 
7 michaesp 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
 
3 michaesp 400
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
401
      if ( vardim(3).eq.1 ) then
402
         do i=1,nx
403
            do j=1,ny
404
               do k=1,nz
405
                  field(i,j,k) = field(i,j,1)
406
               enddo
407
            enddo
408
         enddo
409
      endif
410
 
7 michaesp 411
c     Exit point
3 michaesp 412
      return
7 michaesp 413
 
3 michaesp 414
      end
415
 
416
c     ------------------------------------------------------------
417
c     Close input file
418
c     ------------------------------------------------------------
419
 
420
      subroutine input_close(fid)
421
 
422
c     Close the input file with file identifier <fid>.
423
 
7 michaesp 424
      use netcdf
425
 
3 michaesp 426
      implicit none
427
 
428
c     Declaration of subroutine parameters
429
      integer fid
430
 
431
c     Auxiliary variables
432
      integer ierr
433
 
434
c     Close file
7 michaesp 435
      ierr = NF90_CLOSE(fid)
436
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
437
 
3 michaesp 438
      end