Subversion Repositories lagranto.ecmwf

Rev

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

Rev Author Line No. Line
21 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)
25 michaesp 105
      integer        dimids (nf90_max_var_dims),dimid
21 michaesp 106
      character*80   dimname(nf90_max_var_dims)
25 michaesp 107
      character*80   stdname
21 michaesp 108
      real,allocatable, dimension (:)     :: lon,lat,lev
109
      real,allocatable, dimension (:)     :: times
110
      real,allocatable, dimension (:,:)   :: tmp2
111
      real,allocatable, dimension (:,:,:) :: tmp3
23 michaesp 112
      real,allocatable, dimension (:)     :: aktmp,bktmp
21 michaesp 113
      character*80  units
25 michaesp 114
      character*80  leveltype
115
      integer       nakbktmp
116
      integer       vertical_swap
21 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
 
25 michaesp 132
c     ---- Read data from netCDF file as they are ---------------------
133
 
134
c     Set file identifier
21 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
 
25 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
 
21 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
 
25 michaesp 189
c     Check about the type of vertical levels
190
      varname=dimname(3)
191
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
192
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
193
      ierr = nf90_get_att(cdfid, varid, "standard_name", leveltype)
194
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
195
      if ( (leveltype.ne.'hybrid_sigma_pressure').and.
196
     >     (leveltype.ne.'air_pressure'         ) )
197
     >then
198
         print*,' ERROR: input netCDF data must be on hybrid-sigma'
199
         print*,'        or air pressure levels!',trim(leveltype)
200
         stop
201
      endif
202
 
21 michaesp 203
c     Allocate memory for reading arrays
204
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
205
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
206
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
207
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
208
      allocate(lon(vardim(1)),stat=stat)
209
      if (stat.ne.0) print*,'*** error allocating array lon     ***' 
210
      allocate(lat(vardim(2)),stat=stat)
211
      if (stat.ne.0) print*,'*** error allocating array lat     ***' 
25 michaesp 212
      allocate(lev(vardim(3)),stat=stat)
213
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
21 michaesp 214
      allocate(times(vardim(4)),stat=stat)
215
      if (stat.ne.0) print*,'*** error allocating array times   ***'
25 michaesp 216
      allocate(aktmp(nakbktmp),stat=stat)
23 michaesp 217
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
25 michaesp 218
      allocate(bktmp(nakbktmp),stat=stat)
23 michaesp 219
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
21 michaesp 220
 
25 michaesp 221
c     Get domain longitudes, latitudes and levels
21 michaesp 222
      varname = dimname(1)
223
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
224
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
225
      ierr = nf90_get_var(cdfid,varid,lon)
226
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
227
      varname = dimname(2)
228
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
229
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
230
      ierr = nf90_get_var(cdfid,varid,lat)
231
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
25 michaesp 232
      varname = dimname(3)
233
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
234
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
235
      ierr = nf90_get_var(cdfid,varid,lev)
236
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
21 michaesp 237
 
238
c     Get ak and bk
239
      varname='hyam'
240
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
241
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
23 michaesp 242
      ierr = nf90_get_var(cdfid,varid,aktmp)
21 michaesp 243
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
244
      varname='hybm'
245
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
246
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
23 michaesp 247
      ierr = nf90_get_var(cdfid,varid,bktmp)
21 michaesp 248
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
249
 
250
c     Check that unit of ak is in hPa - if necessary correct it
251
      varname='hyam'
252
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
253
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
254
      ierr = nf90_get_att(cdfid, varid, "units", units)
255
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
256
      if ( units.eq.'Pa' ) then
25 michaesp 257
         do k=1,nakbktmp
23 michaesp 258
            aktmp(k) = 0.01 * aktmp(k)
21 michaesp 259
         enddo
260
      endif
261
 
25 michaesp 262
c     Check that unit of lev is in hPa - if necessary correct it
263
      if ( leveltype.eq.'air_pressure' ) then
264
         varname='lev'
265
         ierr = NF90_INQ_VARID(cdfid,varname,varid)
266
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
267
         ierr = nf90_get_att(cdfid, varid, "units", units)
268
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
269
         if ( units.eq.'Pa' ) then
