Subversion Repositories lagranto.ecmwf

Rev

Rev 36 | Rev 46 | 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
 
34
      implicit none
35
 
36
c     Declaration of subroutine parameters
37
      integer      fid              ! File identifier
38
      character*80 filename         ! Filename
39
 
40
c     Declaration of auxiliary variables
41
      integer      ierr
42
 
44 michaesp 43
c     Open IVE netcdf file
44
      call cdfopn(filename,fid,ierr)
45
      if (ierr.ne.0) goto 900
46
 
3 michaesp 47
c     Exception handling
48
      return
44 michaesp 49
 
50
 900  print*,'Cannot open input file ',trim(filename)
51
      stop
3 michaesp 52
 
53
      end
44 michaesp 54
 
55
 
3 michaesp 56
c     ------------------------------------------------------------
57
c     Read information about the grid
58
c     ------------------------------------------------------------
59
 
60
      subroutine input_grid 
61
     >                   (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
62
     >                    time,pollon,pollat,p3,ps,nz,ak,bk,stagz,
63
     >                    timecheck)
64
 
65
c     Read grid information at <time> from file with identifier <fid>. 
66
c     The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>
67
c     with pole position at <pollon,pollat> and grid dimension <nx,ny>.
68
c     The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either
69
c     on the staggered or unstaggered grid (with <stagz> as the flag).
70
c     The surface pressure is given in <ps(nx,ny)>. If <fid> is negative, 
71
c     only the grid dimensions and grid parameters (xmin...pollat,nz) are 
72
c     determined and returned (this is needed for dynamical allocation of 
73
c     memory).
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
44 michaesp 103
      real         times(10)
104
      integer      ntimes
105
      real         aklay(nz),bklay(nz),aklev(nz),bklev(nz)
3 michaesp 106
      integer      nvars
107
      character*80 vars(100)
108
 
44 michaesp 109
c     Auxiliary varaibles
3 michaesp 110
      integer      ierr       
111
      integer      i,j,k
112
      integer      isok
113
      real         tmp(200)
114
      character*80 varname
115
      real         rtime
44 michaesp 116
      integer      is2d
117
      integer      plev
3 michaesp 118
 
44 michaesp 119
c     Init the flag for 2D variables - assume a 3D field
120
      is2d = 0
25 michaesp 121
 
44 michaesp 122
c     Init the flag for pressure levels (PS is not needed)
123
      plev = 0
124
 
125
c     Inquire dimensions and grid constants if <fid> is negative
3 michaesp 126
      if (fid.lt.0) then
127
 
44 michaesp 128
c        Get grid info for the selected variable
129
         if ( fieldname.eq.'PLEV' ) then
130
            varname = 'PS'
131
            stagz   = 0.
132
            call getdef(-fid,varname,ndim,mdv,vardim,
133
     >                  varmin,varmax,stag,ierr)
134
            if (ierr.ne.0) goto 900
135
            call getcfn(-fid,cstfile,ierr)
136
            if (ierr.ne.0) goto 903
137
            call cdfopn(cstfile,cstid,ierr)
138
            if (ierr.ne.0) goto 903
139
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
140
            if (ierr.ne.0) goto 903
141
            call clscdf(cstid,ierr)
142
            if (ierr.ne.0) goto 903
143
 
144
         elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then       
145
            varname = 'PS'
146
            stagz   = -0.5
147
            call getdef(-fid,varname,ndim,mdv,vardim,
148
     >                  varmin,varmax,stag,ierr)
149
            if (ierr.ne.0) goto 900
150
            call getcfn(-fid,cstfile,ierr)
151
            if (ierr.ne.0) goto 903
152
            call cdfopn(cstfile,cstid,ierr)
153
            if (ierr.ne.0) goto 903
154
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
155
            if (ierr.ne.0) goto 903
156
            call clscdf(cstid,ierr)
157
            if (ierr.ne.0) goto 903
3 michaesp 158
 
44 michaesp 159
         else
160
            varname = fieldname
161
            call getdef(-fid,varname,ndim,mdv,vardim,
162
     >                  varmin,varmax,stag,ierr)
163
            if (ierr.ne.0) goto 900
164
 
165
         endif
