Subversion Repositories lagranto.ecmwf

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
5 michaesp 1
      subroutine wricst(cstnam,datar,aklev,bklev,aklay,bklay,stdate)
2
C------------------------------------------------------------------------
3
 
4
C     Creates the constants file for NetCDF files containing ECMWF
5
C     data. The constants file is compatible with the one created
6
C     for EM data (with subroutine writecst).
7
C
8
C     Input parameters:
9
C
10
C     cstnam    name of constants file
11
C     datar     array contains all required parameters to write file
12
C               datar(1):       number of points along x        
13
C               datar(2):       number of points along y
14
C               datar(3):       maximum latitude of data region (ymax)
15
C               datar(4):       minimum longitude of data region (xmin)
16
C               datar(5):       minimum latitude of data region (ymin)
17
C               datar(6):       maximum longitude of data region (xmax)
18
C               datar(7):       grid increment along x
19
C               datar(8):       grid increment along y
20
C               datar(9):       number of levels        
21
C		datar(10):	data type (forecast or analysis)
22
C		datar(11):	data version
23
C		datar(12):	constants file version
24
C		datar(13):	longitude of pole of coordinate system
25
C		datar(14):	latitude of pole of coordinate system
26
C     aklev     array contains the aklev values
27
C     bklev	array contains the bklev values
28
C     aklay     array contains the aklay values
29
C     bklay     array contains the bklay values
30
C     stdate    array contains date (year,month,day,time,step) of first
31
C               field on file (start-date), dimensionised as stdate(5)
32
C------------------------------------------------------------------------
33
 
34
 
35
      include "netcdf.inc"
36
 
37
      integer   nchar,maxlev
38
 
39
      parameter (nchar=20,maxlev=32)
40
      real	aklev(maxlev),bklev(maxlev)
41
      real      aklay(maxlev),bklay(maxlev)
42
      real	pollat,latmin,latmax
43
      integer   datar(14)
44
      integer	stdate(5)
45
      character*80 cstnam
46
 
47
C     declarations for constants-variables
48
 
49
      integer   nz
50
      integer   dattyp, datver, cstver
51
 
52
C     further declarations
53
 
54
      integer	ierr			! error flag
55
      integer	cdfid			! NetCDF id
56
      integer	xid,yid,zid		! dimension ids
57
      integer	pollonid, pollatid,	! variable ids
58
     >		aklevid, bklevid, aklayid, bklayid,
59
     >		lonminid, lonmaxid, latminid, latmaxid,
60
     >		dellonid, dellatid,
61
     >		startyid, startmid, startdid, starthid, startsid,
62
     >		dattypid, datverid, cstverid
63
 
64
      nz=datar(9)			! number of levels
65
 
66
C     Set data-type and -version, version of cst-file-format
67
 
68
      dattyp=datar(10)
69
      datver=datar(11)
70
      cstver=datar(12)
71
 
72
C     Initially set error to false
73
 
74
      ierr=0
75
 
76
C     Create constants file
77
 
78
      cdfid=nccre(trim(cstnam),NCCLOB,ierr)
79
 
80
C     Define the dimensions
81
 
82
      xid = ncddef (cdfid,'nx',datar(1),ierr)
83
      yid = ncddef (cdfid,'ny',datar(2),ierr)
84
      zid = ncddef (cdfid,'nz',datar(9),ierr)
85
 
86
C     Define integer constants
87
 
88
      pollonid = ncvdef(cdfid,'pollon', NCFLOAT,0,0,ierr)
89
      pollatid = ncvdef(cdfid,'pollat', NCFLOAT,0,0,ierr)
90
 
91
      aklevid = ncvdef (cdfid, 'aklev', NCFLOAT, 1, zid, ierr)
92
      bklevid = ncvdef (cdfid, 'bklev', NCFLOAT, 1, zid, ierr)
93
      aklayid = ncvdef (cdfid, 'aklay', NCFLOAT, 1, zid, ierr)
94
      bklayid = ncvdef (cdfid, 'bklay', NCFLOAT, 1, zid, ierr)
95
 
96
      lonminid = ncvdef (cdfid, 'lonmin', NCFLOAT, 0, 0, ierr)
97
      lonmaxid = ncvdef (cdfid, 'lonmax', NCFLOAT, 0, 0, ierr)
98
      latminid = ncvdef (cdfid, 'latmin', NCFLOAT, 0, 0, ierr)
99
      latmaxid = ncvdef (cdfid, 'latmax', NCFLOAT, 0, 0, ierr)
100
      dellonid = ncvdef (cdfid, 'dellon', NCFLOAT, 0, 0, ierr)
101
      dellatid = ncvdef (cdfid, 'dellat', NCFLOAT, 0, 0, ierr)
102
      startyid = ncvdef (cdfid, 'starty', NCLONG, 0, 0, ierr)
103
      startmid = ncvdef (cdfid, 'startm', NCLONG, 0, 0, ierr)
104
      startdid = ncvdef (cdfid, 'startd', NCLONG, 0, 0, ierr)
105
      starthid = ncvdef (cdfid, 'starth', NCLONG, 0, 0, ierr)
106
      startsid = ncvdef (cdfid, 'starts', NCLONG, 0, 0, ierr)
107
      dattypid = ncvdef (cdfid, 'dattyp', NCLONG, 0, 0, ierr)
108
      datverid = ncvdef (cdfid, 'datver', NCLONG, 0, 0, ierr)
109
      cstverid = ncvdef (cdfid, 'cstver', NCLONG, 0, 0, ierr)
110
 
111
C     Leave define mode
112
 
113
      call ncendf(cdfid,ierr)
114
 
115
C     Store levels
116
      call ncvpt(cdfid, aklevid, 1, nz, aklev, ierr)
117
      call ncvpt(cdfid, bklevid, 1, nz, bklev, ierr)
118
      call ncvpt(cdfid, aklayid, 1, nz, aklay, ierr)
119
      call ncvpt(cdfid, bklayid, 1, nz, bklay, ierr)
120
 
121
C     Store position of pole (trivial for ECMWF data)
122
      call ncvpt1(cdfid, pollonid, 1, real(datar(13))/1000., ierr)
123
      if (datar(14).gt.0) then
124
        pollat=min(real(datar(14))/1000.,90.)
125
      else
126
        pollat=max(real(datar(14))/1000.,-90.)
127
      endif
128
      call ncvpt1(cdfid, pollatid, 1, pollat, ierr)
129
 
130
C     Store horizontal data borders and grid increments
131
      call ncvpt1(cdfid, lonminid, 1, real(datar(4))/1000., ierr)
132
      call ncvpt1(cdfid, lonmaxid, 1, real(datar(6))/1000., ierr)
133
      latmin=max(real(datar(5))/1000.,-90.)
134
      latmax=min(real(datar(3))/1000.,90.)
135
      call ncvpt1(cdfid, latminid, 1, latmin, ierr)
136
      call ncvpt1(cdfid, latmaxid, 1, latmax, ierr)
137
      call ncvpt1(cdfid, dellonid, 1, real(datar(7))/1000., ierr)
138
      call ncvpt1(cdfid, dellatid, 1, real(datar(8))/1000., ierr)
139
 
140
C     Store date of first field on file (start-date)
141
      call ncvpt1(cdfid, startyid, 1, stdate(1), ierr)
142
      call ncvpt1(cdfid, startmid, 1, stdate(2), ierr)
143
      call ncvpt1(cdfid, startdid, 1, stdate(3), ierr)
144
      call ncvpt1(cdfid, starthid, 1, stdate(4), ierr)
145
      call ncvpt1(cdfid, startsid, 1, stdate(5), ierr)
146
 
147
C     Store datatype and version
148
      call ncvpt1(cdfid, dattypid, 1, dattyp, ierr)
149
      call ncvpt1(cdfid, datverid, 1, datver, ierr)
150
 
151
C     Store version of the constants file format
152
      call ncvpt1(cdfid, cstverid, 1, cstver, ierr)
153
 
154
C     Store strings
155
 
156
      call ncclos(cdfid,ierr)
157
      return
158
 
159
      end
160
      subroutine writelmcst(cdfid,nx,ny,nz,pollon,pollat,lonmin,
161
     &lonmax,latmin,latmax,dellon,dellat,dattyp,datver,cstver,
162
     &psref,tstar,tbeta,pintf,p0top,idate)
163
c     ------------------------------------------------------------------
164
 
165
      implicit none
166
 
167
      integer   cdfid
168
 
169
c     deklarationen der constants-variablen
170
      real       pollon,pollat
171
      real       lonmin,lonmax,latmin,latmax,dellon,dellat
172
      integer    idate(5)
173
      integer    nx,ny,nz
