Subversion Repositories pvinversion.ecmwf

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
      PROGRAM p2z
2
 
3
c     Calculate the geopotential and add it to the P file
4
c     Michael Sprenger / Spring 2006
5
 
6
      implicit none
7
 
8
c     ---------------------------------------------------------------
9
c     Declaration of variables
10
c     ---------------------------------------------------------------
11
 
12
c     Variables for input P file : model level
13
      character*80 ml_pfn
14
      real         ml_varmin(4),ml_varmax(4),ml_stag(4)
15
      integer      ml_vardim(4)
16
      real         ml_mdv
17
      integer      ml_ndim
18
      integer      ml_nx,ml_ny,ml_nz
19
      real         ml_xmin,ml_xmax,ml_ymin,ml_ymax,ml_dx,ml_dy
20
      integer      ml_ntimes
21
      real         ml_aklev(500),ml_bklev(500)
22
      real         ml_aklay(500),ml_bklay(500)
23
      real         ml_time
24
      real         ml_pollon,ml_pollat
25
      integer      ml_nvars
26
      character*80 ml_vnam(100)
27
      integer      ml_idate(5)
28
      real,allocatable, dimension (:,:)   :: ml_ps,ml_zb
29
      real,allocatable, dimension (:,:,:) :: ml_t3,ml_q3,ml_p3,ml_tv3
30
      real,allocatable, dimension (:,:,:) :: ml_z3
31
 
32
c     Variables for input Z file : pressure level
33
      character*80 pl_zfn
34
      real         pl_varmin(4),pl_varmax(4),pl_stag(4)
35
      integer      pl_vardim(4)
36
      real         pl_mdv
37
      integer      pl_ndim
38
      integer      pl_nx,pl_ny,pl_nz 
39
      real         pl_xmin,pl_xmax,pl_ymin,pl_ymax,pl_dx,pl_dy
40
      integer      pl_ntimes
41
      real         pl_aklev(500),pl_bklev(500)
42
      real         pl_aklay(500),pl_bklay(500)
43
      real         pl_time
44
      real         pl_pollon,pl_pollat
45
      integer      pl_nvars
46
      character*80 pl_vnam(100)
47
      integer      pl_idate(5)
48
      real,allocatable, dimension (:,:,:) :: pl_z3,pl_p3      
49
 
50
c     Variables for output Z file  : height level
51
      character*80 zl_ofn
52
      real         zl_varmin(4),zl_varmax(4),zl_stag(4)
53
      integer      zl_vardim(4)
54
      real         zl_mdv
55
      integer      zl_ndim
56
      integer      zl_nx,zl_ny,zl_nz 
57
      real         zl_xmin,zl_xmax,zl_ymin,zl_ymax,zl_dx,zl_dy
58
      integer      zl_ntimes
59
      real         zl_aklev(500),zl_bklev(500)
60
      real         zl_aklay(500),zl_bklay(500)
61
      real         zl_time
62
      real         zl_pollon,zl_pollat
63
      integer      zl_idate(5)    
64
      real,allocatable, dimension (:,:,:) :: zl_field
65
      real,allocatable, dimension (:,:,:) :: zl_p
66
 
67
 
68
c     Variables for input P,S file  : model level
69
      character*80 in_ofn
70
      real         in_varmin(4),in_varmax(4),in_stag(4)
71
      integer      in_vardim(4)
72
      real         in_mdv
73
      integer      in_ndim
74
      integer      in_nx,in_ny,in_nz 
75
      real         in_xmin,in_xmax,in_ymin,in_ymax,in_dx,in_dy
76
      integer      in_ntimes
77
      real         in_aklev(500),in_bklev(500)
78
      real         in_aklay(500),in_bklay(500)
79
      real         in_time
80
      real         in_pollon,in_pollat
81
      integer      in_nvars
82
      character*80 in_vnam(100)
83
      integer      in_idate(5)  
84
      real,allocatable, dimension (:,:,:) :: in_field
85
 
86
c     Physical and numerical parameters
87
      real      g
88
      parameter (g=9.80665)
89
      real      eps
90
      parameter (eps=0.01)
91
      real      tzero
92
      parameter (tzero=273.15)
93
      real      kappa
94
      parameter (kappa=0.6078)
95
      real      zerodiv
96
      parameter (zerodiv=0.0000000001)
97
      real      dpmin
98
      parameter (dpmin=10.)
99
      real      rdg
100
      parameter (rdg=29.271)
101
 
102
c     Flag for test mode
103
      integer      test
104
      parameter    (test=0)
105
      character*80 testfile
106
      parameter    (testfile='TEST')
107
 