3 michaesp 166
 
44 michaesp 167
c        Set the grid dimensions and constants - vardim(3) is taken from constants file
168
         nx   = vardim(1)
169
         ny   = vardim(2)
170
         nz   = vardim(3)
171
         xmin = varmin(1)
172
         ymin = varmin(2)
173
         xmax = varmax(1)
174
         ymax = varmax(2)
175
         dx   = (xmax-xmin)/real(nx-1)
176
         dy   = (ymax-ymin)/real(ny-1)
3 michaesp 177
 
44 michaesp 178
c        Get pole position - if no constants file available, set default pole
179
         call getcfn(-fid,cstfile,ierr)
180
         if (ierr.eq.0) then  
181
            call cdfopn(cstfile,cstid,ierr)
182
            if (ierr.ne.0) goto 903
183
            call getpole(cstid,pollon,pollat,ierr)
184
            if (ierr.ne.0) goto 903
185
            call clscdf(cstid,ierr)
186
            if (ierr.ne.0) goto 903
187
         else
188
            pollon = 0.
189
            pollat = 90.
190
         endif
25 michaesp 191
 
44 michaesp 192
c     Get non-constant grid parameters (surface pressure and vertical grid)
193
      else
194
 
195
c        Special handling for fieldname 'P.ML' -> in this case the fields
196
c        P and PS are available on the data file and can be read in. There
197
c        is no need to reconstruct it from PS,AK and BK. This mode is
198
c        used for model-level (2D) trajectories
199
         if ( fieldname.eq.'P.ML' ) then
3 michaesp 200
 
44 michaesp 201
c           Get the right time to read
202
            call gettimes(fid,times,ntimes,ierr)
203
            if (ierr.ne.0) goto 901
204
            isok=0
205
            do i=1,ntimes
206
               if (abs(time-times(i)).lt.eps) then
207
                  isok = 1
208
                  rtime = times(i)
209
               elseif (timecheck.eq.'no') then
210
                  isok = 1
211
                  rtime = times(1)
212
               endif
213
            enddo
29 michaesp 214
 
44 michaesp 215
c           Read surface pressure and 3D pressure
216
            varname='PS'
217
            call getdat(fid,varname,rtime,0,ps,ierr)
218
            if (ierr.ne.0) goto 904
219
            varname='P'
220
            call getdat(fid,varname,rtime,0,p3,ierr)
221
            if (ierr.ne.0) goto 904
19 michaesp 222
 
44 michaesp 223
c           Set MDV to 1050. - otherwise interpolation routines don't work
224
            do i=1,nx
225
              do j=1,ny
226
                do k=1,nz
227
                   if ( p3(i,j,k).lt.0. ) p3(i,j,k) = 1050.
228
                enddo
229
              enddo
230
            enddo
19 michaesp 231
 
44 michaesp 232
c           Don't care about other stuff - finish subroutine
233
            goto 800
19 michaesp 234
 
29 michaesp 235
         endif
236
 
44 michaesp 237
c        Get full grid info - in particular staggering flag; set flag for 2D tracing
238
         if ( fieldname.eq.'PLEV' ) then
239
            varname = 'PS'
240
            stagz   = 0.
241
            call getdef(fid,varname,ndim,mdv,vardim,
242
     >                  varmin,varmax,stag,ierr)
243
            if (ierr.ne.0) goto 900
244
            call getcfn(fid,cstfile,ierr)
245
            if (ierr.ne.0) goto 903
246
            call cdfopn(cstfile,cstid,ierr)
247
            if (ierr.ne.0) goto 903
248
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
249
            if (ierr.ne.0) goto 903
250
            call clscdf(cstid,ierr)
251
            if (ierr.ne.0) goto 903
252
 
253
         elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then       
254
            varname = 'PS'
255
            stagz   = -0.5
256
            call getdef(fid,varname,ndim,mdv,vardim,
257
     >                  varmin,varmax,stag,ierr)
258
            if (ierr.ne.0) goto 900
259
            call getcfn(fid,cstfile,ierr)
260
            if (ierr.ne.0) goto 903
261
            call cdfopn(cstfile,cstid,ierr)
262
            if (ierr.ne.0) goto 903
