Subversion Repositories lagranto.ecmwf

Rev

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

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