108
c     Variables and levels for interpolation onto z levels
109
      character*80 levelfile,varfile
110
      integer      nvar,nlev
111
      character*80 varinp(100),varout(100),varsrc(100)
112
      real         zlev(500)
113
 
114
c     Auxiliary variables
115
      integer      ierr
116
      integer      cdfid,cstid
117
      character*80 cfn
118
      integer      stat
119
      real         time
120
      real         tv1(1000),z1(1000),p1(1000),f1(1000)
121
      real         spline_tv1(1000),spline_f1(1000),spline_z1(1000)
122
      real         pu,po,zu,zo,p,z,dp,p0,tvu,tvo,ff
123
      integer      i,j,k,l
124
      integer      lmin,n
125
      character*80 varname,cdfname
126
      integer      idate(5),stdate(5),datar(14)
127
      integer      isok
128
      character*80 name
129
      real         zmin,dz
130
 
131
c     -----------------------------------------------------------------
132
c     Read input fields
133
c     -----------------------------------------------------------------
134
 
135
      print*,'*********************************************************'
136
      print*,'* p2z                                                   *'
137
      print*,'*********************************************************'
138
 
139
c     Read in the parameter file
140
      open(10,file='fort.10')
141
 
142
       read(10,*) ml_pfn
143
       read(10,*) pl_zfn
144
       read(10,*) zl_ofn
145
 
146
       read(10,*) name,zmin
147
       if ( trim(name).ne.'GEO_ZMIN') stop
148
       read(10,*) name,nlev
149
       if ( trim(name).ne.'GEO_NZ'  ) stop
150
       read(10,*) name,dz
151
       if ( trim(name).ne.'GEO_DZ'  ) stop
152
       do i=1,nlev
153
          zlev(i)=zmin+real(i-1)*dz
154
       enddo
155
 
156
       nvar=1
157
 102   continue
158
        read(10,*,end=103) varinp(nvar),varout(nvar),varsrc(nvar)
159
        nvar=nvar+1
160
        goto 102
161
 103   continue
162
      nvar=nvar-1
163
      print*
164
      do i=1,nvar
165
        write(*,'(a10,a10,a30)') trim(varinp(i)),trim(varout(i)),
166
     >                           trim(varsrc(i))
167
      enddo
168
      print*
169
 
170
      close(10)
171
 
172
 
173
c     Get grid description for P file : model level
174
      call cdfopn(ml_pfn,cdfid,ierr)
175
      if (ierr.ne.0) goto 998
176
      call getcfn(cdfid,cfn,ierr)
177
      if (ierr.ne.0) goto 998
178
      call cdfopn(cfn,cstid,ierr)   
179
      if (ierr.ne.0) goto 998
180
      call getvars(cdfid,ml_nvars,ml_vnam,ierr)
181
      varname='T'
182
      isok=0
183
      call check_varok (isok,varname,ml_vnam,ml_nvars)
184
      if (isok.ne.1) goto 998
185
      call getdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
186
     >            ml_varmin,ml_varmax,ml_stag,ierr)
187
      if (ierr.ne.0) goto 998
188
      ml_nx  =ml_vardim(1)
189
      ml_ny  =ml_vardim(2)
190
      ml_nz  =ml_vardim(3)
191
      ml_xmin=ml_varmin(1)
192
      ml_ymin=ml_varmin(2)
193
      call getlevs(cstid,ml_nz,ml_aklev,ml_bklev,ml_aklay,ml_bklay,ierr)
194
      call getgrid(cstid,ml_dx,ml_dy,ierr)
195
      ml_xmax=ml_xmin+real(ml_nx-1)*ml_dx
196
      ml_ymax=ml_ymin+real(ml_ny-1)*ml_dy
197
      call gettimes(cdfid,ml_time,ml_ntimes,ierr) 
198
      call getstart(cstid,ml_idate,ierr)
199
      call getpole(cstid,ml_pollon,ml_pollat,ierr)
200
      call clscdf(cstid,ierr)
201
      call clscdf(cdfid,ierr) 
202
 
203
c     Get grid description for Z file : pressure level
204
      call cdfopn(pl_zfn,cdfid,ierr)
205
      if (ierr.ne.0) goto 998
206
      call getcfn(cdfid,cfn,ierr)
207
      if (ierr.ne.0) goto 998
208
      call cdfopn(cfn,cstid,ierr)   
209
      if (ierr.ne.0) goto 998
210
      call getvars(cdfid,pl_nvars,pl_vnam,ierr)
211
      varname='Z'
212
      isok=0
213
      call check_varok (isok,varname,pl_vnam,pl_nvars)
214
      if (isok.ne.1) goto 998