174
      integer    dattyp, datver, cstver
175
      real       psref, tstar, tbeta, pintf, p0top
176
 
177
      include 'netcdf.inc'
178
 
179
* netcdf declaration
180
      integer   iret, k
181
* dimension ids
182
      integer  nxdim, nydim, nzdim
183
* variable ids
184
      integer  startyid, startmid, startdid, starthid
185
* variable shapes, corners and edge lengths
186
      integer dims(1), corner(1), edges(1)
187
 
188
* enter define mode
189
      call ncredf(cdfid, iret)
190
 
191
      startyid = ncvdef (cdfid, 'starty', NCLONG, 0, 0, iret)
192
      startmid = ncvdef (cdfid, 'startm', NCLONG, 0, 0, iret)
193
      startdid = ncvdef (cdfid, 'startd', NCLONG, 0, 0, iret)
194
      starthid = ncvdef (cdfid, 'starth', NCLONG, 0, 0, iret)
195
 
196
* store the rest as global attributes
197
* store nx,ny,nz
198
      call ncapt(cdfid,NCGLOBAL,'nx',NCLONG,1,nx,iret)
199
      call ncapt(cdfid,NCGLOBAL,'ny',NCLONG,1,ny,iret)
200
      call ncapt(cdfid,NCGLOBAL,'nz',NCLONG,1,nz,iret)
201
 
202
* store pollon, pollat
203
      call ncapt(cdfid,NCGLOBAL,'pollon',NCFLOAT,1,pollon,iret)
204
      call ncapt(cdfid,NCGLOBAL,'pollat',NCFLOAT,1,pollat,iret)
205
 
206
* store lonmin, etc
207
      call ncapt(cdfid,NCGLOBAL,'lonmin',NCFLOAT,1,lonmin,iret)
208
      call ncapt(cdfid,NCGLOBAL,'lonmax',NCFLOAT,1,lonmax,iret)
209
      call ncapt(cdfid,NCGLOBAL,'latmin',NCFLOAT,1,latmin,iret)
210
      call ncapt(cdfid,NCGLOBAL,'latmax',NCFLOAT,1,latmax,iret)
211
      call ncapt(cdfid,NCGLOBAL,'dellon',NCFLOAT,1,dellon,iret)
212
      call ncapt(cdfid,NCGLOBAL,'dellat',NCFLOAT,1,dellat,iret)
213
 
214
* store data type and version
215
      call ncapt(cdfid,NCGLOBAL,'dattyp',NCLONG,1,dattyp,iret)
216
      call ncapt(cdfid,NCGLOBAL,'datver',NCLONG,1,datver,iret)
217
      call ncapt(cdfid,NCGLOBAL,'cstver',NCLONG,1,cstver,iret)
218
 
219
* store information of lm model vertical grid
220
      call ncapt(cdfid,NCGLOBAL,'psref',NCFLOAT,1,psref,iret)
221
      call ncapt(cdfid,NCGLOBAL,'tstar',NCFLOAT,1,tstar,iret)
222
      call ncapt(cdfid,NCGLOBAL,'tbeta',NCFLOAT,1,tbeta,iret)
223
      call ncapt(cdfid,NCGLOBAL,'pintf',NCFLOAT,1,pintf,iret)
224
      call ncapt(cdfid,NCGLOBAL,'p0top',NCFLOAT,1,p0top,iret)
225
 
226
* leave define mode
227
      call ncendf(cdfid, iret)
228
 
229
* store starty, etc
230
      corner(1) = 1
231
      edges(1) = 1
232
      call ncvpt(cdfid, startyid, corner, edges, idate(1), iret)
233
      call ncvpt(cdfid, startmid, corner, edges, idate(2), iret)
234
      call ncvpt(cdfid, startdid, corner, edges, idate(3), iret)
235
      call ncvpt(cdfid, starthid, corner, edges, idate(4), iret)
236
 
237
      end
238
      subroutine globcst(cdfnam,datar,aklev,bklev,aklay,bklay,stdate)
239
C------------------------------------------------------------------------
240
C+
241
C NAME:
242
C     subroutine globcst
243
C
244
C PURPOSE:
245
C     instead of writing a constants-file (*_cst), the information
246
C     is added to the netCDF file as global variables
247
C     the data format is compatible with the one requested by
248
C     the IVE ETH/MIT version, contact author about details
249
C
250
C CATEGORY:
251
C     model,netCDF
252
C
253
C CALLING SEQUENCE:
254
C     subroutine globcst(cdfnam,datar,aklev,bklev,aklay,bklay,stdate)
255
C
256
C INPUTS:
257
C     cdfnam    name of netCDF file
258
C               The file needs to exist, otherwise an ERROR occurs,
259
C               i.e. nothing is done
260
C     datar     array contains all required parameters to write file
261
C               datar(1):       number of points along x
262
C               datar(2):       number of points along y
263
C               datar(3):       maximum latitude of data region (ymax)
264
C               datar(4):       minimum longitude of data region (xmin)
265
C               datar(5):       minimum latitude of data region (ymin)
266
C               datar(6):       maximum longitude of data region (xmax)
267
C               datar(7):       grid increment along x
268
C               datar(8):       grid increment along y
269
C               datar(9):       number of levels
270
C               datar(10):      data type (forecast or analysis)
271
C               datar(11):      data version
272
C               datar(12):      constants file version
273
C               datar(13):      longitude of pole of coordinate system
274
C               datar(14):      latitude of pole of coordinate system
275
C     aklev     array contains the aklev values
276
C     bklev     array contains the bklev values
277
C     aklay     array contains the aklay values
278
C     bklay     array contains the bklay values
279
C     stdate    array contains date (year,month,day,time,step) of first
280
C               field on file (start-date), dimensionised as stdate(5)
281
C     list    the griblist-ASCII-file
282
C     varno   the GRIB code number
283
C
284
C OUTPUTS:
285
C     Adds cdf-information to EXISTING netCDF-file
286
C
287
C MODIFICATION HISTORY:
288
C
289
C     June  93    Christoph Schaer (ETHZ) created
290
C     Nov   93    Heini Wernli (ETHZ) wricst
291
C     Nov   98    David N. Bresch (MIT) wricst to globcst
292
C-
293
 
294
C     Sun include statement.
295
      include "netcdf.inc"
296
 
297
      integer   nchar,maxlev
298
 
299
      parameter (nchar=20,maxlev=32)
300
      real      aklev(maxlev),bklev(maxlev)
301
      real      aklay(maxlev),bklay(maxlev)
302
      integer   datar(14)
303
      integer   stdate(5)
304
      character*80 cdfnam
305
 
306
C     declarations for constants-variables
307
 
308
      integer   nz
309
      integer   dattyp, datver, cstver
310
 
311
C     further declarations
312
 
313
      integer   ierr                    ! error flag
314
      integer   cdfid                   ! NetCDF id
315
      integer   xid,yid,zid             ! dimension ids
316
      integer   pollonid, pollatid,     ! variable ids
317
     >          aklevid, bklevid, aklayid, bklayid,
318
     >          lonminid, lonmaxid, latminid, latmaxid,
319
     >          dellonid, dellatid,
320
     >          startyid, startmid, startdid, starthid, startsid,
321
     >          dattypid, datverid, cstverid
322
 
323
      nz=datar(9)                       ! number of levels
324
 
325
C     Set data-type and -version, version of cst-file-format
326
 
327
      dattyp=datar(10)
328
      datver=datar(11)
329
      cstver=datar(12)
330
 
331
C     Initially set error to false
332
 
333
      ierr=0
334
 
335
C     open the netCDF-file:
336
 
337
      call cdfwopn(cdfnam,cdfid,ierr)
338
      if (ierr.ne.0) then
339
         print*,'ERROR opening netCDF-file ',cdfnam
340
         return
341
      endif
342
 
343
C     Put file into define mode
344
      call ncredf(cdfid,ierr)
345
      if (ierr.ne.0) then
346
         print*,'ERROR switching to netCDF redefine mode'
347
         return
348
      endif
349
 
350
C     Define the dimensions
351
 
352
      xid = ncddef (cdfid,'nx',datar(1),ierr)
353
      yid = ncddef (cdfid,'ny',datar(2),ierr)
354
      zid = ncddef (cdfid,'nz',datar(9),ierr)
355
 
356
C     Define integer constants
357
 
358
      pollonid = ncvdef(cdfid,'pollon', NCFLOAT,0,0,ierr)
359
      pollatid = ncvdef(cdfid,'pollat', NCFLOAT,0,0,ierr)
360
 
361
      aklevid = ncvdef (cdfid, 'aklev', NCFLOAT, 1, zid, ierr)