270
            do k=1,vardim(3)
271
               lev(k) = 0.01 * lev(k)
272
            enddo
273
         endif
274
      endif
275
 
27 michaesp 276
c     Decide whether to swap vertical levels - highest pressure at index 1
25 michaesp 277
      vertical_swap = 1
278
      if ( leveltype.eq.'hybrid_sigma_pressure') then
279
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
280
     >       (aktmp(2) + bktmp(2) * 1000.) )
281
     >  then
282
          vertical_swap = 0
283
        endif
284
      elseif ( leveltype.eq.'air_pressure') then
285
        if ( lev(1).gt.lev(2) ) then
286
          vertical_swap = 0
287
        endif
288
      endif
289
c      print*,' Vertical SWAP P -> ', vertical_swap
290
 
21 michaesp 291
c     Get time information (check if time is correct)
292
      varname = 'time'
293
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
294
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
295
      ierr = nf90_get_var(cdfid,varid,times)
296
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
297
      isok=0
298
      do i=1,vardim(4)
299
        if (abs(time-times(i)).lt.eps) then
300
               isok = 1
301
               rtime = times(i)
302
        elseif (timecheck.eq.'no') then
303
               isok = 1
304
               rtime = times(1)
305
        endif
306
      enddo
307
      if ( isok.eq.0 ) then
308
         print*,' ERROR: time ',rtime,' not found on netCDF file' 
309
         stop
310
      endif
311
 
312
c     Read surface pressure
313
      varname='PS'
314
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
315
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
316
      ierr = nf90_get_var(cdfid,varid,tmp2)
317
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
318
 
319
c     Check that surface pressure is in hPa
25 michaesp 320
      maxps = -1.e19
321
      minps =  1.e19
21 michaesp 322
      do i=1,vardim(1)
323
        do j=1,vardim(2)
324
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
325
             if (tmp2(i,j).lt.minps) minps = tmp2(i,j)
326
        enddo
327
      enddo
328
      if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then
329
         print*,' ERROR: surface pressre PS must be in hPa'
330
         print*,'       ',maxps,minps
331
         stop
332
      endif
333
 
25 michaesp 334
c     ---- Define output of subroutine --------------------------------
21 michaesp 335
 
27 michaesp 336
c     If not full list of vertical levels, reduce AK,BK arrays
337
      if ( (leveltype.eq.'hybrid_sigma_pressure').and.
338
     >     (nakbktmp.ne.vardim(3) ) )
339
     >then
340
         print*,' WARNING: only subset of vertical levels used...'
341
         do k=1,vardim(3)
342
            if ( vertical_swap.eq.1 ) then
343
               aktmp(k) = aktmp( k+nakbktmp-vardim(3) )
344
               bktmp(k) = bktmp( k+nakbktmp-vardim(3) )
345
            endif
346
         enddo
347
      endif
348
 
21 michaesp 349
c     Set the grid dimensions and constants
350
      nx      = vardim(1)
351
      ny      = vardim(2)
352
      nz      = vardim(3)
353
      xmin    = lon(1)
354
      ymin    = lat(1)
355
      xmax    = lon(nx)
356
      ymax    = lat(ny)
357
      dx      = (xmax-xmin)/real(nx-1)
358
      dy      = (ymax-ymin)/real(ny-1)
359
      pollon  = 0.
360
      pollat  = 90.
361
      stagz   = 0.
362
      delta   = xmax-xmin-360.
363
      if (abs(delta+dx).lt.eps) then
364
          xmax    = xmax + dx
365
          nx      = nx + 1
366
          closear = 1
367
      else
368
          closear = 0
369
      endif
370
 
371
c     Save the output arrays (if fid>0) - close arrays on request
372
      if ( fid.gt.0 ) then
373
 
25 michaesp 374
c        Calculate layer pressures
375
         if (leveltype.eq.'hybrid_sigma_pressure' ) then
376
            do i=1,vardim(1)
377
              do j=1,vardim(2)
378
                 do k=1,vardim(3)
379
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
380
                 enddo
381
              enddo
382
           enddo
383
         elseif (leveltype.eq.'air_pressure' ) then
