Subversion Repositories lagranto.ecmwf

Rev

Rev 31 | Rev 36 | 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
33 michaesp 340
c         print*,' WARNING: only subset of vertical levels used...'
27 michaesp 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)
31 michaesp 421
              else
422
                 ak(k) = aktmp(k)
423
                 bk(k) = bktmp(k)
25 michaesp 424
              endif
425
           enddo
426
         elseif (leveltype.eq.'air_pressure' ) then
427
           do k=1,vardim(3)
428
              if ( vertical_swap.eq.1 ) then
429
                 ak(k) = lev(vardim(3)-k+1)
430
                 bk(k) = 0.
431
              else
432
                ak(k) = lev(k)
433
                bk(k) = 0.
434
              endif
435
           enddo
436
         endif
23 michaesp 437
 
21 michaesp 438
      endif
439
 
25 michaesp 440
 
21 michaesp 441
      return
442
 
443
      end
444
 
445
c     ------------------------------------------------------------
446
c     Read wind information
447
c     ------------------------------------------------------------
448
 
449
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
450
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
451
     >                       timecheck)
452
 
453
c     Read the wind component <fieldname> from the file with identifier
454
c     <fid> and save it in the 3d array <field>. The vertical staggering 
455
c     information is provided in <stagz> and gives the reference to either
456
c     the layer or level field from <input_grid>. A consistency check is
457
c     performed to have an agreement with the grid specified by <xmin,xmax,
458
c     ymin,ymax,dx,dy,nx,ny,nz>.
459
 
460
      use netcdf
461
 
462
      implicit none
463
 
464
c     Declaration of variables and parameters
465
      integer      fid                 ! File identifier
466
      character*80 fieldname           ! Name of the wind field
467
      integer      nx,ny,nz            ! Dimension of fields
468
      real         field(nx,ny,nz)     ! 3d wind field
469
      real         stagz               ! Staggering in the z direction
470
      real         mdv                 ! Missing data flag
471
      real         xmin,xmax,ymin,ymax ! Domain size
472
      real         dx,dy               ! Horizontal resolution
473
      real         time                ! Time
474
      character*80 timecheck           ! Either 'yes' or 'no'
475
 
476
c     Numerical and physical parameters
477
      real        eps                 ! Numerical epsilon
478
      parameter  (eps=0.001)
479
      real        notimecheck         ! 'Flag' for no time check
480
      parameter  (notimecheck=7.26537)
481
 
482
c     Netcdf variables
483
      integer      ierr
484
      character*80 varname
485
      integer      vardim(4)
486
      real         varmin(4),varmax(4)
487
      real         stag(4)
488
      integer      ndim
489
      real         times(10)
490
      integer      ntimes
491
      character*80 cstfile
492
      integer      cstid
493
      real         aklay(200),bklay(200),aklev(200),bklev(200)
494
      real         ps(nx,ny)
495
      integer      dimids (nf90_max_var_dims)
496
      character*80 dimname(nf90_max_var_dims)
497
      integer      varid
498
      integer      cdfid
25 michaesp 499
      real,allocatable, dimension (:)     :: lon,lat,lev
21 michaesp 500
      real,allocatable, dimension (:,:)   :: tmp2
501
      real,allocatable, dimension (:,:,:) :: tmp3
25 michaesp 502
      real,allocatable, dimension (:)     :: aktmp,bktmp
503
      character*80  leveltype
504
      integer       vertical_swap
505
      character*80  units
506
      integer       nakbktmp
507
      integer       dimid
21 michaesp 508
 
509
c     Auxiliary variables
510
      integer      isok
511
      integer      i,j,k
512
      integer      nz1
513
      real         rtime
514
      integer      closear
515
      integer      stat
516
      real         delta
517
 
29 michaesp 518
c     Init mdv
519
      mdv = -999.
520
 
21 michaesp 521
c     Get the number of dimensions -> ndim
522
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
523
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
524
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
525
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
526
 
527
c     Get the dimensions of the arrays -> varid(1:ndim)
528
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
529
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
530
      ierr = nf90_inquire_variable(fid, varid, 
531
     >                                   dimids = dimids(1:ndim))
532
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
533
      do i=1,ndim
534
           ierr = nf90_inquire_dimension(fid, dimids(i), 
535
     >                               name = dimname(i) )
536
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
537
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
538
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
539
      enddo
540
 
541
c     Check whether the list of dimensions is OK
542
      if ( ( dimname(1).ne.'lon'  ).or.
543
     >     ( dimname(2).ne.'lat'  ).or.
544
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
545
     >     ( dimname(4).ne.'time' ) )
546
     >then
547
        print*,' ERROR: the dimensions of the variable are not correct'
548
        print*,'        expected -> lon / lat / lev / time'
549
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
550
        stop
551
      endif
552
 
25 michaesp 553
c     Get dimension of AK,BK
554
      varname = 'nhym'
555
      ierr = NF90_INQ_DIMID(fid,varname,dimid)
556
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
557
      ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
558
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
559
 
560
c     Check about the type of vertical levels
561
      varname=dimname(3)
562
      ierr = NF90_INQ_VARID(fid,varname,varid)
563
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
564
      ierr = nf90_get_att(fid, varid, "standard_name", leveltype)
565
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
566
      if ( (leveltype.ne.'hybrid_sigma_pressure').and.
567
     >     (leveltype.ne.'air_pressure'         ) )
568
     >then
569
         print*,' ERROR: input netCDF data must be on hybrid-sigma'
570
         print*,'        or air pressure levels!',trim(leveltype)
571
         stop
572
      endif
573
 
21 michaesp 574
c     Allocate memory for reading arrays - depending on <closear>
575
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
576
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
577
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
578
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
579
      allocate(lon(vardim(1)),stat=stat)