263
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
264
            if (ierr.ne.0) goto 903
265
            call clscdf(cstid,ierr)
266
            if (ierr.ne.0) goto 903
25 michaesp 267
 
44 michaesp 268
         else
269
            varname=fieldname
270
            call getdef(fid,varname,ndim,mdv,vardim,
271
     >                  varmin,varmax,stag,ierr)
272
            if (ierr.ne.0) goto 900
273
            if (vardim(3).eq.1) is2d = 1
274
         endif
275
 
276
c        Get time information (check if time is correct)
277
         call gettimes(fid,times,ntimes,ierr)
278
         if (ierr.ne.0) goto 901
279
         isok=0
280
         do i=1,ntimes
281
            if (abs(time-times(i)).lt.eps) then
3 michaesp 282
               isok = 1
283
               rtime = times(i)
44 michaesp 284
            elseif (timecheck.eq.'no') then
3 michaesp 285
               isok = 1
286
               rtime = times(1)
44 michaesp 287
            endif
288
         enddo
289
         if ( isok.eq.0) goto 905
3 michaesp 290
 
44 michaesp 291
c        If 2D tracing requested: take dummay values for PS, AKLEV,BKLEV,AKLAY,BKLAY
292
         if ( is2d.eq.1 ) then
293
 
294
            do i=1,nx
295
               do j=1,ny
296
                  ps(i,j) = 1050.
297
               enddo
298
            enddo
299
 
300
            do k=1,nz
301
               aklev(k) = 0.
302
               bklev(k) = real(nz-k)/real(nz-1) 
303
               aklay(k) = 0.
304
               bklay(k) = real(nz-k)/real(nz-1) 
305
            enddo
21 michaesp 306
 
44 michaesp 307
c        3D tracing - read PS, AKLEV,BKLEV,AKLAY;BKLAY
308
         else
3 michaesp 309
 
44 michaesp 310
c           Read the level coefficients from the constants file
311
            call getcfn(fid,cstfile,ierr)
312
            if (ierr.ne.0) goto 903
313
            call cdfopn(cstfile,cstid,ierr)
314
            if (ierr.ne.0) goto 903
315
            call getlevs(cstid,vardim(3),aklev,bklev,aklay,bklay,ierr)
316
            if (ierr.ne.0) goto 903
317
            call clscdf(cstid,ierr)
318
            if (ierr.ne.0) goto 903
27 michaesp 319
 
44 michaesp 320
c           Check whether PS is needed to get the 3d pressure field
321
            plev = 1
322
            do i=1,nz
323
              if ( (abs(stagz).lt.eps).and.(abs(bklev(i)).gt.eps) ) then
324
                plev = 0
325
              endif
326
              if ( (abs(stagz).gt.eps).and.(abs(bklay(i)).gt.eps) ) then
327
                plev = 0
328
              endif
329
            enddo
3 michaesp 330
 
44 michaesp 331
c           Read surface pressure if needed
332
            if ( plev.ne.1 ) then
333
              varname='PS'
334
              call getdat(fid,varname,rtime,0,ps,ierr)
335
              if (ierr.ne.0) goto 904
336
            else
337
              do i=1,nx
338
                do j=1,ny
339
                    ps(i,j) = 1000.
340
                enddo
341
              enddo
342
            endif
3 michaesp 343
 
29 michaesp 344
         endif
25 michaesp 345
 
44 michaesp 346
c        Calculate layer and level pressures
347
         do i=1,nx
348
            do j=1,ny
349
               do k=1,nz
350
                  if ( abs(stagz).lt.eps ) then
