Subversion Repositories lagranto.ecmwf

Rev

Rev 25 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
25 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
 
34
      use netcdf
35
 
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
 
45
c     Open netcdf file
46
      ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
47
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
48
 
49
c     Exception handling
50
      return
51
 
52
      end
53
 
54
c     ------------------------------------------------------------
55
c     Read information about the grid
56
c     ------------------------------------------------------------
57
 
58
      subroutine input_grid 
59
     >                   (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
60
     >                    time,pollon,pollat,p3,ps,nz,ak,bk,stagz,
61
     >                    timecheck)
62
 
63
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>
65
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
67
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, 
69
c     only the grid dimensions and grid parameters (xmin...pollat,nz) are 
70
c     determined and returned (this is needed for dynamical allocation of 
71
c     memory).
72
 
73
      use netcdf
74
 
75
      implicit none
76
 
77
c     Declaration of subroutine parameters 
78
      integer      fid                 ! File identifier
79
      real         xmin,xmax,ymin,ymax ! Domain size
80
      real         dx,dy               ! Horizontal resolution
81
      integer      nx,ny,nz            ! Grid dimensions
82
      real         pollon,pollat       ! Longitude and latitude of pole
83
      real         p3(nx,ny,nz)        ! Staggered levels
84
      real         ps(nx,ny)           ! Surface pressure
85
      real         time                ! Time of the grid information
86
      real         ak(nz),bk(nz)       ! Ak and Bk for layers or levels
87
      real         stagz               ! Vertical staggering (0 or -0.5)
88
      character*80 fieldname           ! Variable from which to take grid info
89
      character*80 timecheck           ! Either 'yes' or 'no'
90
 
91
c     Numerical and physical parameters
92
      real          eps                 ! Numerical epsilon
93
      parameter    (eps=0.001)
94
 
95
c     Netcdf variables
96
      integer      vardim(4)
97
      real         varmin(4),varmax(4)
98
      real         mdv
99
      real         stag(4)
100
      integer      ndim
101
      character*80 cstfile
102
      integer      cstid
103
      integer      nvars
104
      character*80 vars(100)
27 michaesp 105
      integer        dimids (nf90_max_var_dims),dimid
25 michaesp 106
      character*80   dimname(nf90_max_var_dims)
27 michaesp 107
      character*80   stdname
25 michaesp 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
27 michaesp 114
      character*80  leveltype
115
      integer       nakbktmp
116
      integer       vertical_swap
25 michaesp 117
 
118
c     Auxiliary variables
119
      integer      ierr       
120
      integer      i,j,k
121
      integer      isok
122
      real         tmp(200)
123
      character*80 varname
124
      real         rtime
125
      integer      varid
126
      integer      cdfid
127
      integer      stat
128
      real         delta
129
      integer      closear
130
      real         maxps,minps
131
 
27 michaesp 132
c     ---- Read data from netCDF file as they are ---------------------
133
 
134
c     Set file identifier
25 michaesp 135
      if (fid.lt.0) then
136
        cdfid = -fid
137
      else 
138
        cdfid = fid
139
      endif
140
 
141
c     Special handling if 3D pressure is
142
      if ( fieldname.eq.'P' ) then
143
         fieldname = 'U'
144
      endif
145
 
146
c     Get number of dimensions of variable -> ndim
147
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
148
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
149
      ierr = nf90_inquire_variable(cdfid, varid, ndims  = ndim)
150
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
151
      if ( ndim.ne.4 ) then
152
          print*,' ERROR: netCDF variables need to be 4D'
153
          print*,'      ',trim(fieldname)
154
          stop
155
      endif
156
 
157
c     Get dimensions -> vardim(1:ndim),dimname(1:ndim)
158
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
159
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
160
      ierr = nf90_inquire_variable(cdfid, varid, 
161
     >                                   dimids = dimids(1:ndim))
162
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
163
      do i=1,ndim
164
           ierr = nf90_inquire_dimension(cdfid, dimids(i), 
165
     >                               name = dimname(i) )
166
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
167
           ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))
168
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
169
      enddo
170
 
27 michaesp 171
c     Get dimension of AK,BK
172
      varname = 'nhym'