384
           do i=1,vardim(1)
385
              do j=1,vardim(2)
386
                 do k=1,vardim(3)
387
                  tmp3(i,j,k)=lev(k)
388
                 enddo
389
              enddo
390
           enddo
391
         endif
392
 
393
c        Get PS - close array on demand
21 michaesp 394
         do j=1,vardim(2)
395
           do i=1,vardim(1)
396
             ps(i,j) = tmp2(i,j)
397
           enddo
398
           if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)
399
         enddo
400
 
25 michaesp 401
c        Get P3 - close array on demand + vertical swap
21 michaesp 402
         do j=1,vardim(2)
403
           do k=1,vardim(3)
404
             do i=1,vardim(1)
25 michaesp 405
               if ( vertical_swap.eq.1 ) then
406
                  p3(i,j,k) = tmp3(i,j,vardim(3)-k+1)
407
               else
408
                  p3(i,j,k) = tmp3(i,j,k)
409
               endif
21 michaesp 410
             enddo
411
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
412
           enddo
413
         enddo
23 michaesp 414
 
25 michaesp 415
c        Get AK,BK - vertical swap on demand
416
         if ( leveltype.eq.'hybrid_sigma_pressure' ) then
417
           do k=1,vardim(3)
418
              if ( vertical_swap.eq.1 ) then
419
                 ak(k) = aktmp(vardim(3)-k+1)
420
                 bk(k) = bktmp(vardim(3)-k+1)
421
              endif
422
           enddo
423
         elseif (leveltype.eq.'air_pressure' ) then
424
           do k=1,vardim(3)
425
              if ( vertical_swap.eq.1 ) then
426
                 ak(k) = lev(vardim(3)-k+1)
427
                 bk(k) = 0.
428
              else
429
                ak(k) = lev(k)
430
                bk(k) = 0.
431
              endif
432
           enddo
433
         endif
23 michaesp 434
 
21 michaesp 435
      endif
436
 
25 michaesp 437
 
21 michaesp 438
      return
439
 
440
      end
441
 
442
c     ------------------------------------------------------------
443
c     Read wind information
444
c     ------------------------------------------------------------
445
 
446
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
447
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
448
     >                       timecheck)
449
 
450
c     Read the wind component <fieldname> from the file with identifier
451
c     <fid> and save it in the 3d array <field>. The vertical staggering 
452
c     information is provided in <stagz> and gives the reference to either
453
c     the layer or level field from <input_grid>. A consistency check is
454
c     performed to have an agreement with the grid specified by <xmin,xmax,
455
c     ymin,ymax,dx,dy,nx,ny,nz>.
456
 
457
      use netcdf
458
 
459
      implicit none
460
 
461
c     Declaration of variables and parameters
462
      integer      fid                 ! File identifier
463
      character*80 fieldname           ! Name of the wind field
464
      integer      nx,ny,nz            ! Dimension of fields
465
      real         field(nx,ny,nz)     ! 3d wind field
466
      real         stagz               ! Staggering in the z direction
467
      real         mdv                 ! Missing data flag
468
      real         xmin,xmax,ymin,ymax ! Domain size
469
      real         dx,dy               ! Horizontal resolution
470
      real         time                ! Time
471
      character*80 timecheck           ! Either 'yes' or 'no'
472
 
473
c     Numerical and physical parameters
474
      real        eps                 ! Numerical epsilon
475
      parameter  (eps=0.001)
476
      real        notimecheck         ! 'Flag' for no time check
477
      parameter  (notimecheck=7.26537)
478
 
479
c     Netcdf variables
480
      integer      ierr
481
      character*80 varname
482
      integer      vardim(4)
483
      real         varmin(4),varmax(4)
484
      real         stag(4)
485
      integer      ndim
486
      real         times(10)
487
      integer      ntimes
488
      character*80 cstfile
489
      integer      cstid
490
      real         aklay(200),bklay(200),aklev(200),bklev(200)
491
      real         ps(nx,ny)
492
      integer      dimids (nf90_max_var_dims)
493
      character*80 dimname(nf90_max_var_dims)
494
      integer      varid
495
      integer      cdfid