351
                     p3(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
352
                  else
353
                     p3(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
354
                  endif
355
               enddo
356
            enddo
3 michaesp 357
         enddo
358
 
44 michaesp 359
c        Set the ak and bk for the vertical grid
360
         do k=1,nz
361
            if ( abs(stagz).lt.eps ) then
362
               ak(k)=aklev(k)
363
               bk(k)=bklev(k)
364
            else
365
               ak(k)=aklay(k)
366
               bk(k)=bklay(k)
367
            endif
3 michaesp 368
         enddo
23 michaesp 369
 
21 michaesp 370
      endif
44 michaesp 371
 
372
c     Exit point for subroutine
373
 800  continue
3 michaesp 374
      return
375
 
44 michaesp 376
c     Exception handling
377
 900  print*,'Cannot retrieve grid dimension from ',fid
378
      stop
379
 901  print*,'Cannot retrieve grid parameters from ',fid
380
      stop
381
 902  print*,'Grid inconsistency detected for ',fid
382
      stop
383
 903  print*,'Problem with level coefficients from ',trim(cstfile)
384
      stop
385
 904  print*,'Cannot read surface pressure from ',trim(cstfile)
386
      stop
387
 905  print*,'Cannot find time ',time,' on ',fid
388
      stop
389
 906  print*,'Unable to get grid info ',fid
390
      stop
391
 
3 michaesp 392
      end
393
 
394
c     ------------------------------------------------------------
395
c     Read wind information
396
c     ------------------------------------------------------------
397
 
398
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
399
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
400
     >                       timecheck)
401
 
402
c     Read the wind component <fieldname> from the file with identifier
403
c     <fid> and save it in the 3d array <field>. The vertical staggering 
404
c     information is provided in <stagz> and gives the reference to either
405
c     the layer or level field from <input_grid>. A consistency check is
406
c     performed to have an agreement with the grid specified by <xmin,xmax,
407
c     ymin,ymax,dx,dy,nx,ny,nz>.
408
 
409
      implicit none
410
 
411
c     Declaration of variables and parameters
412
      integer      fid                 ! File identifier
413
      character*80 fieldname           ! Name of the wind field
414
      integer      nx,ny,nz            ! Dimension of fields
415
      real         field(nx,ny,nz)     ! 3d wind field
416
      real         stagz               ! Staggering in the z direction
417
      real         mdv                 ! Missing data flag
418
      real         xmin,xmax,ymin,ymax ! Domain size
419
      real         dx,dy               ! Horizontal resolution
420
      real         time                ! Time
421
      character*80 timecheck           ! Either 'yes' or 'no'
422
 
423
c     Numerical and physical parameters
424
      real        eps                 ! Numerical epsilon
425
      parameter  (eps=0.001)
426
      real        notimecheck         ! 'Flag' for no time check
427
      parameter  (notimecheck=7.26537)
428
 
429
c     Netcdf variables
430
      integer      ierr
431
      character*80 varname
432
      integer      vardim(4)
433
      real         varmin(4),varmax(4)
434
      real         stag(4)
435
      integer      ndim
436
      real         times(10)
437
      integer      ntimes
438
      character*80 cstfile
439
      integer      cstid
440
      real         aklay(200),bklay(200),aklev(200),bklev(200)
441
      real         ps(nx,ny)
442
 
443
c     Auxiliary variables
444
      integer      isok
445
      integer      i,j,k
446
      integer      nz1
447
      real         rtime
448
 
44 michaesp 449
c     Read variable definition - for P, PLEV and PLAY: load also ak,bk
450
      if ( ( fieldname.eq.'PLEV' ).or.
451
     >     ( fieldname.eq.'PLAY' ).or.
452
     >     ( fieldname.eq.'P'    ) )
453
     >then
454
         call getcfn(fid,cstfile,ierr)
455
         if (ierr.ne.0) goto 905
456
         call cdfopn(cstfile,cstid,ierr)
457
         if (ierr.ne.0) goto 905
458
         call getlevs(cstid,nz1,aklev,bklev,aklay,bklay,ierr)
459
         if (ierr.ne.0) goto 905
460
         call clscdf(cstid,ierr)
461
         if (ierr.ne.0) goto 905
462
         varname = 'PS'
463
         call getdef(fid,varname,ndim,mdv,vardim,
464
     >               varmin,varmax,stag,ierr)
465
         vardim(3) = nz1
466
         if (ierr.ne.0) goto 900
467
 
36 michaesp 468
      else
44 michaesp 469
         varname = fieldname
470
         call getdef(fid,varname,ndim,mdv,vardim,
471
     >               varmin,varmax,stag,ierr)
472
         if (ierr.ne.0) goto 900
473
         stagz=stag(3)
36 michaesp 474
      endif
29 michaesp 475
 
44 michaesp 476
c     Get time information (set time to first one in the file)
477
      call gettimes(fid,times,ntimes,ierr)
478
      if (ierr.ne.0) goto 902
479
      isok=0
480
      do i=1,ntimes
481
         if (abs(time-times(i)).lt.eps) then
482
            isok = 1
483
            rtime = times(i)
484
         elseif (timecheck.eq.'no') then
485
            isok = 1
486
            rtime = times(1)
487
         endif
3 michaesp 488
      enddo
44 michaesp 489
      if ( isok.eq.0 ) goto 904
3 michaesp 490
 
44 michaesp 491
c     Read  field
492
      if ( ( fieldname.eq.'P' ).or.(fieldname.eq.'PLAY') )  then       ! P, PLAY
493
         stagz   = -0.5
494
         varname = 'PS'
495
         call getdat(fid,varname,rtime,0,ps,ierr)
496
         if (ierr.ne.0) goto 903
497
         do i=1,nx
498
            do j=1,ny
499
               do k=1,nz
500
                  field(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
501
               enddo
502
            enddo
25 michaesp 503
         enddo
44 michaesp 504
 
505
      elseif ( fieldname.eq.'PLEV' )  then                             ! PLEV
506
         stagz   = 0.
507
         varname = 'PS'
508
         call getdat(fid,varname,rtime,0,ps,ierr)
509
         if (ierr.ne.0) goto 903
510
         do i=1,nx
511
            do j=1,ny
512
               do k=1,nz
513
                  field(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
514
               enddo
29 michaesp 515
            enddo
44 michaesp 516
         enddo
29 michaesp 517
 
44 michaesp 518
      else                                                             ! Other fields
519
         varname=fieldname
520
         call getdat(fid,varname,rtime,0,field,ierr)
521
         if (ierr.ne.0) goto 903
522
 
25 michaesp 523
      endif
524
 
3 michaesp 525
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
526
      if ( vardim(3).eq.1 ) then
44 michaesp 527
         do i=1,nx
528
            do j=1,ny
529
               do k=1,nz
530
                  field(i,j,k) = field(i,j,1)
3 michaesp 531
               enddo
532
            enddo
533
         enddo
534
      endif
44 michaesp 535
 
3 michaesp 536
 
537
 
44 michaesp 538
c     Exception handling
3 michaesp 539
      return
44 michaesp 540
 
541
 900  print*,'Cannot retrieve definition for ',trim(varname),'  ',fid
542
      stop
543
 901  print*,'Grid inconsistency detected for ',trim(varname),'  ',fid
544
      stop
545
 902  print*,'Cannot retrieve time for ',trim(varname),'  ',fid
546
      stop
547
 903  print*,'Cannot load wind component ',trim(varname),'  ',fid
548
      stop
549
 904  print*,'Cannot load time ',time,' for ',trim(varname),'  ',fid
550
      stop
551
 905  print*,'Cannot load time vertical grid AK, BK from file  ',fid
552
      stop
553
 
3 michaesp 554
      end
555
 
556
c     ------------------------------------------------------------
557
c     Close input file
558
c     ------------------------------------------------------------
559
 
560
      subroutine input_close(fid)
561
 
562
c     Close the input file with file identifier <fid>.
563
 
564
      implicit none
565
 
566
c     Declaration of subroutine parameters
567
      integer fid
568
 
569
c     Auxiliary variables
570
      integer ierr
571
 
572
c     Close file
44 michaesp 573
      call clscdf(fid,ierr)
574
 
21 michaesp 575
      end
576
 
577
c     ------------------------------------------------------------
578
c     Get a list of variables on netCDF file
579
c     ------------------------------------------------------------
580
 
581
      subroutine input_getvars(fid,vnam,nvars)
582
 
583
c     List of variables on netCDF file
584
 
585
      implicit none
586
 
587
c     Declaration of subroutine parameters
588
      integer      fid
589
      integer      nvars
590
      character*80 vnam(200)
591
 
592
c     Auxiliary variables
593
      integer ierr
594
 
44 michaesp 595
c     Get list and save      
596
      call getvars(fid,nvars,vnam,ierr)
3 michaesp 597
 
598
      end