173
      ierr = NF90_INQ_DIMID(cdfid,varname,dimid)
174
      ierr = nf90_inquire_dimension(cdfid, dimid,len=nakbktmp)
175
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
176
 
25 michaesp 177
c     Check whether the list of dimensions is OK
178
      if ( ( dimname(1).ne.'lon'  ).or.
179
     >     ( dimname(2).ne.'lat'  ).or. 
180
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
181
     >     ( dimname(4).ne.'time' ) )
182
     >then
183
        print*,' ERROR: the dimensions of the variable are not correct'
184
        print*,'        expected -> lon / lat / lev / time'
185
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
186
        stop
187
      endif
188
 
189
c     Allocate memory for reading arrays
190
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
191
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
192
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
193
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
194
      allocate(lon(vardim(1)),stat=stat)
195
      if (stat.ne.0) print*,'*** error allocating array lon     ***' 
196
      allocate(lat(vardim(2)),stat=stat)
197
      if (stat.ne.0) print*,'*** error allocating array lat     ***' 
27 michaesp 198
      allocate(lev(vardim(3)),stat=stat)
199
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
25 michaesp 200
      allocate(times(vardim(4)),stat=stat)
201
      if (stat.ne.0) print*,'*** error allocating array times   ***'
27 michaesp 202
      allocate(aktmp(nakbktmp),stat=stat)
25 michaesp 203
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
27 michaesp 204
      allocate(bktmp(nakbktmp),stat=stat)
25 michaesp 205
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
206
 
27 michaesp 207
c     Get domain longitudes, latitudes and levels
25 michaesp 208
      varname = dimname(1)
209
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
210
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
211
      ierr = nf90_get_var(cdfid,varid,lon)
212
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
213
      varname = dimname(2)
214
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
215
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
216
      ierr = nf90_get_var(cdfid,varid,lat)
217
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
27 michaesp 218
      varname = dimname(3)
219
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
220
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
221
      ierr = nf90_get_var(cdfid,varid,lev)
222
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
25 michaesp 223
 
224
c     Get ak and bk
225
      varname='hyam'
226
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
227
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
228
      ierr = nf90_get_var(cdfid,varid,aktmp)
229
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
230
      varname='hybm'
231
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
232
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
233
      ierr = nf90_get_var(cdfid,varid,bktmp)
234
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
235
 
236
c     Check that unit of ak is in hPa - if necessary correct it
237
      varname='hyam'
238
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
239
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
240
      ierr = nf90_get_att(cdfid, varid, "units", units)
241
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
242
      if ( units.eq.'Pa' ) then
27 michaesp 243
         do k=1,nakbktmp
25 michaesp 244
            aktmp(k) = 0.01 * aktmp(k)
245
         enddo
246
      endif
247
 
27 michaesp 248
c     Decide whether to swap vertical levels - highest pressure at index 1
249
      vertical_swap = 1
250
      if ( (aktmp(1) + bktmp(1) * 1000.).gt.
251
     >       (aktmp(2) + bktmp(2) * 1000.) )
252
     >then
253
          vertical_swap = 0
254
      endif
255
 
25 michaesp 256
c     Get time information (check if time is correct)
257
      varname = 'time'
258
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
259
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
260
      ierr = nf90_get_var(cdfid,varid,times)
261
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
262
      isok=0
263
      do i=1,vardim(4)
264
        if (abs(time-times(i)).lt.eps) then
265
               isok = 1
266
               rtime = times(i)
267
        elseif (timecheck.eq.'no') then
268
               isok = 1
269
               rtime = times(1)
270
        endif
271
      enddo
272
      if ( isok.eq.0 ) then
273
         print*,' ERROR: time ',rtime,' not found on netCDF file' 
274
         stop
275
      endif
276
 
277
c     Read surface pressure
278
      varname='PS'
279
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
280
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
281
      ierr = nf90_get_var(cdfid,varid,tmp2)
282
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
283
 
284
c     Check that surface pressure is in hPa
27 michaesp 285
      maxps = -1.e19
286
      minps =  1.e19
25 michaesp 287
      do i=1,vardim(1)
