Subversion Repositories lagranto.ecmwf

Rev

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

Rev Author Line No. Line
3 michaesp 1
c     ************************************************************
2
c     * This package provides input routines to read the wind    *
3
c     * and other fields from IVE necdf files. The routines are  *
4
c     *                                                          *
5
c     * 1) input_open  : to open a data file                     *
6
c     * 2) input_grid  : to read the grid information, including *
7
c     *                  the vertical levels                     *
8
c     * 3) input_wind  : to read the wind components             *
9
c     * 4) input_close : to close an input file                  *
10
c     *                                                          *
11
c     * The file is characterised by an filename <filename> and  *
12
c     * a file identifier <fid>. The horizontal grid is given by *
13
c     * <xmin,xmax,ymin,ymax,dx,dy,nx,ny> where the pole of the  *
14
c     * rotated grid is given by <pollon,pollat>. The vertical   *
15
c     * grid is characterised by the surface pressure <ps> and   *
16
c     * the pressure at staggered <slev> and unstaggered <ulev>  *
17
c     * levels. The number of levels is given by <nz>. Finally,  *
18
c     * the retrieval of the wind <field> with name <fieldname>  *
19
c     * is characterised by a <time> and a missing data value    *
20
c     * <mdv>.                                                   *
21
c     *                                                          *
22
c     * Author: Michael Sprenger, Autumn 2008                    *
23
c     ************************************************************
24
 
25
c     ------------------------------------------------------------
26
c     Open input file
27
c     ------------------------------------------------------------
28
 
29
      subroutine input_open (fid,filename)
30
 
31
c     Open the input file with filename <filename> and return the
32
c     file identifier <fid> for further reference. 
33
 
21 michaesp 34
      use netcdf
35
 
3 michaesp 36
      implicit none
37
 
38
c     Declaration of subroutine parameters
39
      integer      fid              ! File identifier
40
      character*80 filename         ! Filename
41
 
42
c     Declaration of auxiliary variables
43
      integer      ierr
44
 
21 michaesp 45
c     Open netcdf file
46
      ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
47
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
48
 
3 michaesp 49
c     Exception handling
50
      return
51
 
52
      end
21 michaesp 53
 
3 michaesp 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
 
21 michaesp 73
      use netcdf
74
 
3 michaesp 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
3 michaesp 117
 
21 michaesp 118
c     Auxiliary variables
3 michaesp 119
      integer      ierr       
120
      integer      i,j,k
121
      integer      isok
122
      real         tmp(200)
123
      character*80 varname
124
      real         rtime
21 michaesp 125
      integer      varid
126
      integer      cdfid
127
      integer      stat
128
      real         delta
129
      integer      closear
130
      real         maxps,minps
3 michaesp 131
 
25 michaesp 132
c     ---- Read data from netCDF file as they are ---------------------
133
 
134
c     Set file identifier
3 michaesp 135
      if (fid.lt.0) then
21 michaesp 136
        cdfid = -fid
137
      else 
138
        cdfid = fid
139
      endif
3 michaesp 140
 
21 michaesp 141
c     Special handling if 3D pressure is
142
      if ( fieldname.eq.'P' ) then
143
         fieldname = 'U'
144
      endif
3 michaesp 145
 
21 michaesp 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
3 michaesp 156
 
21 michaesp 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
3 michaesp 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
3 michaesp 188
 
29 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   ***'
19 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)
19 michaesp 249
 
21 michaesp 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
19 michaesp 261
 
29 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
29 michaesp 278
      if ( leveltype.eq.'hybrid_sigma_pressure') then
279
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
25 michaesp 280
     >       (aktmp(2) + bktmp(2) * 1000.) )
29 michaesp 281
     >  then
25 michaesp 282
          vertical_swap = 0
29 michaesp 283
        endif
284
      elseif ( leveltype.eq.'air_pressure') then
285
        if ( lev(1).gt.lev(2) ) then
286
          vertical_swap = 0
287
        endif