215
      call getdef(cdfid,varname,pl_ndim,pl_mdv,pl_vardim,
216
     >            pl_varmin,pl_varmax,pl_stag,ierr)
217
      if (ierr.ne.0) goto 998
218
      pl_nx  =pl_vardim(1)
219
      pl_ny  =pl_vardim(2)
220
      pl_nz  =pl_vardim(3)
221
      pl_xmin=pl_varmin(1)
222
      pl_ymin=pl_varmin(2)
223
      call getlevs(cstid,pl_nz,pl_aklev,pl_bklev,pl_aklay,pl_bklay,ierr)
224
      call getgrid(cstid,pl_dx,pl_dy,ierr)
225
      pl_xmax=pl_xmin+real(pl_nx-1)*pl_dx
226
      pl_ymax=pl_ymin+real(pl_ny-1)*pl_dy
227
      call gettimes(cdfid,pl_time,pl_ntimes,ierr)
228
      call getstart(cstid,pl_idate,ierr)
229
      call getpole(cstid,pl_pollon,pl_pollat,ierr)
230
      call clscdf(cstid,ierr)
231
      call clscdf(cdfid,ierr) 
232
 
233
c     Set grid description for output file : height level
234
      zl_vardim(1) = pl_vardim(1)
235
      zl_vardim(2) = pl_vardim(2)
236
      zl_vardim(3) = nlev
237
      zl_varmin(1) = pl_varmin(1)
238
      zl_varmin(2) = pl_varmin(2)
239
      zl_varmin(3) = zlev(1)
240
      zl_varmax(1) = pl_varmax(1)
241
      zl_varmax(2) = pl_varmax(2)
242
      zl_varmax(3) = zlev(nlev)
243
      do i=1,nlev
244
         zl_aklay(i) = zlev(i)        
245
         zl_bklay(i) = 0.
246
         zl_aklev(i) = zlev(i)
247
         zl_bklev(i) = 0.
248
      enddo
249
      zl_dx       = pl_dx
250
      zl_dy       = pl_dy
251
      zl_time     = pl_time
252
      zl_ntimes   = pl_ntimes
253
      zl_ndim     = 4
254
      zl_mdv      = pl_mdv
255
      zl_nx       = zl_vardim(1)
256
      zl_ny       = zl_vardim(2)
257
      zl_nz       = zl_vardim(3)
258
      zl_xmin     = zl_varmin(1)
259
      zl_ymin     = zl_varmin(2)
260
      zl_pollon   = ml_pollon
261
      zl_pollat   = ml_pollat
262
      do i=1,5
263
         zl_idate(i) = ml_idate(i)
264
      enddo
265
 
266
c     Consitency check for the grids
267
      if ( (ml_nx.ne.pl_nx).or.
268
     >     (ml_ny.ne.pl_ny).or.
269
     >     (abs(ml_xmin-pl_xmin    ).gt.eps).or.
270
     >     (abs(ml_ymin-pl_ymin    ).gt.eps).or.
271
     >     (abs(ml_xmax-pl_xmax    ).gt.eps).or.
272
     >     (abs(ml_ymax-pl_ymax    ).gt.eps).or.
273
     >     (abs(ml_dx  -pl_dx      ).gt.eps).or.
274
     >     (abs(ml_dy  -pl_dy      ).gt.eps).or.
275
     >     (abs(ml_time-pl_time    ).gt.eps).or.
276
     >     (abs(ml_pollon-pl_pollon).gt.eps).or.
277
     >     (abs(ml_pollat-pl_pollat).gt.eps)) then
278
         print*,'Input P and Z grids are not consistent... Stop'
279
         print*
280
         print*,'Xmin   : ',ml_xmin,pl_xmin
281
         print*,'Ymin   : ',ml_ymin,pl_ymin
282
         print*,'Xmax   : ',ml_xmax,pl_xmax
283
         print*,'Ymax   : ',ml_ymax,pl_ymax
284
         print*,'Dx     : ',ml_dx,  pl_dx
285
         print*,'Dy     : ',ml_dy,  pl_dy
286
         print*,'Time   : ',ml_time,pl_time
287
         print*,'Pollon : ',ml_pollon,pl_pollon
288
         print*,'Pollat : ',ml_pollat,pl_pollat
289
         stop
290
      endif
291
 
292
c     Allocate memory for all fields
293
      allocate(ml_ps(ml_nx,ml_ny),stat=stat)
294
      if (stat.ne.0) print*,'*** error allocating array ml_ps ***'
295
      allocate(ml_zb(ml_nx,ml_ny),stat=stat)
296
      if (stat.ne.0) print*,'*** error allocating array ml_zb ***'