288
        do j=1,vardim(2)
289
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
290
             if (tmp2(i,j).lt.minps) minps = tmp2(i,j)
291
        enddo
292
      enddo
293
      if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then
294
         print*,' ERROR: surface pressre PS must be in hPa'
295
         print*,'       ',maxps,minps
296
         stop
297
      endif
298
 
27 michaesp 299
c     ---- Define output of subroutine --------------------------------
300
 
301
c     If not full list of vertical levels, reduce AK,BK arrays
302
      if ( (leveltype.eq.'hybrid_sigma_pressure').and.
303
     >     (nakbktmp.ne.vardim(3) ) )
304
     >then
305
         print*,' WARNING: only subset of vertical levels used...'
306
         do k=1,vardim(3)
307
            if ( vertical_swap.eq.1 ) then
308
               aktmp(k) = aktmp( k+nakbktmp-vardim(3) )
309
               bktmp(k) = bktmp( k+nakbktmp-vardim(3) )
310
            endif
25 michaesp 311
         enddo
27 michaesp 312
      endif
25 michaesp 313
 
314
c     Set the grid dimensions and constants
315
      nx      = vardim(1)
316
      ny      = vardim(2)
317
      nz      = vardim(3)
318
      xmin    = lon(1)
319
      ymin    = lat(1)
320
      xmax    = lon(nx)
321
      ymax    = lat(ny)
322
      dx      = (xmax-xmin)/real(nx-1)
323
      dy      = (ymax-ymin)/real(ny-1)
324
      pollon  = 0.
325
      pollat  = 90.
326
      stagz   = 0.
327
      delta   = xmax-xmin-360.
328
      if (abs(delta+dx).lt.eps) then
329
          xmax    = xmax + dx
330
          nx      = nx + 1
331
          closear = 1
332
      else
333
          closear = 0
334
      endif
335
 
336
c     Save the output arrays (if fid>0) - close arrays on request
337
      if ( fid.gt.0 ) then
338
 
27 michaesp 339
c        Calculate layer pressures
340
         do i=1,vardim(1)
341
              do j=1,vardim(2)
342
                 do k=1,vardim(3)
343
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
344
                 enddo
345
              enddo
346
         enddo
347
 
348
c        Get PS - close array on demand
25 michaesp 349
         do j=1,vardim(2)
350
           do i=1,vardim(1)
351
             ps(i,j) = tmp2(i,j)
352
           enddo
353
           if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)
354
         enddo
355
 
27 michaesp 356
c        Get P3 - close array on demand + vertical swap
25 michaesp 357
         do j=1,vardim(2)
358
           do k=1,vardim(3)
359
             do i=1,vardim(1)
27 michaesp 360
               if ( vertical_swap.eq.1 ) then
361
                  p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)
362
               else
363
                  p3(i,j,k) = tmp3(i,j,k)
364
               endif
25 michaesp 365
             enddo
366
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
367
           enddo
368
         enddo
369
 
27 michaesp 370
c        Get AK,BK - vertical swap on demand
25 michaesp 371
         do k=1,vardim(3)
27 michaesp 372
              if ( vertical_swap.eq.1 ) then
373
                 ak(k) = aktmp(vardim(3)-k+1)
374
                 bk(k) = bktmp(vardim(3)-k+1)
375
              endif
25 michaesp 376
         enddo
377
 
378
      endif
379
 
27 michaesp 380
 
25 michaesp 381
      return
382
 
383
      end
384
 
385
c     ------------------------------------------------------------
386
c     Read wind information
387
c     ------------------------------------------------------------
388
 
389
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
390
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
391
     >                       timecheck)
392
 
393
c     Read the wind component <fieldname> from the file with identifier
394
c     <fid> and save it in the 3d array <field>. The vertical staggering 
395
c     information is provided in <stagz> and gives the reference to either
396
c     the layer or level field from <input_grid>. A consistency check is
397
c     performed to have an agreement with the grid specified by <xmin,xmax,
398
c     ymin,ymax,dx,dy,nx,ny,nz>.
399
 
400
      use netcdf
401
 
402
      implicit none
403
 
