Subversion Repositories lagranto.icon

Rev

Rev 8 | 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 height <zb> and     *
16
c     * the model level height (given on the W grid). The number *
17
c     * of levels is given by <nz>. Finally, the retrieval of the* 
18
c     * wind <field> with name <fieldname>  is characterised by  *
19
c     a <time> and a missing data value <mdv>.                   *
20
c     *                                                          *
21
c     * Author: Michael Sprenger, Spring 2011/2012               *
22
c     ************************************************************
23
 
24
c     ------------------------------------------------------------
25
c     Open input file
26
c     ------------------------------------------------------------
27
 
28
      subroutine input_open (fid,filename)
29
 
30
c     Open the input file with filename <filename> and return the
31
c     file identifier <fid> for further reference. 
32
 
33
      use netcdf
34
      implicit none
35
 
36
c     Declaration of subroutine parameters
37
      integer      fid              ! File identifier
38
      character*80 filename         ! Filename
39
 
40
c     Declaration of auxiliary variables
41
      integer      ierr
42
 
43
c     Open netcdf file
44
      ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
45
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
46
 
47
      end
48
 
49
 
50
c     ------------------------------------------------------------
51
c     Read information about the grid
52
c     ------------------------------------------------------------
53
 
54
      subroutine input_grid 
55
     >             (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
56
     >              time,pollon,pollat,polgam,z3,zb,nz,stagz,timecheck)
57
 
58
c     Read grid information at <time> from file with identifier <fid>. 
59
c     The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>
60
c     with pole position at <pollon,pollat> and grid dimension <nx,ny>.
61
c     The 3d arrays <z3(nx,ny,nz)> gives the vertical coordinates, either
62
c     on the staggered or unstaggered grid (with <stagz> as the flag).
63
c     The surface height is given in <zb(nx,ny)>. If <fid> is negative, 
64
c     only the grid dimensions and grid parameters (xmin...pollat,nz) are 
65
c     determined and returned (this is needed for dynamical allocation of 
66
c     memory).
67
 
68
      use netcdf
69
      implicit none
70
 
71
c     Declaration of subroutine parameters 
72
      integer      fid                  ! File identifier
73
      real         xmin,xmax,ymin,ymax  ! Domain size
74
      real         dx,dy                ! Horizontal resolution
75
      integer      nx,ny,nz             ! Grid dimensions
76
      real         pollon,pollat,polgam ! Longitude and latitude of pole
77
      real         z3(nx,ny,nz)         ! Staggered levels
78
      real         zb(nx,ny)            ! Surface pressure
79
      real         time                 ! Time of the grid information
80
      real         stagz                ! Vertical staggering (0 or -0.5)
81
      character*80 fieldname            ! Variable from which to take grid info
82
      character*80 timecheck            ! Either 'yes' or 'no'
83
 
84
c     Numerical and physical parameters
85
      real          eps                ! Numerical epsilon
86
      parameter    (eps=0.001)
87
 
88
c     Name of the constants file
89
      character*80  constfile
90
      parameter     (constfile = 'ICONCONST')
91
 
92
c     netCDF variables
93
      real         lon(5000)
94
      real         lat(5000)
95
 
96
c     Auxiliary variables
97
      integer      ierr   
98
      integer      varid
99
      integer      dimid
100
      integer      dimids(100)
101
      character*80 dimname(100)   
102
      integer      ndim
103
      integer      i,j,k
104
      real         tmp3(nx,ny,nz)
8 michaesp 105
      real         stag3(nx,ny,nz+1)
3 michaesp 106
      integer      cdfid 
107
 
108
c     fid can be negative, take its absolute value
109
      cdfid = abs(fid)
110
 
111
c     Open ICONCONST      
112
      ierr = NF90_OPEN(constfile,nf90_nowrite, cdfid)
113
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
114
 
115
c     Get <ndim>, <dimids> and <dimname>
116
      ierr = NF90_INQ_VARID(cdfid,'z_mc',varid)
117
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
118
      ierr = nf90_inquire_variable(cdfid, varid,ndims  = ndim)
119
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
120
      ierr = nf90_inquire_variable(cdfid, varid,dimids = dimids(1:ndim))