25 michaesp 288
      endif
29 michaesp 289
c      print*,' Vertical SWAP P -> ', vertical_swap
25 michaesp 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
3 michaesp 300
               isok = 1
301
               rtime = times(i)
21 michaesp 302
        elseif (timecheck.eq.'no') then
3 michaesp 303
               isok = 1
304
               rtime = times(1)
21 michaesp 305
        endif
306
      enddo
307
      if ( isok.eq.0 ) then
308
         print*,' ERROR: time ',rtime,' not found on netCDF file' 
309
         stop
310
      endif
3 michaesp 311
 
21 michaesp 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 --------------------------------
3 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
3 michaesp 370
 
21 michaesp 371
c     Save the output arrays (if fid>0) - close arrays on request
372
      if ( fid.gt.0 ) then
3 michaesp 373
 
25 michaesp 374
c        Calculate layer pressures
29 michaesp 375
         if (leveltype.eq.'hybrid_sigma_pressure' ) then
376
            do i=1,vardim(1)
25 michaesp 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
29 michaesp 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
25 michaesp 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)
3 michaesp 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
3 michaesp 413
         enddo
23 michaesp 414
 
25 michaesp 415
c        Get AK,BK - vertical swap on demand
29 michaesp 416
         if ( leveltype.eq.'hybrid_sigma_pressure' ) then
417
           do k=1,vardim(3)
25 michaesp 418
              if ( vertical_swap.eq.1 ) then
419
                 ak(k) = aktmp(vardim(3)-k+1)
420
                 bk(k) = bktmp(vardim(3)-k+1)
33 michaesp 421
              else
422
                 ak(k) = aktmp(k)
423
                 bk(k) = bktmp(k)
25 michaesp 424
              endif
29 michaesp 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
3 michaesp 439
 
25 michaesp 440
 
3 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
 
21 michaesp 460
      use netcdf
461
 
3 michaesp 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)
21 michaesp 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
3 michaesp 508
 
509
c     Auxiliary variables
510
      integer      isok
511
      integer      i,j,k
512
      integer      nz1
513
      real         rtime
21 michaesp 514
      integer      closear
515
      integer      stat
516
      real         delta
36 michaesp 517
      integer      pressure
3 michaesp 518
 
29 michaesp 519
c     Init mdv
520
      mdv = -999.
36 michaesp 521
 
522
c     Special handling if 3D pressure is
523
      if ( fieldname.eq.'P' ) then
524
         fieldname = 'U'
525
         pressure  = 1
526
      else
527
         pressure = 0
528
      endif
29 michaesp 529
 
21 michaesp 530
c     Get the number of dimensions -> ndim
531
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
532
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
533
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
534
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
19 michaesp 535
 
21 michaesp 536
c     Get the dimensions of the arrays -> varid(1:ndim)
537
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
538
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
539
      ierr = nf90_inquire_variable(fid, varid, 
540
     >                                   dimids = dimids(1:ndim))
541
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
542
      do i=1,ndim
543
           ierr = nf90_inquire_dimension(fid, dimids(i), 
544
     >                               name = dimname(i) )
545
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
546
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
547
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
3 michaesp 548
      enddo
549
 
21 michaesp 550
c     Check whether the list of dimensions is OK
551
      if ( ( dimname(1).ne.'lon'  ).or.
552
     >     ( dimname(2).ne.'lat'  ).or.
553
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
554
     >     ( dimname(4).ne.'time' ) )
555
     >then
556
        print*,' ERROR: the dimensions of the variable are not correct'
557
        print*,'        expected -> lon / lat / lev / time'
558
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
559
        stop
560
      endif
3 michaesp 561
 
25 michaesp 562
c     Get dimension of AK,BK
563
      varname = 'nhym'
564
      ierr = NF90_INQ_DIMID(fid,varname,dimid)
565
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
566
      ierr = nf90_inquire_dimension(fid, dimid,len=nakbktmp)
567
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
568
 
29 michaesp 569
c     Check about the type of vertical levels
570
      varname=dimname(3)
