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 modify_anomaly
2
 
3
c     ******************************************************************************
4
c     * Read the modified and unmodified PV from the R file and                  *
5
c     * perform some non-standard modifications: (a) Change of amplitude,          *
6
c     * (b) Stretching and shrinking along an axis, (c) rotation, or               *
7
c     * (d) change in position.     
8
c     *                                                                            
9
c     # modify_anomaly.sh R_20060116_18 shifty=-10,rot=10,stretchy=1.5,rot=-10,shifty=10
10
c     # modify_anomaly.sh R_20060116_18 fac=2
11
c     # modify_anomaly.sh R_20060116_18 cex=-30,cey=-30,rot=20
12
c     *                                                                            *
13
c     * Michael Sprenger / Spring, Summer 2007                                     *
14
c     ******************************************************************************
15
 
16
      implicit none
17
 
18
c     -----------------------------------------------------------------------------
19
c     Declaration of parameters and variables
20
c     -----------------------------------------------------------------------------
21
 
22
c     Input/output file and command string
23
      character*80   pvsrcfile
24
      character*80   commandstr
25
 
26
c     Grid parameters
27
      integer        nx,ny,nz
28
      real           xmin,ymin,zmin,xmax,ymax,zmax
29
      real           dx,dy,dz
30
      real           mdv
31
 
32
c     3d fields for calculation of qgPV and Ertel's PV
33
      real,allocatable,dimension (:,:,:) :: pv1,pv2,ano
34
 
35
c     Numerical epsilon
36
      real           eps
37
      parameter      (eps=0.01)
38
 
39
c     Parameters for the transformations
40
      real           cex,cey              ! Centre point for rotation
41
      real           angle                ! Rotation angle
42
 
43
c     Auxiliary variables
44
      real           zpos
45
      integer        i,j,k
46
      integer        stat
47
      character*80   varname
48
      integer        n
49
      real           par
50
      character*80   com
51
 
52
c     -----------------------------------------------------------------------------
53
c     Preparations
54
c     -----------------------------------------------------------------------------
55
 
56
      print*,'********************************************************'
57
      print*,'* MODIFY_ANOMALY                                       *'
58
      print*,'********************************************************'
59
 
60
c     Read parameter file
61
      open(10,file='fort.10')
62
       read(10,*) pvsrcfile
63
       read(10,*) commandstr
64
      close(10)
65
      print*
66
      print*,'Input file : ',trim(pvsrcfile)
67
      print*,'Command    : ',trim(commandstr)
68
      print*
69
 
70
c     Get lat/lon gid parameters from input file
71
      call read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv,
72
     >               pvsrcfile)
73
      print*,'Read_Dim: nx,ny,nz         = ',nx,ny,nz
74
      print*,'          dx,dy,dz         = ',dx,dy,dz
75
      print*,'          xmin,ymin,zmin   = ',xmin,ymin,zmin
76
      print*,'          mdv              = ',mdv
77
      print*
78
 
79
c     Count from 0, not from 1: consistent with <inv_cart.f>.
80
      nx=nx-1
81
      ny=ny-1
82
      nz=nz-1
83
 
84
c     Allocate memory for modified and non-modified PV
85
      allocate(pv1 (0:nx,0:ny,0:nz),STAT=stat)
86
      if (stat.ne.0) print*,'error allocating pv1'
87
      allocate(pv2 (0:nx,0:ny,0:nz),STAT=stat)
88
      if (stat.ne.0) print*,'error allocating pv2'
89
      allocate(ano  (0:nx,0:ny,0:nz),STAT=stat)
90
      if (stat.ne.0) print*,'error allocating ano'   
91
 
92
c     Read data from file
93
      varname='PV'