404
c     Declaration of variables and parameters
405
      integer      fid                 ! File identifier
406
      character*80 fieldname           ! Name of the wind field
407
      integer      nx,ny,nz            ! Dimension of fields
408
      real         field(nx,ny,nz)     ! 3d wind field
409
      real         stagz               ! Staggering in the z direction
410
      real         mdv                 ! Missing data flag
411
      real         xmin,xmax,ymin,ymax ! Domain size
412
      real         dx,dy               ! Horizontal resolution
413
      real         time                ! Time
414
      character*80 timecheck           ! Either 'yes' or 'no'
415
 
416
c     Numerical and physical parameters
417
      real        eps                 ! Numerical epsilon
418
      parameter  (eps=0.001)
419
      real        notimecheck         ! 'Flag' for no time check
420
      parameter  (notimecheck=7.26537)
421
 
422
c     Netcdf variables
423
      integer      ierr
424
      character*80 varname
425
      integer      vardim(4)
426
      real         varmin(4),varmax(4)
427
      real         stag(4)
428
      integer      ndim
429
      real         times(10)
430
      integer      ntimes
431
      character*80 cstfile
432
      integer      cstid
433
      real         aklay(200),bklay(200),aklev(200),bklev(200)
434
      real         ps(nx,ny)
435
      integer      dimids (nf90_max_var_dims)
436
      character*80 dimname(nf90_max_var_dims)
437
      integer      varid
438
      integer      cdfid
27 michaesp 439
      real,allocatable, dimension (:)     :: lon,lat,lev
25 michaesp 440
      real,allocatable, dimension (:,:)   :: tmp2
441
      real,allocatable, dimension (:,:,:) :: tmp3
27 michaesp 442
      real,allocatable, dimension (:)     :: aktmp,bktmp
443
      character*80  leveltype
444
      integer       vertical_swap
445
      character*80  units
446
      integer       nakbktmp
447
      integer       dimid
25 michaesp 448
 
449
c     Auxiliary variables
450
      integer      isok
451
      integer      i,j,k
452
      integer      nz1
453
      real         rtime
454
      integer      closear
455
      integer      stat
456
      real         delta
457
 
458
c     Get the number of dimensions -> ndim
459
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
460
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
461
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
462
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
463
 
464
c     Get the dimensions of the arrays -> varid(1:ndim)
465
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
466
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
467
      ierr = nf90_inquire_variable(fid, varid, 
468
     >                                   dimids = dimids(1:ndim))
469
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
470
      do i=1,ndim
471
           ierr = nf90_inquire_dimension(fid, dimids(i), 
472
     >                               name = dimname(i) )
473
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
474
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
475
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
476
      enddo
477
 
478
c     Check whether the list of dimensions is OK
479
      if ( ( dimname(1).ne.'lon'  ).or.
480
     >     ( dimname(2).ne.'lat'  ).or.
481
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
482
     >     ( dimname(4).ne.'time' ) )
483
     >then
484
        print*,' ERROR: the dimensions of the variable are not correct'
485
        print*,'        expected -> lon / lat / lev / time'
486
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
487
        stop
488
      endif
489
 
27 michaesp 490
c     Get dimension of AK,BK
491
      varname = 'nhym'
492
      ierr = NF90_INQ_DIMID(fid,varname,dimid)
493
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
494
      ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
495
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
496
 
25 michaesp 497
c     Allocate memory for reading arrays - depending on <closear>
498
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
499
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
500
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
501
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
502
      allocate(lon(vardim(1)),stat=stat)
503
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
504
      allocate(lat(vardim(2)),stat=stat)
505
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
27 michaesp 506
      allocate(lev(vardim(3)),stat=stat)
507
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
508
      allocate(aktmp(nakbktmp),stat=stat)
509
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
510
      allocate(bktmp(nakbktmp),stat=stat)
511
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
25 michaesp 512
 
27 michaesp 513
c     Get domain boundaries - longitude, latitude, levels
25 michaesp 514
      varname = dimname(1)
515
      ierr = NF90_INQ_VARID(fid,varname,varid)
516
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
517
      ierr = nf90_get_var(fid,varid,lon)