362
      bklevid = ncvdef (cdfid, 'bklev', NCFLOAT, 1, zid, ierr)
363
      aklayid = ncvdef (cdfid, 'aklay', NCFLOAT, 1, zid, ierr)
364
      bklayid = ncvdef (cdfid, 'bklay', NCFLOAT, 1, zid, ierr)
365
 
366
      lonminid = ncvdef (cdfid, 'lonmin', NCFLOAT, 0, 0, ierr)
367
      lonmaxid = ncvdef (cdfid, 'lonmax', NCFLOAT, 0, 0, ierr)
368
      latminid = ncvdef (cdfid, 'latmin', NCFLOAT, 0, 0, ierr)
369
      latmaxid = ncvdef (cdfid, 'latmax', NCFLOAT, 0, 0, ierr)
370
      dellonid = ncvdef (cdfid, 'dellon', NCFLOAT, 0, 0, ierr)
371
      dellatid = ncvdef (cdfid, 'dellat', NCFLOAT, 0, 0, ierr)
372
      startyid = ncvdef (cdfid, 'starty', NCLONG, 0, 0, ierr)
373
      startmid = ncvdef (cdfid, 'startm', NCLONG, 0, 0, ierr)
374
      startdid = ncvdef (cdfid, 'startd', NCLONG, 0, 0, ierr)
375
      starthid = ncvdef (cdfid, 'starth', NCLONG, 0, 0, ierr)
376
      startsid = ncvdef (cdfid, 'starts', NCLONG, 0, 0, ierr)
377
      dattypid = ncvdef (cdfid, 'dattyp', NCLONG, 0, 0, ierr)
378
      datverid = ncvdef (cdfid, 'datver', NCLONG, 0, 0, ierr)
379
      cstverid = ncvdef (cdfid, 'cstver', NCLONG, 0, 0, ierr)
380
 
381
C     Leave define mode
382
 
383
      call ncendf(cdfid,ierr)
384
      if (ierr.ne.0) then
385
         print*,'ERROR exiting define mode'
386
         return
387
      endif
388
 
389
C     Store levels
390
      call ncvpt(cdfid, aklevid, 1, nz, aklev, ierr)
391
      call ncvpt(cdfid, bklevid, 1, nz, bklev, ierr)
392
      call ncvpt(cdfid, aklayid, 1, nz, aklay, ierr)
393
      call ncvpt(cdfid, bklayid, 1, nz, bklay, ierr)
394
 
395
C     Store position of pole (trivial for ECMWF data)
396
      call ncvpt1(cdfid, pollonid, 1, real(datar(13))/1000., ierr)
397
      call ncvpt1(cdfid, pollatid, 1, real(datar(14))/1000., ierr)
398
 
399
C     Store horizontal data borders and grid increments
400
      call ncvpt1(cdfid, lonminid, 1, real(datar(4))/1000., ierr)
401
      call ncvpt1(cdfid, lonmaxid, 1, real(datar(6))/1000., ierr)
402
      call ncvpt1(cdfid, latminid, 1, real(datar(5))/1000., ierr)
403
      call ncvpt1(cdfid, latmaxid, 1, real(datar(3))/1000., ierr)
404
      call ncvpt1(cdfid, dellonid, 1, real(datar(7))/1000., ierr)
405
      call ncvpt1(cdfid, dellatid, 1, real(datar(8))/1000., ierr)
406
 
407
C     Store date of first field on file (start-date)
408
      call ncvpt1(cdfid, startyid, 1, stdate(1), ierr)
409
      call ncvpt1(cdfid, startmid, 1, stdate(2), ierr)
410
      call ncvpt1(cdfid, startdid, 1, stdate(3), ierr)
411
      call ncvpt1(cdfid, starthid, 1, stdate(4), ierr)
412
      call ncvpt1(cdfid, startsid, 1, stdate(5), ierr)
413
 
414
C     Store datatype and version
415
      call ncvpt1(cdfid, dattypid, 1, dattyp, ierr)
416
      call ncvpt1(cdfid, datverid, 1, datver, ierr)
417
 
418
C     Store version of the constants file format
419
      call ncvpt1(cdfid, cstverid, 1, cstver, ierr)
420
 
421
      if (ierr.ne.0) then
422
         print*,'ERROR adding cst-date as global variables'
423
         return
424
      endif
425
 
426
C     Store strings
427
 
428
      call ncclos(cdfid,ierr)
429
      if (ierr.ne.0) then
430
         print*,'ERROR closing netCDF file'
431
      endif
432
 
433
      return
434
      end
435
      subroutine getsdat(cdfid,varnam,time,ix,iy,iz,sx,sy,sz,dat,error)
436
c-----------------------------------------------------------------------
437
c     Purpose:
438
c        This routine is called to read the data within a selected
439
c	 domain of a variable from an IVE-NetCDF file.
440
c        Prior to calling this routine, the file must be opened with
441
c        a call to opncdf (for extension) or crecdf (for creation) or
442
c        readcdf (for readonly).
443
c     Arguments:
444
c        cdfid   int   input   file-identifier
445
c                              (must be obtained by calling routine
446
c                              opncdf,readcdf  or crecdf)
447
c        varnam  char  input   the user-supplied variable name
448
c        time    real  input   the user-supplied time-level of the
449
c                              data to be read from the file (the time-
450
c                              levels stored in the file can be obtained
451
c                              with a call to gettimes).
452
c        ix/y/z  int   input   indices of lower left corner of selected
453
c			       data volume.
454
c	 sx/y/z  int   input   size of selected data volume
455
c        dat     real  output  data-array with dimensions (sx,sy,sz).
456
c        error   int   output  indicates possible errors found in this
457
c                              routine.
458
c                              error = 0   no errors detected.
459
c                              error = 1   the variable is not present on
460
c                                          the file.
461
c                              error = 2   the value of 'time' is not
462
c                                          known.to the file.
463
c			       error = 6,7,8   data volume too large
464
c                              error =10   another error.
465
c     History:
466
c       June  93    Christoph Schaer (ETHZ)  Created getdat
467
c	Nov   93    Heini Wernli (ETHZ)	     Created getsdat
468
c-----------------------------------------------------------------------
469
 
470
      include "netcdf.inc"
471
 
472
C     Declaration of local variables
473
      character*(*) varnam
474
      character*(20) chars
475
      integer cdfid
476
 
477
      integer     ix,iy,iz,sx,sy,sz
478
      real        dat(sx,sy,sz)
479
      real        misdat,varmin(3),varmax(3),stag(3)
480
      real        time, timeval
481
 
482
      integer     corner(4),edgeln(4),didtim,vardim(4),ndims
483
      integer     error, ierr
484
      integer     ntime
485
      integer     idtime,idvar,iflag
486
      integer     i
487
 
488
      call ncpopt(NCVERBOS)
489
 
490
c     access the variable
491
      call getdef (cdfid, trim(varnam), ndims, misdat,
492
     &                           vardim, varmin, varmax, stag, ierr)
493
      if (ierr.ne.0) then
494
        print *,'*ERROR* in getdef in getdat'
495
        error=1
496
        return
497
      endif
498
      idvar=ncvid(cdfid,trim(varnam),ierr)
499
      if (ierr.ne.0) then
500
        print *,'*ERROR* in ncvid in getsdat'
501
        error=1
502
        return
503
      endif
504
 
505
C     Get times-array
506
      didtim=ncdid(cdfid,'time',ierr)
507
      if (ierr.ne.0) then
508
        print *,'*ERROR* didtim in getsdat'
509
        error=10
510
        return
511
      endif
512
      call ncdinq(cdfid,didtim,chars,ntime,ierr)
513
      if (ierr.ne.0) then
514
        print *,'*ERROR* in ncdinq in getsdat'
515
        error=10
516
        return
517
      endif
518
      idtime=ncvid(cdfid,'time',ierr)
519
      if (ierr.ne.0) then
520
        print *,'*ERROR* in ncvid for time in getsdat'
521
        error=10
522
        return
523
      endif
524
c     find appropriate time-index
525
      iflag=0
526
      do i=1,ntime
527
        call ncvgt1(cdfid,idtime,i,timeval,ierr)
528
        if (ierr.ne.0) print *,'*ERROR* in ncvgt1 in getsdat'
529
        if (time.eq.timeval) iflag=i
530
      enddo
531
      if (iflag.eq.0) then
532
        error=2
533
        print *,'Error: Unknown time in getsdat'
534
        print *,time,timeval
535
        return
536
      endif
537
 
538
C     Define data volume to be written (index space)
539
      corner(1)=ix
540
      corner(2)=iy
541
      corner(3)=iz
542
      corner(4)=iflag
543
      edgeln(1)=sx
544
      edgeln(2)=sy