297
      allocate(ml_p3(ml_nx,ml_ny,ml_nz),stat=stat)
298
      if (stat.ne.0) print*,'*** error allocating array ml_p3 ***'
299
      allocate(ml_t3(ml_nx,ml_ny,ml_nz),stat=stat)
300
      if (stat.ne.0) print*,'*** error allocating array ml_t3 ***'
301
      allocate(ml_q3(ml_nx,ml_ny,ml_nz),stat=stat)
302
      if (stat.ne.0) print*,'*** error allocating array ml_q3 ***'
303
      allocate(ml_tv3(ml_nx,ml_ny,ml_nz),stat=stat)
304
      if (stat.ne.0) print*,'*** error allocating array ml_tv3 ***'
305
      allocate(ml_z3(ml_nx,ml_ny,ml_nz),stat=stat)
306
      if (stat.ne.0) print*,'*** error allocating array ml_z3 ***'
307
 
308
      allocate(in_field(ml_nx,ml_ny,ml_nz),stat=stat)
309
      if (stat.ne.0) print*,'*** error allocating array in_field ***'
310
 
311
      allocate(pl_z3(pl_nx,pl_ny,pl_nz),stat=stat)
312
      if (stat.ne.0) print*,'*** error allocating array pl_z3 ***'
313
      allocate(pl_p3(pl_nx,pl_ny,pl_nz),stat=stat)
314
      if (stat.ne.0) print*,'*** error allocating array pl_p3 ***'
315
 
316
      allocate(zl_field(zl_nx,zl_ny,zl_nz),stat=stat)
317
      if (stat.ne.0) print*,'*** error allocating array zl_field ***'
318
      allocate(zl_p(zl_nx,zl_ny,zl_nz),stat=stat)
319
      if (stat.ne.0) print*,'*** error allocating array zl_p ***'
320
 
321
c     Read T, Q, PS from P file
322
      call cdfopn(ml_pfn,cdfid,ierr)
323
      if (ierr.ne.0) goto 998
324
      isok=0
325
      varname='T'
326
      call check_varok (isok,varname, ml_vnam,ml_nvars)
327
      varname='Q'
328
      call check_varok (isok,varname, ml_vnam,ml_nvars)
329
      varname='PS'
330
      call check_varok (isok,varname,ml_vnam,ml_nvars)
331
      if (isok.ne.3) goto 998
332
      print*,'R T  ',trim(ml_pfn)
333
      call getdat(cdfid,'T',ml_time,0,ml_t3,ierr)
334
      print*,'R Q  ',trim(ml_pfn)
335
      call getdat(cdfid,'Q',ml_time,0,ml_q3,ierr)
336
      print*,'R PS ',trim(ml_pfn)
337
      call getdat(cdfid,'PS',ml_time,1,ml_ps,ierr)
338
      call clscdf(cdfid,ierr)
339
 
340
c     Read Z from Z file
341
      call cdfopn(pl_zfn,cdfid,ierr)
342
      if (ierr.ne.0) goto 998
343
      isok=0
344
      varname='Z'
345
      call check_varok (isok,varname,pl_vnam,pl_nvars)
346
      if (isok.ne.1) goto 998
347
      print*,'R Z  ',trim(pl_zfn)
348
      call getdat(cdfid,varname,pl_time,0,pl_z3,ierr)
349
      call clscdf(cdfid,ierr)
350
 
351
c     Set the values for the pressure on the pressure levels
352
      do i=1,pl_nx
353
         do j=1,pl_ny
354
            do k=1,pl_nz
355
               pl_p3(i,j,k)=pl_aklay(k)
356
            enddo
357
         enddo
358
      enddo
359
 
360
c     -----------------------------------------------------------------
361
c     Calculate geopotential on layers
362
c     -----------------------------------------------------------------
363
 
364
c     Calculate 3d pressure field
365
      print*,'C P'
366
      do k=1,ml_nz
367
        do i=1,ml_nx
368
           do j=1,ml_ny
369
              ml_p3(i,j,k)=ml_aklay(k)+ml_bklay(k)*ml_ps(i,j)
370
            enddo
371
         enddo
372
      enddo
373
 
374
c     Calculate 3d virtual temperature
375
      print*,'C TV'
376
      do k=1,ml_nz
377
         do i=1,ml_nx
378
            do j=1,ml_ny
379
               ml_tv3(i,j,k) = (ml_t3(i,j,k)+tzero)*
380
     >                         (1.+kappa*ml_q3(i,j,k))
381
            enddo
382
         enddo
383
      enddo
384
 
385
c     Loop over all grid points
386
      print*,'C HYDROSTATIC EQUATION'
387
      do i=1,ml_nx