518
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
519
      varname = dimname(2)
520
      ierr = NF90_INQ_VARID(fid,varname,varid)
521
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
522
      ierr = nf90_get_var(fid,varid,lat)
523
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
27 michaesp 524
      varname = dimname(3)
525
      ierr = NF90_INQ_VARID(fid,varname,varid)
526
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
527
      ierr = nf90_get_var(fid,varid,lev)
528
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
25 michaesp 529
 
27 michaesp 530
c     Get ak and bk
531
      varname='hyam'
532
      ierr = NF90_INQ_VARID(fid,varname,varid)
533
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
534
      ierr = nf90_get_var(fid,varid,aktmp)
535
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
536
      varname='hybm'
537
      ierr = NF90_INQ_VARID(fid,varname,varid)
538
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
539
      ierr = nf90_get_var(fid,varid,bktmp)
540
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
541
 
542
c     Check that unit of ak is in hPa - if necessary correct it
543
      varname='hyam'
544
      ierr = NF90_INQ_VARID(fid,varname,varid)
545
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
546
      ierr = nf90_get_att(fid, varid, "units", units)
547
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
548
      if ( units.eq.'Pa' ) then
549
         do k=1,nakbktmp
550
            aktmp(k) = 0.01 * aktmp(k)
551
         enddo
552
      endif
553
 
554
c     Decide whether to swap vertical levels
555
      vertical_swap = 1
556
      if ( (aktmp(1) + bktmp(1) * 1000.).gt.
557
     >     (aktmp(2) + bktmp(2) * 1000.) )
558
     >then
559
          vertical_swap = 0
560
      endif
561
 
25 michaesp 562
c     Read data 
563
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
564
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
565
      ierr = nf90_get_var(fid,varid,tmp3)
566
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
567
 
568
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
569
      if ( vardim(3).eq.1 ) then
570
         do i=1,vardim(1)
571
            do j=1,vardim(2)
572
               do k=1,vardim(3)
573
                  tmp3(i,j,k) = tmp3(i,j,1)
574
               enddo
575
            enddo
576
         enddo
577
      endif
578
 
27 michaesp 579
c     Decide whether to close arrays
25 michaesp 580
      delta = varmax(1)-varmin(1)-360.
581
      if (abs(delta+dx).lt.eps) then
582
          closear = 1
583
      else
584
          closear = 0
585
      endif
586
 
27 michaesp 587
c     Save output array - close array and swap on demand
25 michaesp 588
      do j=1,vardim(2)
589
        do k=1,vardim(3)
590
          do i=1,vardim(1)
27 michaesp 591
             if ( vertical_swap.eq.1 ) then
592
                 field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
593
             else
594
                 field(i,j,k) = tmp3(i,j,k)
595
             endif
25 michaesp 596
          enddo
597
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
598
        enddo
599
      enddo
600
 
601
c     Exit point
602
      return
603
 
604
      end
605
 
606
c     ------------------------------------------------------------
607
c     Close input file
608
c     ------------------------------------------------------------
609
 
610
      subroutine input_close(fid)
611
 
612
c     Close the input file with file identifier <fid>.
613
 
614
      use netcdf
615
 
616
      implicit none
617
 
618
c     Declaration of subroutine parameters
619
      integer fid
620
 
621
c     Auxiliary variables
622
      integer ierr
623
 
624
c     Close file
625
      ierr = NF90_CLOSE(fid)
626
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
627
 
628
      end
629
 
630
c     ------------------------------------------------------------
631
c     Get a list of variables on netCDF file
632
c     ------------------------------------------------------------
633
 
634
      subroutine input_getvars(fid,vnam,nvars)
635
 
636
c     List of variables on netCDF file
637
 
638
      use netcdf
639
 
640
      implicit none
641
 
642
c     Declaration of subroutine parameters
643
      integer      fid
644
      integer      nvars
645
      character*80 vnam(200)
646
 
647
c     Auxiliary variables
648
      integer ierr
649
      integer i
650
      integer nDims,nGlobalAtts,unlimdimid
651
 
652
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
653
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
654
 
655
      do i=1,nVars
656
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
657
      enddo
658
 
659
      end