25 michaesp 496
      real,allocatable, dimension (:)     :: lon,lat,lev
21 michaesp 497
      real,allocatable, dimension (:,:)   :: tmp2
498
      real,allocatable, dimension (:,:,:) :: tmp3
25 michaesp 499
      real,allocatable, dimension (:)     :: aktmp,bktmp
500
      character*80  leveltype
501
      integer       vertical_swap
502
      character*80  units
503
      integer       nakbktmp
504
      integer       dimid
21 michaesp 505
 
506
c     Auxiliary variables
507
      integer      isok
508
      integer      i,j,k
509
      integer      nz1
510
      real         rtime
511
      integer      closear
512
      integer      stat
513
      real         delta
514
 
29 michaesp 515
c     Init mdv
516
      mdv = -999.
517
 
21 michaesp 518
c     Get the number of dimensions -> ndim
519
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
520
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
521
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
522
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
523
 
524
c     Get the dimensions of the arrays -> varid(1:ndim)
525
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
526
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
527
      ierr = nf90_inquire_variable(fid, varid, 
528
     >                                   dimids = dimids(1:ndim))
529
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
530
      do i=1,ndim
531
           ierr = nf90_inquire_dimension(fid, dimids(i), 
532
     >                               name = dimname(i) )
533
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
534
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
535
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
536
      enddo
537
 
538
c     Check whether the list of dimensions is OK
539
      if ( ( dimname(1).ne.'lon'  ).or.
540
     >     ( dimname(2).ne.'lat'  ).or.
541
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
542
     >     ( dimname(4).ne.'time' ) )
543
     >then
544
        print*,' ERROR: the dimensions of the variable are not correct'
545
        print*,'        expected -> lon / lat / lev / time'
546
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
547
        stop
548
      endif
549
 
25 michaesp 550
c     Get dimension of AK,BK
551
      varname = 'nhym'
552
      ierr = NF90_INQ_DIMID(fid,varname,dimid)
553
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
554
      ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
555
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
556
 
557
c     Check about the type of vertical levels
558
      varname=dimname(3)
559
      ierr = NF90_INQ_VARID(fid,varname,varid)
560
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
561
      ierr = nf90_get_att(fid, varid, "standard_name", leveltype)
562
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
563
      if ( (leveltype.ne.'hybrid_sigma_pressure').and.
564
     >     (leveltype.ne.'air_pressure'         ) )
565
     >then
566
         print*,' ERROR: input netCDF data must be on hybrid-sigma'
567
         print*,'        or air pressure levels!',trim(leveltype)
568
         stop
569
      endif
570
 
21 michaesp 571
c     Allocate memory for reading arrays - depending on <closear>
572
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
573
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
574
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
575
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
576
      allocate(lon(vardim(1)),stat=stat)
577
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
578
      allocate(lat(vardim(2)),stat=stat)
579
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
25 michaesp 580
      allocate(lev(vardim(3)),stat=stat)
581
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
582
      allocate(aktmp(nakbktmp),stat=stat)
583
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
584
      allocate(bktmp(nakbktmp),stat=stat)
585
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
21 michaesp 586
 
25 michaesp 587
c     Get domain boundaries - longitude, latitude, levels
21 michaesp 588
      varname = dimname(1)
589
      ierr = NF90_INQ_VARID(fid,varname,varid)
590
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
591
      ierr = nf90_get_var(fid,varid,lon)
592
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
593
      varname = dimname(2)
594
      ierr = NF90_INQ_VARID(fid,varname,varid)
595
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
596
      ierr = nf90_get_var(fid,varid,lat)
597
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
25 michaesp 598
      varname = dimname(3)
599
      ierr = NF90_INQ_VARID(fid,varname,varid)
600
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
601
      ierr = nf90_get_var(fid,varid,lev)
602
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
21 michaesp 603
 
25 michaesp 604
c     Get ak and bk
605
      varname='hyam'
606
      ierr = NF90_INQ_VARID(fid,varname,varid)
607
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
608
      ierr = nf90_get_var(fid,varid,aktmp)
609
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
610
      varname='hybm'
611
      ierr = NF90_INQ_VARID(fid,varname,varid)