545
      edgeln(3)=sz
546
      edgeln(4)=1
547
 
548
C     Check if data volume is within data domain
549
 
550
      if (ix+sx-1.gt.vardim(1)) then
551
        error=7
552
        print *,'Error: data volume too large in x-direction'
553
        print *,ix,sx,vardim(1)
554
        return
555
      endif
556
      if (iy+sy-1.gt.vardim(2)) then
557
        error=8
558
        print *,'Error: data volume too large in y-direction'
559
        return
560
      endif
561
      if (iz+sz-1.gt.vardim(3)) then
562
        error=9
563
        print *,'Error: data volume too large in z-direction'
564
        return
565
      endif
566
 
567
C     Read data from NetCDF file
568
 
569
      call ncvgt(cdfid,idvar,corner,edgeln,dat,error)
570
      if (error.ne.0) then
571
        print *, 'corner ',corner(1),corner(2),corner(3)
572
        print *, 'edgeln ',edgeln(1),edgeln(2),edgeln(3)
573
        print *, '*ERROR* in ncvgt in getsdat'
574
        error=10
575
      endif
576
      end
577
      subroutine getlevs(cstid,nlev,aklev,bklev,aklay,bklay,error)
578
c-----------------------------------------------------------------------
579
c     Purpose:
580
c     	This routine is called to get the level arrays aklev and
581
c	bklev from a NetCDF constants file.
582
c     Arguments:
583
c	cstid     int	input   identifier for NetCDF constants file
584
c	nlev	  int	input	number of levels
585
c	aklev     real	output  array contains all aklev values
586
c       bklev     real  output  array contains all bklev values
587
c	aklay	  real  output	array contains all aklay values
588
c	bklay	  real	output	array contains all bklay values
589
c	error	  int	output	error flag
590
c				error = 0   no errors detected
591
c				error = 1   error detected
592
c     History:
593
c	Aug. 93	  Heini Wernli		Created.
594
c-----------------------------------------------------------------------
595
 
596
      integer   error
597
 
598
      integer   cstid
599
      integer   ncdid,ncvid		! NetCDF functions
600
      integer   didz,idak,idbk,idaky,idbky
601
      integer   nlev
602
      real      aklev(nlev),bklev(nlev),aklay(nlev),bklay(nlev)
603
      character*(20) dimnam
604
      integer   i
605
 
606
      didz	=ncdid(cstid,'nz',error)
607
      if (error.ne.0) goto 920
608
      idak	=ncvid(cstid,'aklev',error)
609
      if (error.ne.0) goto 920
610
      idbk	=ncvid(cstid,'bklev',error)
611
      if (error.ne.0) goto 920
612
      idaky     =ncvid(cstid,'aklay',error)
613
      if (error.ne.0) goto 920
614
      idbky     =ncvid(cstid,'bklay',error)
615
      if (error.ne.0) goto 920
616
 
617
      call ncdinq(cstid,didz,dimnam,nlev,error)	! read number of levels
618
      if (error.ne.0) goto 920
619
 
620
      do 10 i=1,nlev
621
        call ncvgt1(cstid,idak,i,aklev(i),error)      ! get aklev
622
        call ncvgt1(cstid,idbk,i,bklev(i),error)      ! get bklev
623
        call ncvgt1(cstid,idaky,i,aklay(i),error)      ! get aklay
624
        call ncvgt1(cstid,idbky,i,bklay(i),error)      ! get bklay
625
        if (error.ne.0) goto 920
626
   10 continue
627
 
628
      return
629
 
630
c     Error exits.
631
  920 write(*,*)'*ERROR*: An error occured in subroutine getlevs'
632
      return
633
 
634
      end
635
      subroutine getntim(cdfid,ntimes,ierr)
636
C------------------------------------------------------------------------
637
C     Purpose:
638
C        Get number of times on the specified NetCDF file
639
C     Arguments:
640
C        cdfid  int  input   identifier for NetCDF file
641
C        ntimes int  output  number of times on the file
642
C        error  int  output  errorflag
643
C     History:
644
C        Heini Wernli, ETHZ
645
C------------------------------------------------------------------------
646
 
647
      include "netcdf.inc"
648
 
649
      integer   ierr
650
      integer didtim,ntimes
651
 
652
      integer   cdfid,idtime
653
      integer   ncopts
654
      character*(20) dimnam
655
 
656
c     Get current value of error options, and make sure netCDF-errors do
657
c     not abort execution
658
      call ncgopt (ncopts)
659
      call ncpopt(NCVERBOS)
660
 
661
      didtim=ncdid(cdfid,'time',ierr)   ! inquire id for time dimension
662
      if (ierr.ne.0) goto 900
663
      idtime=ncvid(cdfid,'time',ierr)   ! inquire id for time array
664
      if (ierr.ne.0) goto 900
665
      call ncdinq(cdfid,didtim,dimnam,ntimes,ierr)      ! inquire # of times
666
      if (ierr.ne.0) goto 900
667
 
668
c     normal exit
669
      call ncpopt (ncopts)
670
      return
671
 
672
c     error exit
673
 900  ntimes=1
674
      call ncpopt (ncopts)
675
      end
676
      subroutine getstart(cdfid,idate,ierr)
677
C------------------------------------------------------------------------
678
C     Purpose:
679
C	Get start date for fields on specified NetCDF file
680
C     Arguments:
681
C	cdfid	int	input	identifier for NetCDF file
682
C	idate	int	output	array contains date (year,month,day,time,step)
683
C				dimensioned as idate(5)
684
C	ierr	int	output	error flag
685
C------------------------------------------------------------------------
686
 
687
      include "netcdf.inc"
688
 
689
c     variable declarations
690
      integer   ierr
691
      integer   idate(5)
692
      integer   cdfid,ncopts,idvar,nvars
693
      integer   ndims,ngatts,recdim,i,vartyp,nvatts,vardim(4)
694
      character*20 vnam(100)
695
 
696
c     Get current value of error options, and make sure NetCDF-errors do
697
c     not abort execution
698
      call ncgopt (ncopts)
699
      call ncpopt (NCVERBOS)
700
 
701
      idvar=ncvid(cdfid,'starty',ierr)
702
      if (ierr.ne.0) goto 930
703
      call ncvgt1(cdfid,idvar,1,idate(1),ierr)
704
      if (ierr.ne.0) goto 920
705
 
706
      idvar=ncvid(cdfid,'startm',ierr)
707
      if (ierr.ne.0) goto 920
708
      call ncvgt1(cdfid,idvar,1,idate(2),ierr)
709
      if (ierr.ne.0) goto 920
710
 
711
      idvar=ncvid(cdfid,'startd',ierr)
712
      if (ierr.ne.0) goto920
713
      call ncvgt1(cdfid,idvar,1,idate(3),ierr)
714
      if (ierr.ne.0) goto 920
715
 
716
      idvar=ncvid(cdfid,'starth',ierr)
717
      if (ierr.ne.0) goto 920
718
      call ncvgt1(cdfid,idvar,1,idate(4),ierr)
719
      if (ierr.ne.0) goto 920
720
 
721
C     Starts is not defined on all files
722
C     Only ask for it if it exists
723
C     Inquire number of dimensions, variables and attributes
724
 
725
      idate(5)=0
726
      call ncinq(cdfid,ndims,nvars,ngatts,recdim,ierr)
727
      do i=1,nvars
728
        call ncvinq(cdfid,i,vnam(i),vartyp,ndims,vardim,nvatts,ierr)
729
        if (vnam(i).eq.'starts') then
730
          idvar=ncvid(cdfid,'starts',ierr)
731
          call ncvgt1(cdfid,idvar,1,idate(5),ierr)
732
          if (ierr.ne.0) goto 920
733
        endif
734
      enddo
735
 
736
c     normal exit
737
      call ncpopt (ncopts)
738
      return
739
 
740
c     error exit
741
 920  continue
742
      write (6, *) 'ERROR: An error occurred while attempting to ',
743
     &             'read the starting-time in subroutine putstart.'
744
 930  continue
745
      call ncpopt (ncopts)
746
 
747
      end
748
      subroutine putstart(cdfid,idate,ierr)
749
C----------------------------------------------------------------------
750
C     Purpose:
751
C        Puts the 'starting-time' on the specified NetCDF file.
752
C     Arguments:
753
C        cdfid   int     input   identifier for NetCDF file
754
C        idate   int     input   array contains date (year,month,day,time,step)
755
C                                dimensioned as idate(5)
756
C        ierr    int     output  error flag
757
C------------------------------------------------------------------------
758
 
759
      include "netcdf.inc"
760
 
761
c     variable declarations
762
      integer   ierr,idate(5),startid(5),cdfid,ncopts,i