388
         do j=1,ml_ny
389
 
390
c           Make the virtual temperature profile available 
391
            do k=1,ml_nz
392
               p1 (ml_nz-k+1)=ml_p3 (i,j,k)
393
               tv1(ml_nz-k+1)=ml_tv3(i,j,k)
394
            enddo
395
            call spline (p1,tv1,ml_nz,1.e30,1.e30,spline_tv1)
396
 
397
c           Loop over all model levels
398
            do k=1,ml_nz
399
 
400
c              Get pressure at the grid point
401
               p = ml_p3(i,j,k)
402
 
403
c              Find nearest pressure level which is above topography
404
               lmin=pl_nz
405
               do l=1,pl_nz
406
                  if ((abs(p-pl_p3(i,j,l))).lt.(abs(p-pl_p3(i,j,lmin)))
407
     >                .and.
408
     >               (pl_p3(i,j,l).lt.ml_ps(i,j)) ) then
409
                     lmin=l
410
                  endif
411
               enddo
412
 
413
c              Integrate hydrostatic equation from this level to the grid point
414
               p0 = pl_p3(i,j,lmin)
415
               n  = nint(abs(p-p0)/dpmin)
416
               if (n.lt.1) n=1
417
               dp = (p-p0)/real(n)
418
 
419
               pu  = p0
420
               z   = pl_z3(i,j,lmin)
421
               call splint(p1,tv1,spline_tv1,ml_nz,pu,tvu)
422
               do l=1,n
423
                  po  = pu+dp
424
                  call splint(p1,tv1,spline_tv1,ml_nz,po,tvo)
425
                  z   = z + rdg*0.5*(tvu+tvo)*alog(pu/po)
426
                  tvu = tvo
427
                  pu  = po
428
               enddo
429
 
430
c              Set the geopotential at the grid point
431
               ml_z3(i,j,k) = z
432
 
433
            enddo
434
 
435
         enddo
436
      enddo
437
 
438
c     -----------------------------------------------------------------
439
c     Calculate height of topography
440
c     -----------------------------------------------------------------
441
 
442
      print*,'C TOPOGRAPHY'
443
      do i=1,ml_nx
444
         do j=1,ml_ny
445
 
446
c           Make the z/p profile available 
447
            do k=1,ml_nz
448
               p1(ml_nz-k+1)=ml_p3(i,j,k)
449
               z1(ml_nz-k+1)=ml_z3(i,j,k)
450
            enddo
451
 
452
c           Cubic spline interpolation
453
            call spline (p1,z1,ml_nz,1.e30,1.e30,spline_z1) 
454
            call splint (p1,z1,spline_z1,ml_nz,ml_ps(i,j),ml_zb(i,j))
455
 
456
         enddo
457
      enddo
458
 
459
c     -----------------------------------------------------------------
460
c     Interpolate pressure onto z levels
461
c     -----------------------------------------------------------------
462
 
463
      do i=1,zl_nx
464
         do j=1,zl_ny
465
 
466
c           Make 1d profile available 
467
            do k=1,ml_nz
468
               z1(k)=ml_z3(i,j,k)
469
               f1(k)=ml_p3(i,j,k)
470
            enddo
471
            call spline (z1,f1,ml_nz,1.e30,1.e30,spline_f1)
472
            do k=1,zl_nz
473
               if (zl_aklay(k).gt.ml_zb(i,j)) then
474
                  call splint(z1,f1,spline_f1,ml_nz,zl_aklay(k),ff)
475
                  zl_p(i,j,k)=ff
476
               else
477
                  zl_p(i,j,k)=zl_mdv
478
               endif
479
            enddo
480
 
481
         enddo
482
      enddo    
483
 
484
c     -----------------------------------------------------------------
485
c     Write output to netcdf file
486
c     -----------------------------------------------------------------
487
 
488
      if (test.eq.1) then
489
 
490
c        Create the output file (grid taken from temperature)
491
         cdfname=testfile
492
 
493
         call cdfopn(ml_pfn,cdfid,ierr)
494
         if (ierr.ne.0) goto 998
495
         call getcfn(cdfid,cfn,ierr)
496
         if (ierr.ne.0) goto 998
497
         call clscdf(cdfid,ierr)
498
 
499
         call crecdf(cdfname,cdfid,
500
     >               ml_varmin,ml_varmax,ml_ndim,cfn,ierr)
501
         if (ierr.ne.0) goto 997
502
 
503
c        Write the definitions of the output variables
504
         varname='PS'
505
         ml_vardim(3)=1
506
         call putdef(cdfid,varname,ml_ndim,ml_mdv,
507
     >               ml_vardim,ml_varmin,ml_varmax,ml_stag,ierr)