121
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
122
      do i=1,ndim
123
         ierr = nf90_inquire_dimension(cdfid, dimids(i), 
124
     >                               name = dimname(i) )
125
         IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
126
      enddo
127
 
128
c     Check that dimensions are OK
129
      if ( ndim.ne.3 ) then
130
        print*,' ERROR: z_mc must be 3D ',ndim
131
        stop
132
      endif
133
      if ( dimname(2).ne.'lat' ) then
134
        print*,' ERROR: dimname(2) must be lat ',trim(dimname(2))
135
        stop
136
      endif
137
      if ( dimname(1).ne.'lon' ) then
138
        print*,' ERROR: dimname(1) must be lon ',trim(dimname(1))
139
        stop
140
      endif
141
 
142
c     Get lon coordinates        
143
      ierr = nf90_inq_dimid(cdfid,'lon', dimid)
144
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
145
      ierr = nf90_inquire_dimension(cdfid, dimid, len = nx)
146
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
147
      ierr = NF90_INQ_VARID(cdfid,'lon',varid)
148
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
149
      ierr = NF90_GET_VAR(cdfid,varid,lon(1:nx))
150
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 
151
 
152
c     Get lat coordinates        
153
      ierr = nf90_inq_dimid(cdfid,'lat', dimid)
154
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
155
      ierr = nf90_inquire_dimension(cdfid, dimid, len = ny)
156
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
157
      ierr = NF90_INQ_VARID(cdfid,'lat',varid)
158
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
159
      ierr = NF90_GET_VAR(cdfid,varid,lat(1:ny))
160
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 
161
 
162
c     Get number of vertical levels
163
      ierr = nf90_inquire_dimension(cdfid, dimids(3), len = nz)
164
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
165
 
8 michaesp 166
c     Handling of vertical staggering 
167
      nz = nz -1
168
 
3 michaesp 169
c     Set parameters
170
      pollon = 0.
171
      pollat = 90.
172
      polgam = 0.
173
      xmin   = lon(1)
174
      ymin   = lat(1)
175
      xmax   = lon(nx)
176
      ymax   = lat(ny)
177
      dx     = ( xmax - xmin ) / real(nx-1) 
178
      dy     = ( ymax - ymin ) / real(ny-1)
179
      if ( xmin.lt.-180.) then
180
        xmin = xmin + 360.
181
        xmax = xmax + 360.
182
      else if ( xmax.gt.360. ) then
183
        xmin = xmin - 360.
184
        xmax = xmax - 360.
185
      endif
186
      stagz = 0.
187
 
188
c     If <fid<0>, nothing further to do - we have the grid dimensions 
189
      if ( fid.lt.0 ) goto 100
190
 
191
c     Read topography 
192
      ierr = NF90_INQ_VARID(cdfid,'TOPOGRAPHY',varid)
193
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
194
      ierr = NF90_GET_VAR(cdfid,varid,zb)
195
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 
196
 
197
c     Read 3D height field
198
      ierr = NF90_INQ_VARID(cdfid,'z_mc',varid)
199
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
8 michaesp 200
      ierr = NF90_GET_VAR(cdfid,varid,stag3)
3 michaesp 201
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 
202
 
8 michaesp 203
c     Handling of vertical staggering - destagger
204
      do i=1,nx
205
         do j=1,ny
206
         	do k=1,nz
207
         	   tmp3(i,j,k) = 0.5 * ( stag3(i,j,k) + stag3(i,j,k+1) )
208
         	enddo
209
         enddo
210
      enddo   	
211
 
3 michaesp 212
c     Check whether vertical axis is descending - vertical flip
213
      if ( tmp3(1,1,1).lt.tmp3(1,1,nz) ) then
214
        print*,' ERROR: vertical axis must be descending... Stop'
215
        stop
216
      endif
217
      do i=1,nx
218
        do j=1,ny
219
           do k=1,nz
220
              z3(i,j,k) = tmp3(i,j,nz-k+1)
221
           enddo
222
        enddo
223
      enddo
224
 
225
c     Exit point
226
 100  continue
227
 
228
c     Close constants file
229
      ierr = NF90_CLOSE(cdfid)
230
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)   
231
 
232
      end