763
 
764
c     Get current value of error options, and make sure NetCDF-errors do
765
c     not abort execution
766
      call ncgopt (ncopts)
767
      call ncpopt (NCVERBOS)
768
 
769
c     define variables
770
      call ncredf(cdfid,ierr)
771
      if (ierr.ne.0) goto 920
772
      startid(1) = ncvdef (cdfid, 'starty', NCLONG, 0, 0, ierr)
773
      if (ierr.ne.0) goto 920
774
      startid(2) = ncvdef (cdfid, 'startm', NCLONG, 0, 0, ierr)
775
      if (ierr.ne.0) goto 920
776
      startid(3) = ncvdef (cdfid, 'startd', NCLONG, 0, 0, ierr)
777
      if (ierr.ne.0) goto 920
778
      startid(4) = ncvdef (cdfid, 'starth', NCLONG, 0, 0, ierr)
779
      if (ierr.ne.0) goto 920
780
      startid(5) = ncvdef (cdfid, 'starts', NCLONG, 0, 0, ierr)
781
      if (ierr.ne.0) goto 920
782
      call ncendf(cdfid, ierr)
783
      if (ierr.ne.0) goto 920
784
 
785
c     store variables
786
      do i=1,5
787
        call ncvpt1(cdfid,startid(i),1,idate(i),ierr)
788
        if (ierr.ne.0) goto 920
789
      enddo
790
 
791
c     synchronyse output to disk, revert to previous error-mode, and exit
792
      call ncsnc (cdfid,ierr)
793
      call ncpopt (ncopts)
794
      return
795
 
796
c     error exit
797
 920  write (6, *) 'ERROR: An error occurred while attempting to ',
798
     &             'write the starting-time in subroutine putstart.'
799
      call ncpopt (ncopts)
800
      call ncclos (cdfid, ierr)
801
 
802
      end
803
      subroutine getgrid(cdfid,dx,dy,ierr)
804
C------------------------------------------------------------------------
805
C     Purpose:
806
C       Get grid increments for fields on specified NetCDF file
807
C     Arguments:
808
C       cdfid   int     input   identifier for NetCDF file
809
C	dx	real	output	grid increment along latitude
810
C	dy	real	output	grid increment along longitude
811
C       ierr    int     output  error flag
812
C------------------------------------------------------------------------
813
 
814
      integer   ierr
815
 
816
      integer   cdfid
817
      integer   ncvid
818
 
819
      integer   idilon,idilat
820
      real	dx,dy
821
 
822
      idilon    =ncvid(cdfid,'dellon',ierr)
823
      if (ierr.ne.0) return
824
      idilat    =ncvid(cdfid,'dellat',ierr)
825
      if (ierr.ne.0) return
826
 
827
      call ncvgt1(cdfid,idilon,1,dx,ierr)
828
      if (ierr.ne.0) return
829
      call ncvgt1(cdfid,idilat,1,dy,ierr)
830
      if (ierr.ne.0) return
831
 
832
      end
833
      subroutine getdattyp(cdfid,typ,ierr)
834
C------------------------------------------------------------------------
835
C     Purpose:
836
C       Get data type for specified NetCDF file
837
C     Arguments:
838
C       cdfid   int     input   identifier for NetCDF file
839
C       typ     int     output  data type: 1 (52) for pressure (theta) coord
840
C       ierr    int     output  error flag
841
C------------------------------------------------------------------------
842
 
843
      integer   ierr
844
 
845
      integer   cdfid
846
      integer   ncvid
847
 
848
      integer   idtyp,typ
849
 
850
      idtyp    =ncvid(cdfid,'dattyp',ierr)
851
      if (ierr.ne.0) return
852
 
853
      call ncvgt1(cdfid,idtyp,1,typ,ierr)
854
      if (ierr.ne.0) return
855
 
856
      end
857
      subroutine getpole(cdfid,pollon,pollat,ierr)
858
C------------------------------------------------------------------------
859
C     Purpose:
860
C       Get physical coordinates of pole of coordinate system
861
C     Arguments:
862
C       cdfid   int     input   identifier for NetCDF file
863
C	pollon	real	output	longitude of pole
864
C	pollat	real	output	latitude of pole
865
C       ierr    int     output  error flag
866
C------------------------------------------------------------------------
867
 
868
      integer   ierr
869
 
870
      integer   cdfid
871
      integer   ncvid
872
 
873
      integer   idplon,idplat
874
      real      pollon,pollat
875
 
876
      idplon    =ncvid(cdfid,'pollon',ierr)
877
      if (ierr.ne.0) return
878
      idplat    =ncvid(cdfid,'pollat',ierr)
879
      if (ierr.ne.0) return
880
 
881
      call ncvgt1(cdfid,idplon,1,pollon,ierr)
882
      if (ierr.ne.0) return
883
      call ncvgt1(cdfid,idplat,1,pollat,ierr)
884
      if (ierr.ne.0) return
885
 
886
      end
887
      subroutine getmc2grid(cdfid,polx,poly,delx,shem,phi0,lam0,ierr)
888
C------------------------------------------------------------------------
889
C     Purpose:
890
C       Get physical coordinates of pole of coordinate system
891
C     Arguments:
892
C       cdfid   int     input   identifier for NetCDF file
893
C       ierr    int     output  error flag
894
C------------------------------------------------------------------------
895
 
896
      integer   ierr
897
 
898
      integer   cdfid
899
      integer   ncvid
900
 
901
      integer   idpolx,idpoly,iddelx,idshem,idphi0,idlam0
902
      real      polx,poly,delx,shem,phi0,lam0
903
 
904
      idpolx    =ncvid(cdfid,'polx',ierr)
905
      if (ierr.ne.0) return
906
      idpoly    =ncvid(cdfid,'poly',ierr)
907
      if (ierr.ne.0) return
908
      iddelx	=ncvid(cdfid,'delx',ierr)
909
      if (ierr.ne.0) return
910
      idshem	=ncvid(cdfid,'shem',ierr)
911
      if (ierr.ne.0) return
912
      idphi0	=ncvid(cdfid,'phi0',ierr)
913
      if (ierr.ne.0) return
914
      idlam0	=ncvid(cdfid,'lam0',ierr)
915
      if (ierr.ne.0) return
916
 
917
      call ncvgt1(cdfid,idpolx,1,polx,ierr)
918
      if (ierr.ne.0) return
919
      call ncvgt1(cdfid,idpoly,1,poly,ierr)
920
      if (ierr.ne.0) return
921
      call ncvgt1(cdfid,iddelx,1,delx,ierr)
922
      if (ierr.ne.0) return
923
      call ncvgt1(cdfid,idshem,1,shem,ierr)
924
      if (ierr.ne.0) return
925
      call ncvgt1(cdfid,idphi0,1,phi0,ierr)
926
      if (ierr.ne.0) return
927
      call ncvgt1(cdfid,idlam0,1,lam0,ierr)
928
      if (ierr.ne.0) return
929
 
930
      end
931
      subroutine getcfn(cdfid,cfn,ierr)
932
C------------------------------------------------------------------------
933
C     Purpose:
934
C       Get name of constants file
935
C     Arguments:
936
C       cdfid   int     input   identifier for NetCDF file
937
C       cfn     char    output  name of constants file
938
C       ierr    int     output  error flag
939
C------------------------------------------------------------------------
940
 
941
      include "netcdf.inc"
942
 
943
      integer   ierr
944
      integer   cdfid,lenstr
945
      character*80 cfn
946
 
947
      lenstr=80
948
      call ncagtc(cdfid,NCGLOBAL,"constants_file_name",cfn,lenstr,ierr)
949
      if (ierr.ne.0) write(*,*)'error in SR getcfn'
950
 
951
      end
952
      subroutine gettype(cdfid,dattyp,datver,cstver,ierr)
953
C------------------------------------------------------------------------
954
C     Purpose:
955
C       Get data type information from constants file
956
C     Arguments:
957
C       cdfid   int     input   identifier for NetCDF file
958
C       dattyp  int	output  data type
959
C       datver  int	output  data version
960
C       cstver	int     output  constants file version
961
C------------------------------------------------------------------------
962
 
963
      integer   ierr
964
 
965
      integer   cdfid
966
      integer   ncvid
967
 
968
      integer   idtyp,idver,idcstv
969
      integer	dattyp,datver,cstver
970
 
971
      idtyp	=ncvid(cdfid,'dattyp',ierr)
972
      if (ierr.ne.0) return
973
      idver	=ncvid(cdfid,'datver',ierr)
974
      if (ierr.ne.0) return
975
      idcstv    =ncvid(cdfid,'cstver',ierr)
976
      if (ierr.ne.0) return
977
 