508
         if (ierr.ne.0) goto 997
509
         varname='ORO'
510
         call putdef(cdfid,varname,ml_ndim,ml_mdv,
511
     >               ml_vardim,ml_varmin,ml_varmax,ml_stag,ierr)
512
         if (ierr.ne.0) goto 997
513
         ml_vardim(3)=ml_nz
514
         varname='Z'
515
         call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
516
     >               ml_varmin,ml_varmax,ml_stag,ierr)
517
         if (ierr.ne.0) goto 997
518
 
519
c        Write output fields
520
         print*,'W PS TEST'
521
         varname='PS'
522
         call putdat(cdfid,varname,ml_time,1,ml_ps,ierr)
523
         if (ierr.ne.0) goto 997
524
         print*,'W ZB TEST'
525
         varname='ORO'
526
         call putdat(cdfid,varname,ml_time,1,ml_zb,ierr)
527
         if (ierr.ne.0) goto 997
528
         print*,'W Z  TEST'
529
         varname='Z'
530
         call putdat(cdfid,varname,ml_time,0,ml_z3,ierr)
531
         if (ierr.ne.0) goto 997
532
 
533
c        Close cdf file
534
         call clscdf(cdfid,ierr)
535
 
536
      endif
537
 
538
c     -----------------------------------------------------------------
539
c     Interpolate fields onto a stack of z levels and write new file
540
c     -----------------------------------------------------------------
541
 
542
c     Create output file
543
      cfn=trim(zl_ofn)//'_cst'
544
      call crecdf(trim(zl_ofn),cdfid,zl_varmin,zl_varmax,
545
     >            zl_ndim,trim(cfn),ierr)
546
      if (ierr.ne.0) goto 995
547
 
548
c     Write topography
549
      print*,'W ZB ',trim(zl_ofn)
550
      varname='ORO'
551
      zl_vardim(3)=1
552
      call putdef(cdfid,varname,zl_ndim,zl_mdv,zl_vardim,
553
     >            zl_varmin,zl_varmax,zl_stag,ierr)         
554
      zl_vardim(3)=zl_nz
555
      if (ierr.ne.0) goto 995
556
      call putdat(cdfid,varname,zl_time,1,ml_zb,ierr)
557
      if (ierr.ne.0) goto 995
558
 
559
c     Write pressure
560
      print*,'W  P  ',trim(zl_ofn)
561
      varname='P'
562
      call putdef(cdfid,varname,4,zl_mdv,zl_vardim,
563
     >            zl_varmin,zl_varmax,zl_stag,ierr)         
564
      if (ierr.ne.0) goto 995
565
      call putdat(cdfid,varname,zl_time,0,zl_p,ierr)
566
      if (ierr.ne.0) goto 995
567
 
568
c     Close output file
569
      call clscdf(cdfid,ierr)
570
 
571
c     Loop over all variables
572
      do l=1,nvar
573
 
574
         print*,'I ',trim(varinp(l))
575
 
576
c        Get grid description for variable 
577
         call cdfopn(varsrc(l),cdfid,ierr)
578
         if (ierr.ne.0) goto 996
579
         call getcfn(cdfid,cfn,ierr)
580
         if (ierr.ne.0) goto 996
581
         call cdfopn(cfn,cstid,ierr)   
582
         if (ierr.ne.0) goto 996
583
         call getvars(cdfid,in_nvars,in_vnam,ierr)
584
         isok=0
585
         varname=varinp(l)
586
         call check_varok (isok,varname, in_vnam,in_nvars)
587
         if (isok.ne.1) goto 996
588
         call getdef(cdfid,varinp(l),in_ndim,in_mdv,in_vardim,
589
     >               in_varmin,in_varmax,in_stag,ierr)
590
         if (ierr.ne.0) goto 996
591
         in_nx  =in_vardim(1)
592
         in_ny  =in_vardim(2)
593
         in_nz  =in_vardim(3)
594
         in_xmin=in_varmin(1)
595
         in_ymin=in_varmin(2)
596
         call getlevs(cstid,in_nz,in_aklev,in_bklev,
597
     >                in_aklay,in_bklay,ierr)
598
         call getgrid(cstid,in_dx,in_dy,ierr)
599
         in_xmax=in_xmin+real(in_nx-1)*in_dx
600
         in_ymax=in_ymin+real(in_ny-1)*in_dy
601
         call gettimes(cdfid,in_time,in_ntimes,ierr) 
602
         call getstart(cstid,in_idate,ierr)
603
         call getpole(cstid,in_pollon,in_pollat,ierr)
604
         call clscdf(cstid,ierr)