580
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
581
      allocate(lat(vardim(2)),stat=stat)
582
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
25 michaesp 583
      allocate(lev(vardim(3)),stat=stat)
584
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
585
      allocate(aktmp(nakbktmp),stat=stat)
586
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
587
      allocate(bktmp(nakbktmp),stat=stat)
588
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
21 michaesp 589
 
25 michaesp 590
c     Get domain boundaries - longitude, latitude, levels
21 michaesp 591
      varname = dimname(1)
592
      ierr = NF90_INQ_VARID(fid,varname,varid)
593
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
594
      ierr = nf90_get_var(fid,varid,lon)
595
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
596
      varname = dimname(2)
597
      ierr = NF90_INQ_VARID(fid,varname,varid)
598
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
599
      ierr = nf90_get_var(fid,varid,lat)
600
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
25 michaesp 601
      varname = dimname(3)
602
      ierr = NF90_INQ_VARID(fid,varname,varid)
603
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
604
      ierr = nf90_get_var(fid,varid,lev)
605
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
31 michaesp 606
      varmin(1) = lon(1)
607
      varmax(1) = lon( vardim(1) )
608
      varmin(2) = lat(1)
609
      varmax(2) = lat( vardim(2) )
21 michaesp 610
 
25 michaesp 611
c     Get ak and bk
612
      varname='hyam'
613
      ierr = NF90_INQ_VARID(fid,varname,varid)
614
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
615
      ierr = nf90_get_var(fid,varid,aktmp)
616
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
617
      varname='hybm'
618
      ierr = NF90_INQ_VARID(fid,varname,varid)
619
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
620
      ierr = nf90_get_var(fid,varid,bktmp)
621
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
622
 
623
c     Check that unit of ak is in hPa - if necessary correct it
624
      varname='hyam'
625
      ierr = NF90_INQ_VARID(fid,varname,varid)
626
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
627
      ierr = nf90_get_att(fid, varid, "units", units)
628
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
629
      if ( units.eq.'Pa' ) then
630
         do k=1,nakbktmp
631
            aktmp(k) = 0.01 * aktmp(k)
632
         enddo
633
      endif
634
 
635
c     Check that unit of lev is in hPa - if necessary correct it
636
      if ( leveltype.eq.'air_pressure' ) then
637
         varname='lev'
638
         ierr = NF90_INQ_VARID(fid,varname,varid)
639
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
640
         ierr = nf90_get_att(fid, varid, "units", units)
641
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
642
         if ( units.eq.'Pa' ) then
643
            do k=1,vardim(3)
644
               lev(k) = 0.01 * lev(k)
645
            enddo
646
         endif
647
      endif
648
 
649
c     Decide whether to swap vertical levels
650
      vertical_swap = 1
651
      if ( leveltype.eq.'hybrid_sigma_pressure') then
652
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
653
     >       (aktmp(2) + bktmp(2) * 1000.) )
654
     >  then
655
          vertical_swap = 0
656
        endif
657
      elseif ( leveltype.eq.'air_pressure') then
658
        if ( lev(1).gt.lev(2) ) then
659
          vertical_swap = 0
660
        endif
29 michaesp 661
      endif
25 michaesp 662
 
21 michaesp 663
c     Read data 
664
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
665
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
666
      ierr = nf90_get_var(fid,varid,tmp3)
667
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
668
 
669
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
670
      if ( vardim(3).eq.1 ) then
671
         do i=1,vardim(1)
672
            do j=1,vardim(2)
673
               do k=1,vardim(3)
674
                  tmp3(i,j,k) = tmp3(i,j,1)
675
               enddo
676
            enddo
677
         enddo
678
      endif
679
 
25 michaesp 680
c     Decide whether to close arrays
21 michaesp 681
      delta = varmax(1)-varmin(1)-360.
682
      if (abs(delta+dx).lt.eps) then
683
          closear = 1
684
      else
685
          closear = 0
686
      endif
687
 
25 michaesp 688
c     Save output array - close array and swap on demand
21 michaesp 689
      do j=1,vardim(2)
690
        do k=1,vardim(3)
691
          do i=1,vardim(1)
25 michaesp 692
             if ( vertical_swap.eq.1 ) then
693
                 field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
694
             else
695
                 field(i,j,k) = tmp3(i,j,k)
696
             endif
21 michaesp 697
          enddo
698
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
699
        enddo
700
      enddo
701
 
702
c     Exit point
703
      return
704
 
705
      end
706
 
707
c     ------------------------------------------------------------
708
c     Close input file
709
c     ------------------------------------------------------------
710
 
711
      subroutine input_close(fid)
712
 
713
c     Close the input file with file identifier <fid>.
714
 
715
      use netcdf
716
 
717
      implicit none
718
 
719
c     Declaration of subroutine parameters
720
      integer fid
721
 
722
c     Auxiliary variables
723
      integer ierr
724
 
725
c     Close file
726
      ierr = NF90_CLOSE(fid)
727
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
728
 
729
      end
730
 
731
c     ------------------------------------------------------------
732
c     Get a list of variables on netCDF file
733
c     ------------------------------------------------------------
734
 
735
      subroutine input_getvars(fid,vnam,nvars)
736
 
737
c     List of variables on netCDF file
738
 
739
      use netcdf
740
 
741
      implicit none
742
 
743
c     Declaration of subroutine parameters
744
      integer      fid
745
      integer      nvars
746
      character*80 vnam(200)
747
 
748
c     Auxiliary variables
749
      integer ierr
750
      integer i
751
      integer nDims,nGlobalAtts,unlimdimid
752
 
753
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
754
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
755
 
756
      do i=1,nVars
757
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
758
      enddo
759
 
760
      end