978
      call ncvgt1(cdfid,idtyp,1,dattyp,ierr)
979
      if (ierr.ne.0) return
980
      call ncvgt1(cdfid,idver,1,datver,ierr)
981
      if (ierr.ne.0) return
982
      call ncvgt1(cdfid,idcstv,1,cstver,ierr)
983
      if (ierr.ne.0) return
984
 
985
      end
986
      subroutine getvars(cdfid,nvars,vnam,ierr)
987
C------------------------------------------------------------------------
988
 
989
C     Opens the NetCDF file 'filnam' and returns its identifier cdfid.
990
 
991
C     filnam    char    input   name of NetCDF file to open
992
C     nvars     int     output  number of variables on file
993
C     vnam	char	output  array with variable names
994
C     ierr      int     output  error flag
995
C------------------------------------------------------------------------
996
 
997
      include "netcdf.inc"
998
 
999
      integer   cdfid,ierr,nvars
1000
      character*(*) vnam(*)
1001
 
1002
      integer	ndims,ngatts,recdim,i,vartyp,nvatts,vardim(4)
1003
 
1004
      call ncpopt(NCVERBOS)
1005
 
1006
C     Inquire number of dimensions, variables and attributes
1007
 
1008
      call ncinq(cdfid,ndims,nvars,ngatts,recdim,ierr)
1009
 
1010
C     Inquire variable names from NetCDF file
1011
 
1012
      do i=1,nvars
1013
        call ncvinq(cdfid,i,vnam(i),vartyp,ndims,vardim,nvatts,ierr)
1014
      enddo
1015
 
1016
      return
1017
      end
1018
 
1019
      subroutine cdfopn(filnam,cdfid,ierr)
1020
C------------------------------------------------------------------------
1021
 
1022
C     Opens the NetCDF file 'filnam' and returns its identifier cdfid.
1023
 
1024
C     filnam    char    input   name of NetCDF file to open
1025
C     cdfid     int     output  identifier of NetCDF file
1026
C     ierr	int	output  error flag
1027
C------------------------------------------------------------------------
1028
 
1029
      include "netcdf.inc"
1030
 
1031
      integer 	cdfid,ierr
1032
      character*(*) filnam
1033
 
1034
      call ncpopt(NCVERBOS)
1035
      cdfid=ncopn(trim(filnam),NCNOWRIT,ierr)
1036
 
1037
      return
1038
      end
1039
      subroutine cdfwopn(filnam,cdfid,ierr)
1040
C------------------------------------------------------------------------
1041
 
1042
C     Opens the NetCDF file 'filnam' and returns its identifier cdfid.
1043
 
1044
C     filnam    char    input   name of NetCDF file to open
1045
C     cdfid     int     output  identifier of NetCDF file
1046
C     ierr      int     output  error flag
1047
C------------------------------------------------------------------------
1048
 
1049
      include "netcdf.inc"
1050
 
1051
      integer   cdfid,ierr
1052
      character*(*) filnam
1053
 
1054
      call ncpopt(NCVERBOS)
1055
      cdfid=ncopn(trim(filnam),NCWRITE,ierr)
1056
 
1057
      return
1058
      end
1059
      subroutine gettra(cdfid,varnam,ix,iy,iz,ntimes,array,ierr)
1060
C------------------------------------------------------------------------
1061
C
1062
C     Reads the time-evolution for one grid-point of the variable
1063
C     indicated by varnam.
1064
C
1065
C     cdfid     int     input   identifier for NetCDF file
1066
C     varnam    char    input   name of variable
1067
C     ix        int     input   x-index for values to read
1068
C     iy        int     input   y-index for values to read
1069
C     iz        int     input   z-index for values to read
1070
C     ntimes    int     input   number of time-indices to read
1071
C     array     real    output  array contains the readed values
1072
C     ierr      int     output  error flag
1073
C------------------------------------------------------------------------
1074
 
1075
C     Declaration of attributes
1076
 
1077
      integer   cdfid
1078
      character*(*) varnam
1079
      integer   ix,iy,iz
1080
      integer	ntimes
1081
      real      array(ntimes)
1082
 
1083
C     Declaration of local variables
1084
 
1085
      integer   corner(4),edgeln(4)
1086
      integer   idvar,ierr
1087
      integer	ncvid
1088
 
1089
      corner(1)=ix
1090
      corner(2)=iy
1091
      corner(3)=iz
1092
      corner(4)=1
1093
      edgeln(1)=1
1094
      edgeln(2)=1
1095
      edgeln(3)=1
1096
      edgeln(4)=ntimes
1097
 
1098
      idvar =ncvid(cdfid,varnam,ierr)
1099
      call ncvgt(cdfid,idvar,corner,edgeln,array,ierr)
1100
      if (ierr.ne.0) goto 991
1101
 
1102
      return
1103
  991 stop 'Variable not found on NetCDF file in SR gettra'
1104
      end
1105
      subroutine new_gettra(cdfid,varnam,ix,ntimes,array,ierr)
1106
C------------------------------------------------------------------------
1107
C
1108
C     Reads the time-evolution for one grid-point of the variable
1109
C     indicated by varnam.
1110
C
1111
C     cdfid     int     input   identifier for NetCDF file
1112
C     varnam    char    input   name of variable
1113
C     ix        int     input   index for trajectory to read
1114
C     ntimes    int     input   number of time-indices to read
1115
C     array     real    output  array contains the readed values
1116
C     ierr      int     output  error flag
1117
C------------------------------------------------------------------------
1118
 
1119
C     Declaration of attributes
1120
 
1121
      integer   cdfid
1122
      character*(*) varnam
1123
      integer   ix
1124
      integer   ntimes
1125
      real      array(ntimes)
1126
 
1127
C     Declaration of local variables
1128
 
1129
      integer   corner(4),edgeln(4)
1130
      integer   idvar,ierr
1131
      integer   ncvid
1132
 
1133
      corner(1)=ix
1134
      corner(2)=1
1135
      corner(3)=1
1136
      corner(4)=1
1137
      edgeln(1)=1
1138
      edgeln(2)=1
1139
      edgeln(3)=1
1140
      edgeln(4)=ntimes
1141
 
1142
      idvar =ncvid(cdfid,trim(varnam),ierr)
1143
      call ncvgt(cdfid,idvar,corner,edgeln,array,ierr)
1144
      if (ierr.ne.0) goto 991
1145
 
1146
      return
1147
  991 stop 'Variable not found on NetCDF file in SR new_gettra'
1148
      end
1149
      subroutine puttra(cdfid,varnam,ix,ntimes,array,ierr)
1150
C------------------------------------------------------------------------
1151
C
1152
C     Writes the time-evolution for one grid-point of the variable
1153
C     indicated by varnam.
1154
C
1155
C     cdfid     int     input   identifier for NetCDF file
1156
C     varnam    char    input   name of variable
1157
C     ix        int     input   index for trajectory to read
1158
C     ntimes    int     input   number of time-indices to read
1159
C     array     real    output  array contains the readed values
1160
C     ierr      int     output  error flag
1161
C------------------------------------------------------------------------
1162
 
1163
C     Declaration of attributes
1164
 
1165
      integer   cdfid
1166
      character*(*) varnam
1167
      integer   ix
1168
      integer   ntimes
1169
      real      array(ntimes)
1170
 
1171
C     Declaration of local variables
1172
 
1173
      integer   corner(4),edgeln(4)
1174
      integer   idvar,ierr
1175
      integer   ncvid
1176
 
1177
      corner(1)=1
1178
      corner(2)=1
1179
      corner(3)=1
1180
      corner(4)=ix
1181
      edgeln(1)=ntimes
1182
      edgeln(2)=1
1183
      edgeln(3)=1
1184
      edgeln(4)=1
1185
 
1186
      idvar =ncvid(cdfid,varnam,ierr)
1187
      call ncvpt(cdfid,idvar,corner,edgeln,array,ierr)
1188
      if (ierr.ne.0) goto 991
1189
 
1190
      return
1191
  991 stop 'Could not write data on NetCDF file in SR puttra'
1192
      end
1193
      subroutine getakbk(nlev,flev,akbk,nn,aklev,bklev,aklay,bklay)
1194
C------------------------------------------------------------------------
1195
C
1196
C     Defines the level- and layer-arrays given the number of levels nlev.
1197
C
1198
C     nlev      int     input   number of levels/layers wanted
1199
C     akbk	real	input	array contains ak/bk values from grib (zsec2)
1200
C     nn	int	input	number of elements in array akbk
1201
C     aklev     real    output  array contains ak values for levels
1202
C     bklev     real    output  array contains bk values for levels
1203
C     aklay     real    output  array contains ak values for layers
1204
C     bklay     real    output  array contains bk values for layers
1205
C------------------------------------------------------------------------
1206
 