605
         call clscdf(cdfid,ierr) 
606
 
607
c        Check whether this grid is consistent with P grid 
608
         if ( (ml_nx.ne.in_nx).or.
609
     >        (ml_ny.ne.in_ny).or.
610
     >        (ml_nz.ne.in_nz).or.
611
     >        (abs(ml_xmin-in_xmin      ).gt.eps).or.
612
     >        (abs(ml_ymin-in_ymin      ).gt.eps).or.
613
     >        (abs(ml_xmax-in_xmax      ).gt.eps).or.
614
     >        (abs(ml_ymax-in_ymax      ).gt.eps).or.
615
     >        (abs(ml_dx  -in_dx        ).gt.eps).or.
616
     >        (abs(ml_dy  -in_dy        ).gt.eps).or.
617
     >        (abs(ml_time-in_time      ).gt.eps).or.
618
     >        (abs(ml_stag(1)-in_stag(1)).gt.eps).or.
619
     >        (abs(ml_stag(2)-in_stag(2)).gt.eps).or.
620
     >        (abs(ml_stag(3)-in_stag(3)).gt.eps).or.
621
     >        (abs(ml_pollon-in_pollon  ).gt.eps).or.
622
     >        (abs(ml_pollat-in_pollat  ).gt.eps)) then
623
            print*,trim(varinp(l)),
624
     >            ' :Input P and FIELD grids are not consistent... Stop'
625
            print*,'Stag: ',ml_stag, in_stag
626
            print*,'Pol:  ',ml_pollon,in_pollon,ml_pollat,in_pollat
627
            stop
628
         endif
629
 
630
c        Read variable from file
631
         call cdfopn(varsrc(l),cdfid,ierr)
632
         if (ierr.ne.0) goto 996
633
         call getdat(cdfid,varinp(l),in_time,0,in_field,ierr)
634
         if (ierr.ne.0) goto 996
635
         call clscdf(cdfid,ierr) 
636
 
637
c        Write the constants file
638
         if (l.eq.1) then
639
 
640
            datar(1)=zl_nx
641
            datar(2)=zl_ny  
642
            datar(3)=int(1000.*zl_varmax(2))
643
            datar(4)=int(1000.*zl_varmin(1))
644
            datar(5)=int(1000.*zl_varmin(2))
645
            datar(6)=int(1000.*zl_varmax(1))
646
            datar(7)=int(1000.*zl_dx)
647
            datar(8)=int(1000.*zl_dy)
648
            datar(9)=zl_nz
649
            datar(10)=1
650
            datar(11)=1      
651
            datar(12)=0      
652
            datar(13)=int(1000.*zl_pollon) 
653
            datar(14)=int(1000.*zl_pollat) 
654
 
655
            cfn=trim(zl_ofn)//'_cst'
656
            call wricst(cfn,datar,zl_aklev,zl_bklev,
657
     >                  zl_aklay,zl_bklay,zl_idate)
658
 
659
         endif
660
 
661
c        Do the interpolation
662
         do i=1,zl_nx
663
            do j=1,zl_ny
664
 
665
c              Make 1d profile available 
666
               do k=1,ml_nz
667
                  z1(k)=ml_z3   (i,j,k)
668
                  f1(k)=in_field(i,j,k)
669
               enddo
670
               call spline (z1,f1,ml_nz,1.e30,1.e30,spline_f1)
671
               do k=1,zl_nz
672
                  if (zl_aklay(k).gt.ml_zb(i,j)) then
673
                     call splint(z1,f1,spline_f1,ml_nz,zl_aklay(k),ff)
674
                     zl_field(i,j,k)=ff
675
                  else
676
                     zl_field(i,j,k)=zl_mdv
677
                  endif
678
               enddo
679
 
680
            enddo
681
         enddo
682
 
683
c        Write the output field
684
         print*,'W ',trim(varout(l)),' ',trim(zl_ofn)
685
         call cdfwopn(trim(zl_ofn),cdfid,ierr)
686
         varname=trim(varout(l))
687
         call putdef(cdfid,varname,4,zl_mdv,zl_vardim,
688
     >               zl_varmin,zl_varmax,zl_stag,ierr)         
689
         if (ierr.ne.0) goto 995
690
         call putdat(cdfid,varname,zl_time,0,zl_field,ierr)
691
         if (ierr.ne.0) goto 995
692
         call clscdf(cdfid,ierr)
693
 
694
      enddo
695
 
696
c     -----------------------------------------------------------------
697
c     Write topography and geopotential also to the input P file
698
c     -----------------------------------------------------------------
699
 
700
c     Open the P file
701
      call cdfwopn(trim(ml_pfn),cdfid,ierr)     