233
 
234
c     ------------------------------------------------------------
235
c     Read wind information
236
c     ------------------------------------------------------------
237
 
238
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
239
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
240
     >                       timecheck)
241
 
242
c     Read the wind component <fieldname> from the file with identifier
243
c     <fid> and save it in the 3d array <field>. The vertical staggering 
244
c     information is provided in <stagz> and gives the reference to either
245
c     the layer or level field from <input_grid>. A consistency check is
246
c     performed to have an agreement with the grid specified by <xmin,xmax,
247
c     ymin,ymax,dx,dy,nx,ny,nz>.
248
 
249
      use netcdf
250
      implicit none
251
 
252
c     Declaration of variables and parameters
253
      integer      fid                 ! File identifier
254
      character*80 fieldname           ! Name of the wind field
255
      integer      nx,ny,nz            ! Dimension of fields
256
      real         field(nx,ny,nz)     ! 3d wind field
257
      real         stagz               ! Staggering in the z direction
258
      real         mdv                 ! Missing data flag
259
      real         xmin,xmax,ymin,ymax ! Domain size
260
      real         dx,dy               ! Horizontal resolution
261
      real         time                ! Time
262
      character*80 timecheck           ! Either 'yes' or 'no'
263
 
264
c     Auxiliary variables
265
      integer      i,j,k
266
      integer      nlon,nlat,nlev
267
      real         tmp2(nx,ny)
268
      real         tmp3(nx,ny,nz)
8 michaesp 269
      real         stag3(nx,ny,nz+1)
3 michaesp 270
      integer      ierr
271
      integer      varid
272
      integer      dimid
273
      integer      dimids(100)
274
      character*80 dimname(100)   
275
      integer      ndim
9 michaesp 276
 
277
c     ----- Get dimension and check whether OK --------------------------------- 
3 michaesp 278
 
279
c     Get <ndim>, <dimids> and <dimname>
280
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
281
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
282
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
283
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
284
      ierr = nf90_inquire_variable(fid, varid,dimids = dimids(1:ndim))
285
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
286
      do i=1,ndim
287
         ierr = nf90_inquire_dimension(fid, dimids(i), 
288
     >                               name = dimname(i) )
289
         IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
290
      enddo
9 michaesp 291
 
292
      if ( (ndim.ne.3).and.(ndim.ne.4) ) then
293
        print*,' Data must be either 3D or 4D'
294
        stop
295
      endif
3 michaesp 296
 
9 michaesp 297
c     ----- Read 4D data ------------------------------------------------------ 
298
 
299
      if ( ndim.ne.4 ) goto 100
300
 
3 michaesp 301
c     Check that dimensions are OK
302
      if ( dimname(4).ne.'time' ) then
303
        print*,' ERROR: dimname(4) must be time '
304
        print*,trim(fieldname),' / ',trim(dimname(4))
305
        stop
306
      endif
307
      if ( dimname(2).ne.'lat' ) then
308
        print*,' ERROR: dimname(2) must be lat ' 
309
        print*,trim(fieldname),' / ',trim(dimname(2))
310
        stop
311
      endif    
312
      if ( dimname(1).ne.'lon' ) then
313
        print*,' ERROR: dimname(1) must be lon ' 
314
        print*,trim(fieldname),' / ',trim(dimname(2))
315
        stop
316
      endif   
317
 
318
c     Check grid dimensions       
319
      ierr = nf90_inq_dimid(fid,dimname(1), dimid)
320
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
321
      ierr = nf90_inquire_dimension(fid, dimid, len = nlon)
322
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
323
      ierr = nf90_inq_dimid(fid,dimname(2), dimid)
324
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
325
      ierr = nf90_inquire_dimension(fid, dimid, len = nlat)
326
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
327
      ierr = nf90_inq_dimid(fid,dimname(3), dimid)
328
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
329
      ierr = nf90_inquire_dimension(fid, dimid, len = nlev)
330
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
331
      if ( (nlon.ne.nx).or.(nlat.ne.ny).or.
8 michaesp 332
     >  (nlev.ne.nz).and.(nlev.ne.1).and.(nlev.ne.(nz+1)) )
3 michaesp 333
     >then