1207
      integer   nn,nz,nlev,k
1208
      real      aklev(100),bklev(100),    ! level coefficients
1209
     >          aklay(100),bklay(100),    ! layer coefficients
1210
     >		akbk(nn)
1211
      real	ak(100),bk(100)
1212
      real	flev
1213
 
1214
C     Determine number of levels in array akbk
1215
      do k=1,nn
1216
        if (akbk(k).eq.1.0) nz=(k-12)/2
1217
      enddo 
1218
c      print*,nlev,nz
1219
 
1220
      do k=1,nz+1
1221
        ak(k)=akbk(k+10)/100.
1222
        bk(k)=akbk(k+11+nz)
1223
      enddo
1224
 
1225
      do k=1,nz
1226
        aklay(k)=(ak(nz+2-k)+ak(nz+1-k))/2.
1227
        bklay(k)=(bk(nz+2-k)+bk(nz+1-k))/2.
1228
        aklev(k)=ak(nz+1-k)
1229
        bklev(k)=bk(nz+1-k)
1230
c        if (k.eq.2) print*,'bugfix ',bklev(2)
1231
      enddo
1232
 
1233
c      do k=1,nz
1234
c        print*,k,flev,bk(nz+1-k),aklev(k),aklay(k),bklev(k),bklay(k)
1235
c      enddo
1236
 
1237
      return
1238
      end
1239
      subroutine modlevs(nlev,aklev,bklev,aklay,bklay)
1240
C------------------------------------------------------------------------
1241
C
1242
C     Defines the level- and layer-arrays given the number of levels nlev.
1243
C
1244
C     nlev	int	input	number of levels/layers
1245
C     aklev	real	output	array contains ak values for levels
1246
C     bklev     real    output  array contains bk values for levels
1247
C     aklay     real    output  array contains ak values for layers
1248
C     bklay     real    output  array contains bk values for layers
1249
C------------------------------------------------------------------------
1250
 
1251
      integer   n19,n31,n50,nlev,k
1252
      parameter(n19=20,n31=32,n50=51)           ! number of model levels
1253
      real      aklev(nlev+1),bklev(nlev+1),    ! level coefficients
1254
     >          aklay(nlev+1),bklay(nlev+1)     ! layer coefficients
1255
 
1256
      real      ak19(n19),bk19(n19),            ! 19 level version
1257
     >          ak31(n31),bk31(n31),            ! 31 level version
1258
     >		ak50(n50),bk50(n50)             ! 50 level version
1259
 
1260
C     Modell level specification for 19 level version
1261
      DATA AK19/0,20,40,60,83,106,128,146,158,161,153,136,111,
1262
     >        82,52,26,8,0,0,0/
1263
      DATA BK19/0,0,0,0,.004,.014,.035,.072,.127,.202,.296,.405,
1264
     >        .524,.645,.759,.856,.929,.973,.992,1./
1265
 
1266
 
1267
C     Modell level specification for 31 level version
1268
      DATA AK31/
1269
     >   0.000000,  20.00000000,  40.00000000,  60.00000000,
1270
     >  80.000000,  99.76135361, 118.20539617, 134.31393926,
1271
     > 147.363569, 156.89207458, 162.66610500, 164.65005734,
1272
     > 162.976193, 157.91598604, 149.85269630, 139.25517858,
1273
     > 126.652916, 112.61228878,  97.71406290,  82.53212096,
1274
     >  67.613413,  53.45914240,  40.50717678,  29.11569385,
1275
     >  19.548052,  11.95889791,   6.38148911,   2.71626545,
1276
     >    .720635,   0.00000000,   0.00000000,   0.00000000/
1277
 
1278
      DATA BK31/
1279
     >   0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,
1280
     >   0.0000000000, 0.0003908582, 0.0029197006, 0.0091941320,
1281
     >   0.0203191555, 0.0369748598, 0.0594876397, 0.0878949492,
1282
     >   0.1220035886, 0.1614415235, 0.2057032385, 0.2541886223,
1283
     >   0.3062353873, 0.3611450218, 0.4182022749, 0.4766881754,
1284
     >   0.5358865832, 0.5950842740, 0.6535645569, 0.7105944258,
1285
     >   0.7654052430, 0.8171669567, 0.8649558510, 0.9077158297,
1286
     >   0.9442132326, 0.9729851852, 0.9922814815, 1.0000000000/
1287
 
1288
C     Modell level specification for 50 level version
1289
      DATA AK50/
1290
     >     0.0000,    .200061,    .432978,
1291
     >    .753462,   1.150821,   1.618974,   2.158969,
1292
     >   2.780058,   3.501381,   4.355622,   5.396513,
1293
     >   6.686154,   8.283989,  10.263669,  12.716445,
1294
     >  15.755378,  19.520544,  24.185498,  29.965266,
1295
     >  37.126262,  45.998554,  56.991132,  69.983867,
1296
     >  85.074101, 101.817070, 118.830898, 134.429140,
1297
     > 147.363554, 156.892070, 162.666093, 164.650039,
1298
     > 162.976210, 157.915976, 149.852695, 139.255195,
1299
     > 126.652968, 112.612304,  97.714062,  82.532109,
1300
     >  67.613398,  53.459179,  40.507187,  29.115703,
1301
     >  19.548046,  11.958906,   6.381484,   2.716250,
1302
     >    .720625,   0.000000,   0.000000,   0.000000/
1303
 
1304
      DATA BK50/
1305
     >   0.0000000000, 0.0000000000, 0.0000000000,
1306
     >   0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,
1307
     >   0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,
1308
     >   0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,
1309
     >   0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,
1310
     >   0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,
1311
     >   0.0001003604, 0.0006727143, 0.0031633405, 0.0092923380,
1312
     >   0.0203191563, 0.0369748585, 0.0594876409, 0.0878949761,
1313
     >   0.1220036149, 0.1614415050, 0.2057032585, 0.2541885972,
1314
     >   0.3062353730, 0.3611450195, 0.4182022810, 0.4766881466,
1315
     >   0.5358865857, 0.5950842500, 0.6535645723, 0.7105944157,
1316
     >   0.7654052377, 0.8171669841, 0.8649558425, 0.9077158570,
1317
     >   0.9442132115, 0.9729852080, 0.9922814965, 1.0000000000/
1318
 
1319
      do k=1,nlev
1320
        if (nlev.eq.19) then
1321
          aklay(k)=(ak19(nlev+2-k)+ak19(nlev+1-k))/2.   ! layi=(levi+levi+1)/2
1322
          bklay(k)=(bk19(nlev+2-k)+bk19(nlev+1-k))/2.
1323
          aklev(k)=ak19(nlev+1-k)       ! reverse order of coeffs for IVE
1324
          bklev(k)=bk19(nlev+1-k)
1325
        elseif (nlev.eq.31) then
1326
          aklay(k)=(ak31(nlev+2-k)+ak31(nlev+1-k))/2.   ! layi=(levi+levi+1)/2
1327
          bklay(k)=(bk31(nlev+2-k)+bk31(nlev+1-k))/2.
1328
          aklev(k)=ak31(nlev+1-k)       ! reverse order of coeffs for IVE
1329
          bklev(k)=bk31(nlev+1-k)
1330
        elseif (nlev.eq.50) then
1331
          aklay(k)=(ak50(nlev+2-k)+ak50(nlev+1-k))/2.   ! layi=(levi+levi+1)/2
1332
          bklay(k)=(bk50(nlev+2-k)+bk50(nlev+1-k))/2.
1333
          aklev(k)=ak50(nlev+1-k)       ! reverse order of coeffs for IVE
1334
          bklev(k)=bk50(nlev+1-k)
1335
        else
1336
          stop'*** invalid number of modellevels ***'
1337
        endif
1338
      enddo
1339
 
1340
      if (nlev.eq.19) then
1341
        aklay(nlev+1)=ak19(1)/2.
1342
        bklay(nlev+1)=bk19(1)/2.
1343
        aklev(nlev+1)=ak19(1)
1344
        bklev(nlev+1)=bk19(1)
1345
      elseif (nlev.eq.31) then
1346
        aklay(nlev+1)=ak31(1)/2.
1347
        bklay(nlev+1)=bk31(1)/2.
1348
        aklev(nlev+1)=ak31(1)
1349
        bklev(nlev+1)=bk31(1)
1350
      elseif (nlev.eq.50) then
1351
        aklay(nlev+1)=ak50(1)/2.
1352
        bklay(nlev+1)=bk50(1)/2.
1353
        aklev(nlev+1)=ak50(1)
1354
        bklev(nlev+1)=bk50(1)