94
      call read_inp (pv1,varname,pvsrcfile,
95
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
96
      varname='PV_FILT'
97
      call read_inp (pv2,varname,pvsrcfile,
98
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
99
 
100
c     Define the anomaly
101
      do i=0,nx
102
         do j=0,ny
103
            do k=0,nz
104
               if ( (abs(pv1(i,j,k)-mdv).gt.eps).and.
105
     >              (abs(pv2(i,j,k)-mdv).gt.eps) ) then
106
                  ano(i,j,k)=pv1(i,j,k)-pv2(i,j,k)
107
               else
108
                  ano(i,j,k)=0.
109
               endif
110
            enddo
111
         enddo
112
      enddo
113
 
114
 
115
c     -------------------------------------------------------------------------------
116
c     Modifications
117
c     -------------------------------------------------------------------------------
118
 
119
c     Set the default values for parameters
120
      cex   = 0.
121
      cey   = 0.
122
      angle = 0.
123
 
124
c     Set the counter for the command string
125
      n=1
126
 
127
c     Loop over all commands
128
 100  continue
129
 
130
c     Extract new command/parameter pair; exit if no new command
131
      call next_command(commandstr,n,com,par)
132
      if (com.eq.'nil') goto 200
133
      print*,trim(com),par
134
 
135
c     Multiply the anomaly by a constant factor
136
      if (com.eq.'fac') then
137
         call mod_factor (ano,par,nx,ny,nz,mdv)
138
      endif
139
 
140
c     Set the centre point (needed for rotations)
141
      if (com.eq.'cex') then
142
         cex=par
143
      endif
144
      if (com.eq.'cey') then
145
         cey=par
146
      endif
147
 
148
c     Rotation around the centrre point
149
      if (com.eq.'rot') then
150
         call mod_rotation (ano,par,cex,cey,
151
     >                       nx,ny,nz,xmin,ymin,dx,dy,mdv)
152
      endif
153
 
154
c     Shift in x or in y direction
155
      if (com.eq.'shiftx') then
156
         call  mod_shift (ano,par,0.,
157
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
158
      endif
159
      if (com.eq.'shifty') then
160
         call  mod_shift (ano,0.,par,
161
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
162
      endif      
163
 
164
c     Stretch/shrink in x or in y direction
165
      if (com.eq.'stretchx') then
166
         call  mod_stretch (ano,par,1.,
167
     >                      nx,ny,nz,xmin,ymin,dx,dy,mdv)
168
      endif
169
      if (com.eq.'stretchy') then
170
         call  mod_stretch (ano,1.,par,
171
     >                      nx,ny,nz,xmin,ymin,dx,dy,mdv)
172
      endif      
173
 
174
 
175
c     Goto next command/parameter pair
176
      goto 100
177
 
178
c     All commands handled
179
 200  continue
180
 
181
c     --------------------------------------------------------------------------------
182
c     Write modified fields
183
c     --------------------------------------------------------------------------------
184
 
185
c     Build the modified PV distribution from the anomaly
186
      do i=0,nx
187
         do j=0,ny
188
            do k=0,nz
189
               if ( (abs(pv1(i,j,k)-mdv).gt.eps).and.
190
     >              (abs(ano(i,j,k)-mdv).gt.eps) ) then
191
                  pv2(i,j,k)=pv1(i,j,k)-ano(i,j,k)
192
               else
193
                  pv2(i,j,k)=mdv
194
               endif
195
            enddo
196
         enddo
197
      enddo
198
 
199
c     Write result to netcdf file
200
      varname='PV_FILT'
201
      call write_inp (pv2,varname,pvsrcfile,nx,ny,nz)
202
 
203
c     Write the modified anomaly also to the file
204
      varname='PV_ANOM'
205
      call write_inp (ano,varname,pvsrcfile,nx,ny,nz)
206
 
207
      end
208
 
209
c     ********************************************************************************
210
c     * Command handling and transformations                                         *
211
c     ********************************************************************************
212
 
213
c     --------------------------------------------------------------------------------
214
c     Extract next command from command string
215
c     --------------------------------------------------------------------------------
216
 
217
      subroutine next_command (commandstr,n,com,par)
218
 
219
c     Given the command string <commandstr>, extract the next command/parameter pair
220
c     <com,par> starting at position <n> of the command string.
221
 
222
      implicit none
223
 
224
c     Declaration of subroutine parameters
225
      character*80 commandstr
226
      character*80 com
227
      real         par
228
      integer      n
229
 
230
c     Auxiliary variables
231
      integer      i,j,k
232
 
233
c     Check whether end of command line reached
234
      if (n.ge.80) then
235
         com='nil'
236
         par=0.
237
         goto 120
238
      endif
239
 
240
c     Set indices to next next command/parameter pair
241
      i=n
242
      j=n
243
 100  if ( (commandstr(j:j).ne.'=').and.(j.lt.80) ) then
244
        j=j+1
245
        goto 100
246
      endif
247
      k=j+1
248
 110  if ( (commandstr(k:k).ne.',' ).and.
249
     >     (commandstr(k:k).ne.';' ).and.
250
     >     (commandstr(k:k).ne.' ' ).and.
251
     >     (k              .lt.80  ) ) then
252
        k=k+1
253
        goto 110
254
      endif
255
 
256
c     Check whether this is a valid command/parameter pair
257
      if ( ((j-1).lt.i).or.((k-1).lt.(j+1)) ) then
258
         com='nil'
259
         par=0.
260
         goto 120
261
      endif
262
 
263
c     Extract tzhe command and the parameter
264
      com=commandstr(i:j-1)
265
      read(commandstr(j+1:k-1),*) par
266
 
267
c     Set the counter to the next command/parameter pair
268
      n=k+1
269
 
270
c     Exit point
271
 120  continue
272
 
273
      end
274
 
275
c     --------------------------------------------------------------------------------
276
c     Multiply the anomaly by a factor
277
c     --------------------------------------------------------------------------------
278
 
279
      subroutine mod_factor (field,factor,nx,ny,nz,mdv)
280
 
281
c     Multiply the anomaly <field>  by a constant factor <factor>. The grid and the
282
c     missing data value are given by <nx,ny,nz,mdv>.
283
 
284
      implicit none
285
 
286
c     Declaration of subroutine parameters
287
      integer    nx,ny,nz
288
      real       field(0:nx,0:ny,0:nz)
289
      real       mdv
290
      real       factor
291
 
292
c     Parameters
293
      real       eps
294
      parameter  (eps=0.01)
295
 
296
c     Auxiliary variables
297
      integer    i,j,k
298
 
299
c     Do the transformation
300
      do i=0,nx
301
         do j=0,ny
302
            do k=0,nz
303
               if ( (abs(field(i,j,k)-mdv).gt.eps) ) then 
304
                  field(i,j,k)=factor*field(i,j,k)
305
               endif
306
            enddo
307
         enddo
308
      enddo
309
 
310
      end
311
 
312
c     --------------------------------------------------------------------------------
313
c     Rotate the anomaly
314
c     --------------------------------------------------------------------------------
315
 
316
      subroutine mod_rotation (field,angle,cex,cey,
317
     >                         nx,ny,nz,xmin,ymin,dx,dy,mdv)
318
 
319
c     Rotate the anomaly <field> in the horizontal by the angle <angle>. The centre
320
c     for the rotation is given by <cex,cey>, expressed in rotated longitude, latitude.
321
c     The grid and the missing data value are given by <nx,ny,nz,xmin,ymin,dx,dy,mdv>.
322
 
323
 
324
      implicit none
325
 
326
c     Declaration of subroutine parameters
327
      integer    nx,ny,nz
328
      real       field(0:nx,0:ny,0:nz)
329
      real       xmin,ymin,dx,dy
330
      real       mdv
331
      real       angle
332
      real       cex,cey
333
 
334
c     Parameters
335
      real       eps
336
      parameter  (eps=0.01)
337
      real       pi180
338
      parameter  (pi180=3.14159/180.)
339
 
340
c     Auxiliary variables
341
      integer    i,j,k
342
      real       lon,lat,rlon,rlat
343
      real       xmax,ymax
344
      real       tmp1(0:nx,0:ny),tmp2(0:nx,0:ny)
345
 
346
c     Externals 
347
      real       int2d
348
      external   int2d
349
 
350
c     Maximal grid extension
351
      xmax=xmin+real(nx-1)*dx
352
      ymax=ymin+real(ny-1)*dy
353
 
354
c     Do the Transformation
355
      do k=0,nz
356
 
357
c        Copy level to 2d array
358
         do i=0,nx
359
            do j=0,ny
360
               tmp1(i,j)=field(i,j,k)
361
            enddo
362
         enddo
363
 
364
c        Rotate each grid point
365
         do i=0,nx
366
            do j=0,ny
367
 
368
c              Lon/lat coordinates of grid point
369
               lon=xmin+real(i)*dx-cex
370
               lat=ymin+real(j)*dy-cey
371
 
372
c              Rotation
373
               rlon = cex + lon*cos(angle*pi180) + lat*sin(angle*pi180)
374
               rlat = cey - lon*sin(angle*pi180) + lat*cos(angle*pi180)
375
 
376
c              Do the interpolation
377
               if ( (rlon.gt.xmin).and.(rlon.lt.xmax).and.
378
     >              (rlat.gt.ymin).and.(rlat.lt.ymax) ) 
379
     >         then
380
                  tmp2(i,j)=int2d(tmp1,rlon,rlat,
381
     >                            dx,dy,xmin,ymin,nx,ny,mdv)
382
               else
383
                  tmp2(i,j)=0.
384
               endif
385
 
386
            enddo
387
         enddo
388
 
389
c        Copy 2d array to level
390
         do i=0,nx
391
            do j=0,ny
392
               field(i,j,k)=tmp2(i,j)
393
            enddo
394
         enddo
395
 
396
      enddo
397
 
398
      end
399
 
400
c     --------------------------------------------------------------------------------
401
c     Shift the anomaly in x and y direction
402
c     --------------------------------------------------------------------------------
403
 
404
      subroutine mod_shift (field,shiftx,shifty,
405
     >                      nx,ny,nz,xmin,ymin,dx,dy,mdv)
406
 
407
c     Shift the anomaly <field> in the horizontal by the vector <shiftx,shifty>. 
408
c     The grid and the missing data value are given by <nx,ny,nz,xmin,ymin,dx,dy,mdv>.
409
 
410
 
411
      implicit none
412
 
413
c     Declaration of subroutine parameters
414
      integer    nx,ny,nz
415
      real       field(0:nx,0:ny,0:nz)
416
      real       xmin,ymin,dx,dy
417
      real       mdv
418
      real       shiftx,shifty
419
 
420
c     Parameters
421
      real       eps
422
      parameter  (eps=0.01)
423
 
424
c     Auxiliary variables
425
      integer    i,j,k
426
      real       lon,lat,rlon,rlat
427
      real       xmax,ymax
428
      real       tmp1(0:nx,0:ny),tmp2(0:nx,0:ny)
429
 
430
c     Externals 
431
      real       int2d
432
      external   int2d
433
 
434
c     Maximal grid extension
435
      xmax=xmin+real(nx-1)*dx
436
      ymax=ymin+real(ny-1)*dy
437
 
438
c     Do the Transformation
439
      do k=0,nz
440
 
441
c        Copy level to 2d array
442
         do i=0,nx
443
            do j=0,ny
444
               tmp1(i,j)=field(i,j,k)
445
            enddo
446
         enddo
447
 
448
c        Rotate each grid point
449
         do i=0,nx
450
            do j=0,ny
451
 
452
c              Lon/lat coordinates of grid point
453
               lon=xmin+real(i)*dx
454
               lat=ymin+real(j)*dy
455
 
456
c              shifted coordinates
457
               rlon = lon - shiftx
458
               rlat = lat - shifty
459
 
460
c              Do the interpolation
461
               if ( (rlon.gt.xmin).and.(rlon.lt.xmax).and.
462
     >              (rlat.gt.ymin).and.(rlat.lt.ymax) ) 
463
     >         then
464
                  tmp2(i,j)=int2d(tmp1,rlon,rlat,
465
     >                            dx,dy,xmin,ymin,nx,ny,mdv)
466
               else
467
                  tmp2(i,j)=0.
468
               endif
469
 
470
            enddo
471
         enddo
472
 
473
c        Copy 2d array to level
474
         do i=0,nx
475
            do j=0,ny
476
               field(i,j,k)=tmp2(i,j)
477
            enddo
478
         enddo
479
 
480
      enddo
481
 
482
      end
483
 
484
c     --------------------------------------------------------------------------------
485
c     Stretch/shrink the anomaly in x and y direction
486
c     --------------------------------------------------------------------------------
487
 
488
      subroutine mod_stretch (field,stretchx,stretchy,
489
     >                        nx,ny,nz,xmin,ymin,dx,dy,mdv)
490
 
491
c     Stretch the anomaly <field> in the horizontal by the fatcors <stretchx,stretchy>. 
492
c     The grid and the missing data value are given by <nx,ny,nz,xmin,ymin,dx,dy,mdv>.
493
 
494
 
495
      implicit none
496
 
497
c     Declaration of subroutine parameters
498
      integer    nx,ny,nz
499
      real       field(0:nx,0:ny,0:nz)
500
      real       xmin,ymin,dx,dy
501
      real       mdv
502
      real       stretchx,stretchy
503
 
504
c     Parameters
505
      real       eps
506
      parameter  (eps=0.01)
507
 
508
c     Auxiliary variables
509
      integer    i,j,k
510
      real       lon,lat,rlon,rlat
511
      real       xmax,ymax
512
      real       tmp1(0:nx,0:ny),tmp2(0:nx,0:ny)
513
 
514
c     Externals 
515
      real       int2d
516
      external   int2d
517
 
518
c     Maximal grid extension
519
      xmax=xmin+real(nx-1)*dx
520
      ymax=ymin+real(ny-1)*dy
521
 
522
c     Do the Transformation
523
      do k=0,nz
524
 
525
c        Copy level to 2d array
526
         do i=0,nx
527
            do j=0,ny
528
               tmp1(i,j)=field(i,j,k)
529
            enddo
530
         enddo
531
 
532
c        Rotate each grid point
533
         do i=0,nx
534
            do j=0,ny
535
 
536
c              Lon/lat coordinates of grid point
537
               lon=xmin+real(i)*dx
538
               lat=ymin+real(j)*dy
539
 
540
c              shifted coordinates
541
               rlon = 1./stretchx * lon
542
               rlat = 1./stretchy * lat
543
 
544
c              Do the interpolation
545
               if ( (rlon.gt.xmin).and.(rlon.lt.xmax).and.
546
     >              (rlat.gt.ymin).and.(rlat.lt.ymax) ) 
547
     >         then
548
                  tmp2(i,j)=int2d(tmp1,rlon,rlat,
549
     >                            dx,dy,xmin,ymin,nx,ny,mdv)
550
               else
551
                  tmp2(i,j)=0.
552
               endif
553
 
554
            enddo
555
         enddo
556
 
557
c        Copy 2d array to level
558
         do i=0,nx
559
            do j=0,ny
560
               field(i,j,k)=tmp2(i,j)
561
            enddo
562
         enddo
563
 
564
      enddo
565
 
566
      end
567
 
568
c     ------------------------------------------------------------------------
569
c     Two-dimensional interpolation         
570
c     ------------------------------------------------------------------------
571
 
572
      real function int2d(ar,lon,lat,
573
     >                    dx,dy,xmin,ymin,nx,ny,mdv)
574
 
575
c     Interpolate the field <ar(nx,ny)> to the position <lon,lat>. The
576
c     grid is specified by <dx,dy,xmin,ymin,nx,ny>.
577
 
578
      implicit none
579
 
580
c     Declaration of subroutine paramters
581
      integer   nx,ny
582
      real      ar(0:nx,0:ny)
583
      real      lon,lat
584
      real      dx,dy,xmin,ymin
585
      real      mdv
586
 
587
c     Parameters
588
      real      eps
589
      parameter (eps=0.01)
590
 
591
c     Auxiliary variables
592
      real      ri,rj
593
      integer   i,j,ir,ju     
594
      real      frac0i,frac0j,frac1i,frac1j
595
      real      tmp,nor
596
 
597
c     Get index
598
      ri=(lon-xmin)/dx
599
      rj=(lat-ymin)/dy
600
      i=int(ri)
601
      j=int(rj)
602
      if ((i.lt.0).or.(i.gt.nx).or.
603
     >    (j.lt.0).or.(j.gt.ny)) then
604
         print*,'lat/lon interpolation not possible....'
605
         stop
606
      endif
607
 
608
c     Get inidices of left and upper neighbours
609
      ir=i+1
610
      if (ir.gt.nx) ir=nx
611
      ju=j+1
612
      if (ju.gt.ny) ju=ny
613
 
614
c     Get the weights for the bilinear interpolation
615
      frac0i=ri-float(i)
616
      frac0j=rj-float(j)
617
      frac1i=1.-frac0i
618
      frac1j=1.-frac0j
619
 
620
c     Bilinear interpolation with missing data check
621
      tmp=0.
622
      nor=0.
623
      if ( (abs(ar(i ,j )-mdv).gt.eps) ) then
624
         tmp = tmp + ar(i ,j ) * frac1i * frac1j
625
         nor = nor + frac1i * frac1j
626
      endif
627
      if ( (abs(ar(i ,ju)-mdv).gt.eps) ) then
628
         tmp = tmp + ar(i ,ju) * frac1i * frac0j
629
         nor = nor + frac1i * frac0j
630
      endif
631
      if ( (abs(ar(ir,j )-mdv).gt.eps) ) then
632
         tmp = tmp + ar(ir,j ) * frac0i * frac1j
633
         nor = nor + frac0i * frac1j
634
      endif
635
      if ( (abs(ar(ir,ju)-mdv).gt.eps) ) then
636
         tmp = tmp + ar(ir,ju) * frac0i * frac0j
637
         nor = nor + frac0i * frac0j
638
      endif
639
 
640
c     Return result
641
      int2d=tmp/nor
642
 
643
      end
644
 
645
c     ********************************************************************************
646
c     * NETCDF INPUT AND OUTPUT                                                      *
647
c     ********************************************************************************
648
 
649
c     --------------------------------------------------------------------------------
650
c     Write input field to netcdf
651
c     --------------------------------------------------------------------------------
652
 
653
      SUBROUTINE write_inp (field,fieldname,pvsrcfile,nx,ny,nz)
654
 
655
c     Read <fieldname> from netcdf file <pvsrcfile> into <field>. The grid is specified 
656
c     by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed whether the input 
657
c     files are consitent with this grid. 
658
 
659
      implicit   none
660
 
661
c     Declaration of subroutine parameters
662
      integer              nx,ny,nz
663
      real                 field (0:nx,0:ny,0:nz)
664
      character*80         fieldname
665
      character*80         pvsrcfile
666
 
667
c     Auxiliary variables
668
      integer              cdfid,stat
669
      integer              vardim(4)
670
      real                 misdat
671
      real                 varmin(4),varmax(4),stag(4)
672
      integer              ndimin,outid,i,j,k
673
      real                 tmp(0:nx,0:ny,0:nz)
674
      integer              ntimes
675
      real                 time(200)       
676
      integer              nvars
677
      character*80         vnam(100),varname
678
      integer              isok
679
 
680
c     Get grid parameters from PV
681
      call cdfopn(pvsrcfile,cdfid,stat)
682
      if (stat.ne.0) goto 998
683
      call getvars(cdfid,nvars,vnam,stat)
684
      if (stat.ne.0) goto 998
685
      isok=0
686
      varname='PV'
687
      call check_varok(isok,varname,vnam,nvars)
688
      if (isok.eq.0) goto 998
689
      call getdef(cdfid,varname,ndimin,misdat,vardim,
690
     >            varmin,varmax,stag,stat)    
691
      if (stat.ne.0) goto 998
692
      time(1)=0.
693
      call gettimes(cdfid,time,ntimes,stat)  
694
      if (stat.ne.0) goto 998
695
      call clscdf(cdfid,stat)
696
      if (stat.ne.0) goto 998
697
 
698
c     Save variables (write definition, if necessary)
699
      call cdfwopn(pvsrcfile,cdfid,stat)
700
      if (stat.ne.0) goto 998
701
      isok=0
702
      varname=fieldname
703
      call check_varok(isok,varname,vnam,nvars)
704
      if (isok.eq.0) then
705
         call putdef(cdfid,varname,ndimin,misdat,vardim,
706
     >               varmin,varmax,stag,stat)
707
         if (stat.ne.0) goto 998
708
      endif
709
      call putdat(cdfid,varname,time(1),0,field,stat)
710
      print*,'W ',trim(varname),' ',trim(pvsrcfile)
711
      if (stat.ne.0) goto 998
712
 
713
c     Close input netcdf file
714
      call clscdf(cdfid,stat)
715
      if (stat.ne.0) goto 998
716
 
717
      return
718
 
719
c     Exception handling
720
 998  print*,'Write_Inp: Problem with input netcdf file... Stop'
721
      stop
722
 
723
      end
724
 
725
 
726
c     --------------------------------------------------------------------------------
727
c     Read input fields for reference profile
728
c     --------------------------------------------------------------------------------
729
 
730
      SUBROUTINE read_inp (field,fieldname,pvsrcfile,
731
     >                     nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
732
 
733
c     Read <fieldname> from netcdf file <pvsrcfile> into <field>. The grid is specified 
734
c     by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed whether the input 
735
c     files are consitent with this grid. The missing data value is set to <mdv>.
736
 
737
      implicit   none
738
 
739
c     Declaration of subroutine parameters
740
      integer              nx,ny,nz
741
      real                 field(0:nx,0:ny,0:nz)
742
      character*80         fieldname
743
      character*80         pvsrcfile
744
      real                 dx,dy,dz,xmin,ymin,zmin
745
      real                 mdv
746
 
747
c     Numerical and physical parameters
748
      real                 eps
749
      parameter            (eps=0.01)
750
 
751
c     Auxiliary variables
752
      integer              cdfid,stat,cdfid99
753
      integer              vardim(4)
754
      real                 misdat
755
      real                 varmin(4),varmax(4),stag(4)
756
      integer              ndimin,outid,i,j,k
757
      real                 tmp(nx,ny,nz)
758
      integer              ntimes
759
      real                 time(200)       
760
      integer              nvars
761
      character*80         vnam(100),varname
762
      integer              isok
763
 
764
c     Open the input netcdf file
765
      call cdfopn(pvsrcfile,cdfid,stat)
766
      if (stat.ne.0) goto 998
767
 
768
c     Check whether needed variables are on file
769
      call getvars(cdfid,nvars,vnam,stat)
770
      if (stat.ne.0) goto 998
771
      isok=0
772
      varname=trim(fieldname)
773
      call check_varok(isok,varname,vnam,nvars)
774
      if (isok.eq.0) goto 998
775
 
776
c     Get the grid parameters from theta     
777
      call getdef(cdfid,varname,ndimin,misdat,vardim,
778
     >            varmin,varmax,stag,stat)    
779
      if (stat.ne.0) goto 998
780
      time(1)=0.
781
      call gettimes(cdfid,time,ntimes,stat)  
782
      if (stat.ne.0) goto 998
783
 
784
c     Check whether grid parameters are consistent
785
      if ( (vardim(1).ne.(nx+1)).or.
786
     >     (vardim(2).ne.(ny+1)).or.
787
     >     (vardim(3).ne.(nz+1)).or.
788
     >     (abs(varmin(1)-xmin).gt.eps).or.
789
     >     (abs(varmin(2)-ymin).gt.eps).or.
790
     >     (abs(varmin(3)-zmin).gt.eps).or.
791
     >     (abs((varmax(1)-varmin(1))/real(vardim(1)-1)-dx).gt.eps).or.
792
     >     (abs((varmax(2)-varmin(2))/real(vardim(2)-1)-dy).gt.eps).or.
793
     >     (abs((varmax(3)-varmin(3))/real(vardim(3)-1)-dz).gt.eps) ) 
794
     >then
795
         print*,'Input grid inconsitency...'
796
         print*,'  Nx      = ',vardim(1),nx+1
797
         print*,'  Ny      = ',vardim(2),ny+1
798
         print*,'  Nz      = ',vardim(3),nz+1
799
         print*,'  Varminx = ',varmin(1),xmin
800
         print*,'  Varminy = ',varmin(2),ymin
801
         print*,'  Varminz = ',varmin(3),zmin
802
         print*,'  Dx      = ',(varmax(1)-varmin(1))/real(nx-1),dx
803
         print*,'  Dy      = ',(varmax(2)-varmin(2))/real(ny-1),dy
804
         print*,'  Dz      = ',(varmax(3)-varmin(3))/real(nz-1),dz
805
         goto 998
806
      endif
807
 
808
c     Load variables
809
      call getdef(cdfid,varname,ndimin,misdat,vardim,
810
     >            varmin,varmax,stag,stat)
811
      if (stat.ne.0) goto 998
812
      call getdat(cdfid,varname,time(1),0,field,stat)
813
      print*, 'R ',trim(varname),' ',trim(pvsrcfile)
814
      if (stat.ne.0) goto 998
815
 
816
c     Close input netcdf file
817
      call clscdf(cdfid,stat)
818
      if (stat.ne.0) goto 998
819
 
820
c     Set missing data value to <mdv>
821
      do i=1,nx
822
         do j=1,ny
823
            do k=1,nz
824
               if (abs(field(i,j,k)-misdat).lt.eps) then
825
                  field(i,j,k)=mdv
826
               endif
827
            enddo
828
         enddo
829
      enddo
830
 
831
      return
832
 
833
c     Exception handling
834
 998  print*,'Read_Inp: Problem with input netcdf file... Stop'
835
      stop
836
 
837
      end
838
 
839
c     --------------------------------------------------------------------------------
840
c     Check whether variable is found on netcdf file
841
c     --------------------------------------------------------------------------------
842
 
843
      subroutine check_varok (isok,varname,varlist,nvars)
844
 
845
c     Check whether the variable <varname> is in the list <varlist(nvars)>. If this is 
846
C     the case, <isok> is incremented by 1. Otherwise <isok> keeps its value.
847
 
848
      implicit none
849
 
850
c     Declaraion of subroutine parameters
851
      integer      isok
852
      integer      nvars
853
      character*80 varname
854
      character*80 varlist(nvars)
855
 
856
c     Auxiliary variables
857
      integer      i
858
 
859
c     Main
860
      do i=1,nvars
861
         if (trim(varname).eq.trim(varlist(i))) isok=isok+1
862
      enddo
863
 
864
      end
865
 
866
c     --------------------------------------------------------------------------------
867
c     Get grid parameters
868
c     --------------------------------------------------------------------------------
869
 
870
      subroutine read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv,
871
     >                     pvsrcfile)
872
 
873
c     Get the grid parameters from the variable <THETA> on the input file <pvsrcfile>.
874
c     The grid parameters are
875
c        nx,ny,nz                : Number of grid points in x, y and z direction
876
c        xmin,ymin,zmin          : Minimum domain coordinates in x, y and z direction
877
c        xmax,ymax,zmax          : Maximal domain coordinates in x, y and z direction
878
c        dx,dy,dz                : Horizontal and vertical resolution
879
c     Additionally, it is checked whether the vertical grid is equally spaced. If ok,
880
c     the grid paramters are transformed from lon/lat to distance (in meters)
881
 
882
      implicit none
883
 
884
c     Declaration of subroutine parameters
885
      character*80   pvsrcfile
886
      integer        nx,ny,nz
887
      real           dx,dy,dz
888
      real           xmin,ymin,zmin,xmax,ymax,zmax
889
      real           mdv
890
 
891
c     Numerical epsilon and other physical/geoemtrical parameters
892
      real           eps
893
      parameter      (eps=0.01)
894
 
895
c     Auxiliary variables
896
      integer        cdfid,cstid
897
      integer        ierr
898
      character*80   vnam(100),varname
899
      integer        nvars
900
      integer        isok
901
      integer        vardim(4)
902
      real           misdat
903
      real           varmin(4),varmax(4),stag(4)
904
      real           aklev(1000),bklev(1000),aklay(1000),bklay(1000)
905
      real           dh
906
      character*80   csn
907
      integer        ndim
908
      integer        i
909
 
910
c     Get all grid parameters
911
      call cdfopn(pvsrcfile,cdfid,ierr)
912
      if (ierr.ne.0) goto 998
913
      call getvars(cdfid,nvars,vnam,ierr)
914
      if (ierr.ne.0) goto 998
915
      isok=0
916
      varname='PV'
917
      call check_varok(isok,varname,vnam,nvars)
918
      if (isok.eq.0) goto 998
919
      call getcfn(cdfid,csn,ierr)
920
      if (ierr.ne.0) goto 998
921
      call cdfopn(csn,cstid,ierr)
922
      if (ierr.ne.0) goto 998
923
      call getdef(cdfid,varname,ndim,misdat,vardim,varmin,varmax,
924
     >            stag,ierr)
925
      if (ierr.ne.0) goto 998
926
      nx=vardim(1)
927
      ny=vardim(2)
928
      nz=vardim(3)
929
      xmin=varmin(1)
930
      ymin=varmin(2)
931
      zmin=varmin(3)
932
      call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)
933
      if (ierr.ne.0) goto 998
934
      call getgrid(cstid,dx,dy,ierr)
935
      if (ierr.ne.0) goto 998
936
      xmax=varmax(1)
937
      ymax=varmax(2)
938
      zmax=varmax(3)
939
      dz=(zmax-zmin)/real(nz-1)
940
      call clscdf(cstid,ierr)
941
      if (ierr.ne.0) goto 998
942
      call clscdf(cdfid,ierr)
943
      if (ierr.ne.0) goto 998
944
 
945
c     Check whether the grid is equallay spaced in the vertical
946
      do i=1,nz-1
947
         dh=aklev(i+1)-aklev(i)
948
         if (abs(dh-dz).gt.eps) then
949
            print*,'Aklev: Vertical grid must be equally spaced... Stop'
950
            print*,(aklev(i),i=1,nz)
951
            stop
952
         endif
953
         dh=aklay(i+1)-aklay(i)
954
         if (abs(dh-dz).gt.eps) then
955
            print*,'Aklay: Vertical grid must be equally spaced... Stop'
956
            print*,(aklay(i),i=1,nz)
957
            stop
958
         endif
959
      enddo
960
 
961
c     Set missing data value
962
      mdv=misdat
963
 
964
      return
965
 
966
c     Exception handling
967
 998  print*,'Read_Dim: Problem with input netcdf file... Stop'
968
      stop
969
 
970
      end