8 michaesp 334
         print*,' ERROR: grid mismatch between ICONCONST and P file'
3 michaesp 335
         stop
336
      endif
337
 
8 michaesp 338
c     Get varid
3 michaesp 339
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
340
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
8 michaesp 341
 
342
c     Read 2D variable
343
      if ( nlev.eq.1 ) then
3 michaesp 344
        ierr = NF90_GET_VAR(fid,varid,tmp2)
345
        IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 
346
        do i=1,nx
347
            do j=1,ny
348
               do k=1,nz
349
                  tmp3(i,j,k) = tmp2(i,j)
350
               enddo
351
            enddo
352
         enddo
8 michaesp 353
      endif
354
 
355
c     Read 3D variable, destagger if necessary      
356
      if ( nlev.eq.nz ) then
3 michaesp 357
        ierr = NF90_GET_VAR(fid,varid,tmp3)
358
        IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 
8 michaesp 359
      elseif ( nlev.eq.(nz+1) ) then  
360
        ierr = NF90_GET_VAR(fid,varid,stag3)
361
        IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 
362
        do i=1,nx
363
          do j=1,ny
364
             do k=1,nz
365
                tmp3(i,j,k) = 0.5 * ( stag3(i,j,k) + stag3(i,j,k+1) )
366
             enddo
367
          enddo
368
        enddo  
3 michaesp 369
      endif
370
 
371
c     Flip vertically 
372
      do i=1,nx
373
        do j=1,ny
374
           do k=1,nz
375
              field(i,j,k) = tmp3(i,j,nz-k+1)
376
           enddo
377
        enddo
378
      enddo
379
 
9 michaesp 380
c     Exit point for reading 4D data
381
 100  continue
382
 
383
c     ----- Read 3D data ------------------------------------------------------ 
384
 
385
      if ( ndim.ne.3 ) goto 200
386
 
387
c     Check that dimensions are OK
388
      if ( dimname(3).ne.'time' ) then
389
        print*,' ERROR: dimname(3) must be time '
390
        print*,trim(fieldname),' / ',trim(dimname(3))
391
        stop
392
      endif
393
      if ( dimname(2).ne.'lat' ) then
394
        print*,' ERROR: dimname(2) must be lat ' 
395
        print*,trim(fieldname),' / ',trim(dimname(2))
396
        stop
397
      endif    
398
      if ( dimname(1).ne.'lon' ) then
399
        print*,' ERROR: dimname(1) must be lon ' 
400
        print*,trim(fieldname),' / ',trim(dimname(2))
401
        stop
402
      endif   
403
 
404
c     Check grid dimensions       
405
      ierr = nf90_inq_dimid(fid,dimname(1), dimid)
406
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
407
      ierr = nf90_inquire_dimension(fid, dimid, len = nlon)
408
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
409
      ierr = nf90_inq_dimid(fid,dimname(2), dimid)
410
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
411
      ierr = nf90_inquire_dimension(fid, dimid, len = nlat)
412
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
413
      if ( (nlon.ne.nx).or.(nlat.ne.ny) ) then
414
         print*,' ERROR: grid mismatch between ICONCONST and P file'
415
         stop
416
      endif
417
 
418
c     Get varid
419
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
420
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
421
 
422
c     Read 2D variable
423
      ierr = NF90_GET_VAR(fid,varid,tmp2)
424
      IF ( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 
425
      do i=1,nx
426
          do j=1,ny
427
             do k=1,nz
428
                field(i,j,k) = tmp2(i,j)
429
             enddo
430
          enddo
431
      enddo
432
 
433
c     Exit point for reading 4D data
434
 200  continue
3 michaesp 435
 
436
      end
437
 
438
c     ------------------------------------------------------------
439
c     Close input file
440
c     ------------------------------------------------------------
441
 
442
      subroutine input_close(fid)
443
 
444
c     Close the input file with file identifier <fid>.
445
 
446
      use netcdf
447
      implicit none
448
 
449
c     Declaration of subroutine parameters
450
      integer fid
451
 
452
c     Auxiliary variables
453
      integer ierr
454
 
455
c     Close file
456
      ierr = NF90_CLOSE(fid)
457
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
458
 
459
      end