1355
      else
1356
        stop'*** invalid number of modellevels ***'
1357
      endif
1358
*     print*,aklev(1),aklev(2),aklev(3),aklev(4),aklev(5),aklev(6)
1359
 
1360
      return
1361
      end
1362
 
1363
      subroutine prelevs(nlev,level,aklev,bklev,aklay,bklay)
1364
C------------------------------------------------------------------------
1365
C
1366
C     Defines the (dummy-) ak- and bk-arrays given the array that
1367
C     contains all pressure levels.
1368
C
1369
C     nlev	int	input	number of pressure levels
1370
C     level	real	input	pressure levels
1371
C     aklev     real    output  array contains ak values for levels
1372
C     bklev     real    output  array contains bk values for levels
1373
C     aklay     real    output  array contains ak values for layers
1374
C     bklay     real    output  array contains bk values for layers
1375
C------------------------------------------------------------------------
1376
 
1377
      integer   nlev,k
1378
      real      aklev(nlev),bklev(nlev),	! level coefficients
1379
     >          aklay(nlev),bklay(nlev),	! layer coefficients
1380
     >		level(nlev+1)
1381
 
1382
      do k=1,nlev
1383
        aklay(k)=level(k)
1384
        bklay(k)=0.
1385
        if (nlev.eq.1) then
1386
          aklev(k)=level(k)
1387
        else
1388
          aklev(k)=0.5*(level(k)+level(k+1))
1389
        endif
1390
        bklev(k)=0.
1391
      enddo
1392
 
1393
      return
1394
      end
1395
 
1396
 
1397
 
1398
      subroutine cpp_cdfwopn(filnam,filnam_len,cdfid,ierr)
1399
C------------------------------------------------------------------------
1400
C     Purpose:
1401
C        allows to call cdfopn from c++
1402
C     Arguments: 
1403
C        see crecdf
1404
C        additionally: filnam_len, the length of the 
1405
C           string
1406
C        
1407
C        
1408
C     History:
1409
C        981221  Mark A. Liniger ETHZ
1410
C        
1411
C     Note:
1412
C        
1413
C        
1414
C------------------------------------------------------------------------
1415
      integer        filnam_len,cdfid,ierr
1416
      character *(*) filnam
1417
 
1418
 
1419
      call cdfwopn(filnam(1:filnam_len),cdfid,ierr)
1420
 
1421
      end
1422
      subroutine getdim (cdfid, varnam, nx, ny, nz, error)
1423
c-------------------------------------------------------------------------
1424
c     Purpose:
1425
c        This routine is called to get the dimensions of
1426
c        a variable from an IVE-NetCDF file for use with the IVE plotting
1427
c        package. Prior to calling this routine, the file must be opened
1428
c        with a call to opncdf.
1429
c     Arguments:
1430
c        cdfid   int   input   file-identifier
1431
c                              (can be obtained by calling routine
1432
c                              opncdf)
1433
c        varnam  char  input   the user-supplied variable name.
1434
c                              (can be obtained by calling routine
1435
c                              opncdf)
1436
c        nx      int   output  the zonal dimension of the variable.
1437
c        ny      int   output  the meridional dimension of the variable.
1438
c        nz      int   output  the vertical dimension of the variable.
1439
c
1440
c        error   int   output  indicates possible errors found in this
1441
c                              routine.
1442
c                              error = 0   no errors detected.
1443
c                              error = 1   the variable is not on the file.
1444
c                              error =10   other errors.
1445
c     History:
1446
c        March 2000    Heini Wernli (ETHZ)     Created.
1447
c-----------------------------------------------------------------------
1448
 
1449
      include "netcdf.inc"
1450
 
1451
c     Argument declarations.
1452
      character *(*) varnam
1453
      integer        vardim(4), ndim, error, cdfid
1454
      integer        nx,ny,nz
1455
 
1456
c     Local variable declarations.
1457
      character *(20) dimnam(MAXNCDIM),vnam
1458
      integer         id,i,k
1459
      integer         ndims,nvars,ngatts,recdim,dimsiz(MAXNCDIM)
1460
      integer         vartyp,nvatts, ncopts
1461
 
1462
c     Get current value of error options.
1463
      call ncgopt (ncopts)
1464
 
1465
c     make sure NetCDF-errors do not abort execution
1466
      call ncpopt(NCVERBOS)
1467
 
1468
c     Initially set error to indicate no errors.
1469
      error = 0
1470
 
1471
c     inquire for number of dimensions
1472
      call ncinq(cdfid,ndims,nvars,ngatts,recdim,error)
1473
      if (error.eq.1) goto 920
1474
 
1475
c     read dimension-table
1476
      do i=1,ndims
1477
        call ncdinq(cdfid,i,dimnam(i),dimsiz(i),error)
1478
        if (error.gt.0) goto 920
1479
      enddo
1480
 
1481
c     get id of the variable
1482
      id=ncvid(cdfid,varnam,error)
1483
      if (error.eq.1) goto 910
1484
 
1485
c     inquire about variable
1486
      call ncvinq(cdfid,id,vnam,vartyp,ndim,vardim,nvatts,error)
1487
      if (vartyp.ne.NCFLOAT) error=1
1488
      if (error.gt.0) goto 920
1489
 
1490
c     get dimensions from dimension-table
1491
      do k=1,ndim
1492
        vardim(k)=dimsiz(vardim(k))
1493
      enddo
1494
 
1495
      nx=vardim(1)
1496
      ny=vardim(2)
1497
      nz=vardim(3)
1498
 
1499
c     normal exit
1500
      call ncpopt (ncopts)
1501
      return
1502
 
1503
c     Error exits.
1504
 910  write (6, *) '*ERROR*: The selected variable could not be found ',
1505
     &             'in the file by getdim.'
1506
      call ncpopt (ncopts)
1507
      call ncclos (cdfid, error)
1508
      return
1509
 
1510
 920  write (6, *) '*ERROR*: An error occurred while attempting to ',
1511
     &             'read the data file in subroutine getcdf.'
1512
      call ncpopt (ncopts)
1513
      call ncclos (cdfid, error)
1514
      end
1515
      subroutine rvarfile(vnam,gribnr,levty,unit,factor,bias,
1516
     >                    lnum,stg,tdep,p,lval,varcnt,ierr)
1517
C     =======================================================
1518
C     Variablen-File in Arrays einlesen
1519
 
1520
      integer   maxvar
1521
      parameter(maxvar=100)
1522
 
1523
      character*(15) vnam(maxvar)
1524
      character*(13) unit(maxvar)
1525
      character*(1)  flag
1526
      integer   gribnr(maxvar),levty(maxvar),lnum(maxvar),
1527
     >          stg(maxvar),tdep(maxvar),p(maxvar),lval(maxvar)
1528
      real      factor(maxvar),bias(maxvar)
1529
 
1530
      integer   i,varcnt,ierr,nt
1531
 
1532
      nt=14             ! number of tape
1533
      i=1               ! initialize var-counter
1534
 
1535
C     Read first character of row and decide if it is comment or not
1536
 
1537
  100 read(nt,10,err=123,end=126) flag
1538
      if (flag.eq."#") goto 100         ! don't bother about comments
1539
      backspace nt
1540
  121 read(nt,122, err=123, end=126) vnam(i), gribnr(i), levty(i),
1541
     & unit(i), factor(i), bias(i), lnum(i), stg(i), tdep(i), p(i),
1542
     & lval(i)
1543
      i=i+1
1544
*     goto 100
1545
      goto 121
1546
 
1547
   10 format(a1)
1548
  122 format(a14,i3,i11,a17,f7.5,f9.2,i7,i4,i6,i3,i5)
1549
* 123 print *,'*ERROR* in subroutine rvarfile'
1550
  123 goto 121
1551
  126 continue
1552
      varcnt=i-1        ! # of variables in varfile_i
1553
 
1554
C     Check some things
1555
 
1556
      ierr=0            ! initialize error flag
1557
      do i=1,varcnt
1558
        if ((lnum(i).ne.1).and.(lnum(i).ne.2).and.(lnum(i).ne.3)
1559
     >      .and.(lnum(i).ne.4)) ierr=11
1560
        if ((stg(i).ne.0).and.(stg(i).ne.1).and.(stg(i).ne.10).and.
1561
     >      (stg(i).ne.11)) ierr=12
1562
        if ((tdep(i).ne.0).and.(tdep(i).ne.1)) ierr=13
1563
        if ((p(i).ne.0).and.(p(i).ne.1).and.(p(i).ne.2)) ierr=14
1564
        if ((lval(i).lt.0).or.(lval(i).gt.1050)) ierr=15
1565
      enddo
1566
 
1567
      return
1568
      end