702
 
703
c     Write topography
704
      varname='ORO'
705
      print*,'W ',trim(varname),' ',trim(ml_pfn)
706
      isok=0
707
      call check_varok (isok,varname,ml_vnam,ml_nvars)
708
      if (isok.eq.0) then
709
         ml_vardim(3)=1
710
         call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
711
     >               ml_varmin,ml_varmax,ml_stag,ierr)         
712
         ml_vardim(3)=ml_nz
713
         if (ierr.ne.0) goto 997
714
      endif
715
      call putdat(cdfid,varname,ml_time,1,ml_zb,ierr)
716
      if (ierr.ne.0) goto 997
717
 
718
c     Write geopotential height
719
      varname='Z'
720
      print*,'W ',trim(varname),' ',trim(ml_pfn)
721
      isok=0
722
      call check_varok (isok,varname,ml_vnam,ml_nvars)
723
      if (isok.eq.0) then
724
         call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
725
     >               ml_varmin,ml_varmax,ml_stag,ierr)         
726
         if (ierr.ne.0) goto 997
727
      endif
728
      call putdat(cdfid,varname,ml_time,0,ml_z3,ierr)
729
      if (ierr.ne.0) goto 997
730
 
731
c     Close P file
732
      call clscdf(cdfid,ierr)
733
 
734
c     -----------------------------------------------------------------
735
c     Exception handling and format specs
736
c     -----------------------------------------------------------------
737
 
738
      stop
739
 
740
 998  print*,'Z: Problems with input from m level'
741
      stop
742
 
743
 997  print*,'Z: Problems with output on m level'
744
      stop
745
 
746
 996  print*,'F: Problems with input from m level'
747
      stop
748
 
749
 995  print*,'F: Problems with output on z level'
750
      stop
751
 
752
 
753
      end
754
 
755
c     -------------------------------------------------------------
756
c     Natural cubic spline
757
c     -------------------------------------------------------------
758
 
759
      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
760
      INTEGER n,NMAX
761
      REAL yp1,ypn,x(n),y(n),y2(n)
762
      PARAMETER (NMAX=500)
763
      INTEGER i,k
764
      REAL p,qn,sig,un,u(NMAX)
765
      if (yp1.gt..99e30) then
766
        y2(1)=0.
767
        u(1)=0.
768
      else
769
        y2(1)=-0.5
770
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
771
      endif
772
      do 11 i=2,n-1
773
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
774
        p=sig*y2(i-1)+2.
775
        y2(i)=(sig-1.)/p
776
        u(i)=(6.*((y(i+1)-y(i))/(x(i+
777
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
778
     *u(i-1))/p
779
11    continue
780
      if (ypn.gt..99e30) then
781
        qn=0.
782
        un=0.
783
      else
784
        qn=0.5
785
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
786
      endif
787
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
788
      do 12 k=n-1,1,-1
789
        y2(k)=y2(k)*y2(k+1)+u(k)
790
12    continue
791
      return
792
      END
793
 
794
 
795
      SUBROUTINE splint(xa,ya,y2a,n,x,y)
796
      INTEGER n
797
      REAL x,y,xa(n),y2a(n),ya(n)
798
      INTEGER k,khi,klo
799
      REAL a,b,h
800
      klo=1
801
      khi=n
802
1     if (khi-klo.gt.1) then
803
        k=(khi+klo)/2
804
        if(xa(k).gt.x)then
805
          khi=k
806
        else
807
          klo=k
808
        endif
809
      goto 1
810
      endif
811
      h=xa(khi)-xa(klo)
812
      if (h.eq.0.) pause 'bad xa input in splint'
813
      a=(xa(khi)-x)/h
814
      b=(x-xa(klo))/h
815
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
816
     *2)/6.
817
      return
818
      END
819
 
820
c     ----------------------------------------------------------------
821
c     Check whether variable is found on netcdf file
822
c     ----------------------------------------------------------------
823
 
824
      subroutine check_varok (isok,varname,varlist,nvars)
825
 
826
c     Check whether the variable <varname> is in the list <varlist(nvars)>. 
827
c     If this is the case, <isok> is incremented by 1. Otherwise <isok> 
828
c     keeps its value.
829
 
830
      implicit none
831
 
832
c     Declaraion of subroutine parameters
833
      integer      isok
834
      integer      nvars
835
      character*80 varname
836
      character*80 varlist(nvars)
837
 
838
c     Auxiliary variables
839
      integer      i
840
 
841
c     Main
842
      do i=1,nvars
843
         if (trim(varname).eq.trim(varlist(i))) isok=isok+1
844
      enddo
845
 
846
      end