612
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
613
      ierr = nf90_get_var(fid,varid,bktmp)
614
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
615
 
616
c     Check that unit of ak is in hPa - if necessary correct it
617
      varname='hyam'
618
      ierr = NF90_INQ_VARID(fid,varname,varid)
619
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
620
      ierr = nf90_get_att(fid, varid, "units", units)
621
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
622
      if ( units.eq.'Pa' ) then
623
         do k=1,nakbktmp
624
            aktmp(k) = 0.01 * aktmp(k)
625
         enddo
626
      endif
627
 
628
c     Check that unit of lev is in hPa - if necessary correct it
629
      if ( leveltype.eq.'air_pressure' ) then
630
         varname='lev'
631
         ierr = NF90_INQ_VARID(fid,varname,varid)
632
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
633
         ierr = nf90_get_att(fid, varid, "units", units)
634
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
635
         if ( units.eq.'Pa' ) then
636
            do k=1,vardim(3)
637
               lev(k) = 0.01 * lev(k)
638
            enddo
639
         endif
640
      endif
641
 
642
c     Decide whether to swap vertical levels
643
      vertical_swap = 1
644
      if ( leveltype.eq.'hybrid_sigma_pressure') then
645
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
646
     >       (aktmp(2) + bktmp(2) * 1000.) )
647
     >  then
648
          vertical_swap = 0
649
        endif
650
      elseif ( leveltype.eq.'air_pressure') then
651
        if ( lev(1).gt.lev(2) ) then
652
          vertical_swap = 0
653
        endif
29 michaesp 654
      endif
25 michaesp 655
 
21 michaesp 656
c     Read data 
657
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
658
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
659
      ierr = nf90_get_var(fid,varid,tmp3)
660
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
661
 
662
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
663
      if ( vardim(3).eq.1 ) then
664
         do i=1,vardim(1)
665
            do j=1,vardim(2)
666
               do k=1,vardim(3)
667
                  tmp3(i,j,k) = tmp3(i,j,1)
668
               enddo
669
            enddo
670
         enddo
671
      endif
672
 
25 michaesp 673
c     Decide whether to close arrays
21 michaesp 674
      delta = varmax(1)-varmin(1)-360.
675
      if (abs(delta+dx).lt.eps) then
676
          closear = 1
677
      else
678
          closear = 0
679
      endif
680
 
25 michaesp 681
c     Save output array - close array and swap on demand
21 michaesp 682
      do j=1,vardim(2)
683
        do k=1,vardim(3)
684
          do i=1,vardim(1)
25 michaesp 685
             if ( vertical_swap.eq.1 ) then
686
                 field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
687
             else
688
                 field(i,j,k) = tmp3(i,j,k)
689
             endif
21 michaesp 690
          enddo
691
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
692
        enddo
693
      enddo
694
 
695
c     Exit point
696
      return
697
 
698
      end
699
 
700
c     ------------------------------------------------------------
701
c     Close input file
702
c     ------------------------------------------------------------
703
 
704
      subroutine input_close(fid)
705
 
706
c     Close the input file with file identifier <fid>.
707
 
708
      use netcdf
709
 
710
      implicit none
711
 
712
c     Declaration of subroutine parameters
713
      integer fid
714
 
715
c     Auxiliary variables
716
      integer ierr
717
 
718
c     Close file
719
      ierr = NF90_CLOSE(fid)
720
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
721
 
722
      end
723
 
724
c     ------------------------------------------------------------
725
c     Get a list of variables on netCDF file
726
c     ------------------------------------------------------------
727
 
728
      subroutine input_getvars(fid,vnam,nvars)
729
 
730
c     List of variables on netCDF file
731
 
732
      use netcdf
733
 
734
      implicit none
735
 
736
c     Declaration of subroutine parameters
737
      integer      fid
738
      integer      nvars
739
      character*80 vnam(200)
740
 
741
c     Auxiliary variables
742
      integer ierr
743
      integer i
744
      integer nDims,nGlobalAtts,unlimdimid
745
 
746
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
747
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
748
 
749
      do i=1,nVars
750
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
751
      enddo
752
 
753
      end