571
      ierr = NF90_INQ_VARID(fid,varname,varid)
572
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
573
      ierr = nf90_get_att(fid, varid, "standard_name", leveltype)
574
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
575
      if ( (leveltype.ne.'hybrid_sigma_pressure').and.
576
     >     (leveltype.ne.'air_pressure'         ) )
577
     >then
578
         print*,' ERROR: input netCDF data must be on hybrid-sigma'
579
         print*,'        or air pressure levels!',trim(leveltype)
580
         stop
581
      endif
582
 
21 michaesp 583
c     Allocate memory for reading arrays - depending on <closear>
584
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
585
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
586
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
587
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
588
      allocate(lon(vardim(1)),stat=stat)
589
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
590
      allocate(lat(vardim(2)),stat=stat)
591
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
25 michaesp 592
      allocate(lev(vardim(3)),stat=stat)
593
      if (stat.ne.0) print*,'*** error allocating array lev     ***'
594
      allocate(aktmp(nakbktmp),stat=stat)
595
      if (stat.ne.0) print*,'*** error allocating array aktmp   ***'
596
      allocate(bktmp(nakbktmp),stat=stat)
597
      if (stat.ne.0) print*,'*** error allocating array bktmp   ***'
3 michaesp 598
 
25 michaesp 599
c     Get domain boundaries - longitude, latitude, levels
21 michaesp 600
      varname = dimname(1)
601
      ierr = NF90_INQ_VARID(fid,varname,varid)
602
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
603
      ierr = nf90_get_var(fid,varid,lon)
604
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
605
      varname = dimname(2)
606
      ierr = NF90_INQ_VARID(fid,varname,varid)
607
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
608
      ierr = nf90_get_var(fid,varid,lat)
609
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
25 michaesp 610
      varname = dimname(3)
611
      ierr = NF90_INQ_VARID(fid,varname,varid)
612
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
613
      ierr = nf90_get_var(fid,varid,lev)
614
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
31 michaesp 615
      varmin(1) = lon(1)
616
      varmax(1) = lon( vardim(1) )
617
      varmin(2) = lat(1)
618
      varmax(2) = lat( vardim(2) )
3 michaesp 619
 
25 michaesp 620
c     Get ak and bk
621
      varname='hyam'
622
      ierr = NF90_INQ_VARID(fid,varname,varid)
623
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
624
      ierr = nf90_get_var(fid,varid,aktmp)
625
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
626
      varname='hybm'
627
      ierr = NF90_INQ_VARID(fid,varname,varid)
628
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
629
      ierr = nf90_get_var(fid,varid,bktmp)
630
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
631
 
632
c     Check that unit of ak is in hPa - if necessary correct it
633
      varname='hyam'
634
      ierr = NF90_INQ_VARID(fid,varname,varid)
635
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
636
      ierr = nf90_get_att(fid, varid, "units", units)
637
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
638
      if ( units.eq.'Pa' ) then
639
         do k=1,nakbktmp
640
            aktmp(k) = 0.01 * aktmp(k)
641
         enddo
642
      endif
643
 
29 michaesp 644
c     Check that unit of lev is in hPa - if necessary correct it
645
      if ( leveltype.eq.'air_pressure' ) then
646
         varname='lev'
647
         ierr = NF90_INQ_VARID(fid,varname,varid)
648
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
649
         ierr = nf90_get_att(fid, varid, "units", units)
650
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
651
         if ( units.eq.'Pa' ) then
652
            do k=1,vardim(3)
653
               lev(k) = 0.01 * lev(k)
654
            enddo
655
         endif
656
      endif
657
 
25 michaesp 658
c     Decide whether to swap vertical levels
659
      vertical_swap = 1
29 michaesp 660
      if ( leveltype.eq.'hybrid_sigma_pressure') then
661
        if ( (aktmp(1) + bktmp(1) * 1000.).gt.
662
     >       (aktmp(2) + bktmp(2) * 1000.) )
663
     >  then
25 michaesp 664
          vertical_swap = 0
29 michaesp 665
        endif
666
      elseif ( leveltype.eq.'air_pressure') then
667
        if ( lev(1).gt.lev(2) ) then
668
          vertical_swap = 0
669
        endif
25 michaesp 670
      endif
671
 
36 michaesp 672
c     Read data /or set 3D pressure field) 
673
      if ( pressure.eq.0 ) then
674
         ierr = NF90_INQ_VARID(fid,fieldname,varid)
675
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
676
         ierr = nf90_get_var(fid,varid,tmp3)
677
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
678
      else
679
         if (leveltype.eq.'hybrid_sigma_pressure' ) then
680
            do i=1,vardim(1)
681
              do j=1,vardim(2)
682
                 do k=1,vardim(3)
683
                  tmp3(i,j,k)=aktmp(k)+bktmp(k)*tmp2(i,j)
684
                 enddo
685
              enddo
686
           enddo
687
         elseif (leveltype.eq.'air_pressure' ) then
688
           do i=1,vardim(1)
689
              do j=1,vardim(2)
690
                 do k=1,vardim(3)
691
                  tmp3(i,j,k)=lev(k)
692
                 enddo
693
              enddo
694
           enddo
695
         endif
696
      endif
21 michaesp 697
 
3 michaesp 698
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
699
      if ( vardim(3).eq.1 ) then
21 michaesp 700
         do i=1,vardim(1)
701
            do j=1,vardim(2)
702
               do k=1,vardim(3)
703
                  tmp3(i,j,k) = tmp3(i,j,1)
3 michaesp 704
               enddo
705
            enddo
706
         enddo
707
      endif
708
 
25 michaesp 709
c     Decide whether to close arrays
21 michaesp 710
      delta = varmax(1)-varmin(1)-360.
711
      if (abs(delta+dx).lt.eps) then
712
          closear = 1
713
      else
714
          closear = 0
715
      endif
3 michaesp 716
 
25 michaesp 717
c     Save output array - close array and swap on demand
21 michaesp 718
      do j=1,vardim(2)
719
        do k=1,vardim(3)
720
          do i=1,vardim(1)
25 michaesp 721
             if ( vertical_swap.eq.1 ) then
722
                 field(i,j,k) = tmp3(i,j,vardim(3)-k+1)
723
             else
724
                 field(i,j,k) = tmp3(i,j,k)
725
             endif
21 michaesp 726
          enddo
727
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
728
        enddo
729
      enddo
730
 
731
c     Exit point
3 michaesp 732
      return
21 michaesp 733
 
3 michaesp 734
      end
735
 
736
c     ------------------------------------------------------------
737
c     Close input file
738
c     ------------------------------------------------------------
739
 
740
      subroutine input_close(fid)
741
 
742
c     Close the input file with file identifier <fid>.
743
 
21 michaesp 744
      use netcdf
745
 
3 michaesp 746
      implicit none
747
 
748
c     Declaration of subroutine parameters
749
      integer fid
750
 
751
c     Auxiliary variables
752
      integer ierr
753
 
754
c     Close file
21 michaesp 755
      ierr = NF90_CLOSE(fid)
756
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
757
 
758
      end
759
 
760
c     ------------------------------------------------------------
761
c     Get a list of variables on netCDF file
762
c     ------------------------------------------------------------
763
 
764
      subroutine input_getvars(fid,vnam,nvars)
765
 
766
c     List of variables on netCDF file
767
 
768
      use netcdf
769
 
770
      implicit none
771
 
772
c     Declaration of subroutine parameters
773
      integer      fid
774
      integer      nvars
775
      character*80 vnam(200)
776
 
777
c     Auxiliary variables
778
      integer ierr
779
      integer i
780
      integer nDims,nGlobalAtts,unlimdimid
781
 
782
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
783
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
784
 
785
      do i=1,nVars
786
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
787
      enddo
3 michaesp 788
 
789
      end