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 z2s
2
 
3
c     Calculate secondary fields on z levels
4
c     Michael Sprenger / Summer 2006
5
 
6
      implicit none
7
 
8
c     ---------------------------------------------------------------
9
c     Declaration of variables
10
c     ---------------------------------------------------------------
11
 
12
c     Variables for output Z file  : height level
13
      character*80 cfn
14
      real         varmin(4),varmax(4),stag(4)
15
      integer      vardim(4)
16
      real         mdv
17
      integer      ndim
18
      integer      nx,ny,nz
19
      real         xmin,xmax,ymin,ymax,dx,dy
20
      integer      ntimes
21
      real         aklev(1000),bklev(1000)
22
      real         aklay(1000),bklay(1000)
23
      real         time
24
      real         pollon,pollat
25
      integer      idate(5)
26
      integer      nfields
27
      character*80 field_nam(100)
28
      real,allocatable, dimension (:,:,:,:) :: field_dat
29
      real,allocatable, dimension (:,:,:)   :: z3
30
      real,allocatable, dimension (:,:)     :: x2,y2,f2
31
      real,allocatable, dimension (:,:,:)   :: out
32
      real,allocatable, dimension (:,:,:)   :: inp
33
      integer      nvars
34
      character*80 vnam(100),varname
35
      integer      isok
36
      integer      cdfid,cstid
37
 
38
c     Parameter file
39
      character*80 fieldname
40
      integer      nfilt
41
      character*80 ofn,gri
42
 
43
c     Auxiliary variables
44
      integer      ierr,stat
45
      integer      i,j,k,n
46
      real,allocatable, dimension (:,:)   :: out2
47
      character*80 vnam2(100)
48
      integer      nvars2
49
 
50
c     ---------------------------------------------------------------
51
c     Preparations
52
c     ---------------------------------------------------------------
53
 
54
      print*,'*********************************************************'
55
      print*,'* z2s                                                   *'
56
      print*,'*********************************************************'
57
 
58
c     Read parameter file
59
      open(10,file='fort.10')
60
       read(10,*) fieldname
61
       read(10,*) ofn
62
       read(10,*) gri
63
       read(10,*) nfilt
64
      close(10)
65
 
66
c     Get grid description for Z file : height level
67
      call cdfopn(ofn,cdfid,ierr)
68
      if (ierr.ne.0) goto 998
69
      call getcfn(cdfid,cfn,ierr)
70
      if (ierr.ne.0) goto 998
71
      call cdfopn(cfn,cstid,ierr)
72
      if (ierr.ne.0) goto 998
73
      call getdef(cdfid,'T',ndim,mdv,vardim,
74
     >            varmin,varmax,stag,ierr)
75
      if (ierr.ne.0) goto 998
76
      nx  =vardim(1)
77
      ny  =vardim(2)
78
      nz  =vardim(3)
79
      xmin=varmin(1)
80
      ymin=varmin(2)
81
      call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)
82
      if (ierr.ne.0) goto 998
83
      call getgrid(cstid,dx,dy,ierr)
84
      if (ierr.ne.0) goto 998
85
      xmax=xmin+real(nx-1)*dx
86
      ymax=ymin+real(ny-1)*dy
87
      call gettimes(cdfid,time,ntimes,ierr)
88
      if (ierr.ne.0) goto 998
89
      call getstart(cstid,idate,ierr)
90
      if (ierr.ne.0) goto 998
91
      call getpole(cstid,pollon,pollat,ierr)
92
      if (ierr.ne.0) goto 998
93
      call getvars(cdfid,nvars,vnam,ierr)
94
      if (ierr.ne.0) goto 998
95
      call clscdf(cstid,ierr)
96
      if (ierr.ne.0) goto 998
97
      call clscdf(cdfid,ierr)
98
      if (ierr.ne.0) goto 998
99
 
100
c     Get a list of all variables on the GRID file
101
      if (trim(gri).ne.trim(ofn)) then
102
         call cdfopn(gri,cdfid,ierr)
103
         if (ierr.ne.0) goto 998
104
         call getvars(cdfid,nvars2,vnam2,ierr)
105
         if (ierr.ne.0) goto 998
106
         do i=1,nvars2
107
            vnam(nvars+i)=vnam2(i)
108
         enddo
109
         nvars=nvars+nvars2
110
      endif
111
 
112
c     Check whether calculation of <fieldname> is supported
113
      if ( (fieldname.ne.'TH' ).and.
114
     >     (fieldname.ne.'NSQ')  .and.
115
     >     (fieldname.ne.'PV' )   .and.
116
     >     (fieldname.ne.'RHO')) then
117
         print*,'Variable not supported ',trim(fieldname)
118
         stop
119
      endif
120
 
121
c     Set dependencies
122
      if (fieldname.eq.'TH') then
123
         nfields=2
124
         field_nam(1)='T'
125
         field_nam(2)='P'
126
      else if (fieldname.eq.'RHO') then
127
         nfields=2
128
         field_nam(1)='T'
129
         field_nam(2)='P'
130
      else if (fieldname.eq.'NSQ') then
131
         nfields=2
132
         field_nam(1)='T'
133
         field_nam(2)='Q'
134
      else if (fieldname.eq.'PV') then
135
         nfields=4
136
         field_nam(1)='T'
137
         field_nam(2)='P'
138
         field_nam(3)='U'
139
         field_nam(4)='V'
140
      endif
141
 
142
c     Allocate memory
143
      allocate(field_dat(nfields,nx,ny,nz),stat=stat)
144
      if (stat.ne.0) print*,'error allocating field_dat'
145
      allocate(out(nx,ny,nz),stat=stat)
146
      if (stat.ne.0) print*,'error allocating out'
147
      allocate(inp(nx,ny,nz),stat=stat)
148
      if (stat.ne.0) print*,'error allocating inp'
149
      allocate(z3(nx,ny,nz),stat=stat)
150
      if (stat.ne.0) print*,'error allocating z3'
151
      allocate(x2(nx,ny),stat=stat)
152
      if (stat.ne.0) print*,'error allocating x2'
153
      allocate(y2(nx,ny),stat=stat)
154
      if (stat.ne.0) print*,'error allocating y2'
155
      allocate(f2(nx,ny),stat=stat)
156
      if (stat.ne.0) print*,'error allocating f2'
157
 
158
c     Allocate auxiliary fields
159
      allocate(out2(nx,ny),stat=stat)
160
      if (stat.ne.0) print*,'error allocating out2'
161
 
162
c     Read X grid
163
      varname='X'
164
      isok=0
165
      call check_varok (isok,varname,vnam,nvars)
166
      if (isok.eq.0) then
167
         print*,'Variable ',trim(varname),' missing... Stop'
168
         stop
169
      endif
170
      call cdfopn(gri,cdfid,ierr)
171
      if (ierr.ne.0) goto 998
172
      call getdat(cdfid,varname,time,0,x2,ierr)
173
      if (ierr.ne.0) goto 998
174
      print*,'R ',trim(varname),' ',trim(gri)
175
      call clscdf(cdfid,ierr)
176
      if (ierr.ne.0) goto 998
177
 
178
c     Read Y grid
179
      varname='Y'
180
      isok=0
181
      call check_varok (isok,varname,vnam,nvars)
182
      if (isok.eq.0) then
183
         print*,'Variable ',trim(varname),' missing... Stop'
184
         stop
185
      endif
186
      call cdfopn(gri,cdfid,ierr)
187
      if (ierr.ne.0) goto 998
188
      call getdat(cdfid,varname,time,0,y2,ierr)
189
      if (ierr.ne.0) goto 998
190
      print*,'R ',trim(varname),' ',trim(gri)
191
      call clscdf(cdfid,ierr)
192
      if (ierr.ne.0) goto 998
193
 
194
c     Read Coriolis parametzer
195
      varname='CORIOL'
196
      isok=0
197
      call check_varok (isok,varname,vnam,nvars)
198
      if (isok.eq.0) then
199
         print*,'Variable ',trim(varname),' missing... Stop'
200
         stop
201
      endif
202
      call cdfopn(gri,cdfid,ierr)
203
      if (ierr.ne.0) goto 998
204
      call getdat(cdfid,varname,time,0,f2,ierr)
205
      if (ierr.ne.0) goto 998
206
      print*,'R ',trim(varname),' ',trim(gri)
207
      call clscdf(cdfid,ierr)
208
      if (ierr.ne.0) goto 998      
209
 
210
c     Init height levels
211
      do i=1,nx
212
         do j=1,ny
213
            do k=1,nz
214
               z3(i,j,k)=aklay(k)
215
            enddo
216
         enddo
217
      enddo
218
 
219
c     Load needed variables
220
      do n=1,nfields
221
 
222
c        Check whether variable is available on file
223
         varname=field_nam(n)
224
         isok=0
225
         call check_varok (isok,varname,vnam,nvars)
226
 
227
c        Variable is available on file
228
         if (isok.eq.1) then
229
 
230
            call cdfopn(ofn,cdfid,ierr)
231
            if (ierr.ne.0) goto 998
232
            call getdat(cdfid,varname,time,0,inp,ierr)
233
            if (ierr.ne.0) goto 998
234
            print*,'R ',trim(varname),' ',trim(ofn)
235
            call clscdf(cdfid,ierr)
236
            if (ierr.ne.0) goto 998
237
 
238
            do i=1,nx
239
               do j=1,ny
240
                  do k=1,nz
241
                     field_dat(n,i,j,k)=inp(i,j,k)
242
                  enddo
243
               enddo
244
            enddo
245
 
246
         else            
247
            print*,'Variable missing : ',trim(varname)
248
            stop
249
         endif
250
 
251
      enddo
252
 
253
c     Change unit of pressure from hPa to Pa
254
      n=0
255
      do i=1,nfields
256
         if (field_nam(i).eq.'P') n=i
257
      enddo
258
      if (n.ne.0) then
259
         do i=1,nx
260
            do j=1,ny
261
               do k=1,nz
262
                  field_dat(n,i,j,k)=100.*field_dat(n,i,j,k)
263
               enddo
264
            enddo
265
         enddo
266
      endif
267
 
268
c     ----------------------------------------------------------------
269
c     Calculations
270
c     ----------------------------------------------------------------
271
 
272
c     Call to the defining routines
273
      if (fieldname.eq.'RHO') then
274
 
275
         print*,'C ',trim(fieldname)
276
         call calc_rho (out,                                 ! RHO
277
     >                  field_dat(1,:,:,:),                  ! T
278
     >                  field_dat(2,:,:,:),                  ! P
279
     >                  nx,ny,nz,mdv)
280
 
281
      else if (fieldname.eq.'TH') then
282
 
283
         print*,'C ',trim(fieldname)
284
         call calc_theta (out,                               ! TH
285
     >                    field_dat(1,:,:,:),                ! T
286
     >                    field_dat(2,:,:,:),                ! P
287
     >                    nx,ny,nz,mdv)
288
 
289
       else if (fieldname.eq.'NSQ') then
290
 
291
          print*,'C ',trim(fieldname)
292
          call calc_nsq (out,                                 ! NSQ
293
     >                   field_dat(1,:,:,:),                  ! T
294
     >                   field_dat(2,:,:,:),                  ! Q
295
     >                   z3,                                  ! Z      
296
     >                   nx,ny,nz,mdv)
297
 
298
       else if (fieldname.eq.'PV') then
299
 
300
          print*,'C ',trim(fieldname)
301
          call calc_pv (out,                                  ! PV
302
     >                  field_dat(1,:,:,:),                   ! T
303
     >                  field_dat(2,:,:,:),                   ! P
304
     >                  field_dat(3,:,:,:),                   ! U
305
     >                  field_dat(4,:,:,:),                   ! V
306
     >                  z3,x2,y2,f2,                          ! Z,X,Y,CORIOL
307
     >                  nx,ny,nz,mdv)
308
 
309
       endif
310
 
311
c      Horizontal filtering of the resulting fields
312
       print*,'F ',trim(fieldname)
313
       do k=1,nz
314
 
315
          do i=1,nx
316
             do j=1,ny
317
                out2(i,j)=out(i,j,k)
318
             enddo
319
          enddo
320
          do n=1,nfilt     
321
             call filt2d (out2,out2,nx,ny,1.,mdv,0,0,0,0)
322
          enddo
323
 
324
          do i=1,nx
325
             do j=1,ny
326
                out(i,j,k)=out2(i,j)
327
             enddo
328
          enddo
329
 
330
       enddo
331
 
332
c     ----------------------------------------------------------------
333
c     Save result onto netcdf file
334
c     ----------------------------------------------------------------
335
 
336
      call cdfwopn(ofn,cdfid,ierr)
337
      if (ierr.ne.0) goto 998
338
      isok=0
339
      varname=fieldname
340
      call check_varok(isok,varname,vnam,nvars)
341
      if (isok.eq.0) then
342
         call putdef(cdfid,varname,ndim,mdv,vardim,
343
     >               varmin,varmax,stag,ierr)
344
         if (ierr.ne.0) goto 998
345
      endif
346
      call putdat(cdfid,varname,time,0,out,ierr)
347
      if (ierr.ne.0) goto 998
348
      print*,'W ',trim(varname),' ',trim(ofn)
349
      call clscdf(cdfid,ierr)
350
      if (ierr.ne.0) goto 998
351
 
352
 
353
c     ----------------------------------------------------------------
354
c     Exception handling
355
c     ----------------------------------------------------------------
356
 
357
      stop
358
 
359
 998  print*,'Problem with netcdf file'
360
      stop
361
 
362
      end
363
 
364
 
365
 
366
c     ****************************************************************
367
c     * SUBROUTINE SECTION: AUXILIARY ROUTINES                       *
368
c     ****************************************************************     
369
 
370
c     ----------------------------------------------------------------
371
c     Check whether variable is found on netcdf file
372
c     ----------------------------------------------------------------
373
 
374
      subroutine check_varok (isok,varname,varlist,nvars)
375
 
376
c     Check whether the variable <varname> is in the list <varlist(nvars)>. 
377
c     If this is the case, <isok> is incremented by 1. Otherwise <isok> 
378
c     keeps its value.
379
 
380
      implicit none
381
 
382
c     Declaraion of subroutine parameters
383
      integer      isok
384
      integer      nvars
385
      character*80 varname
386
      character*80 varlist(nvars)
387
 
388
c     Auxiliary variables
389
      integer      i
390
 
391
c     Main
392
      do i=1,nvars
393
         if (trim(varname).eq.trim(varlist(i))) isok=isok+1
394
      enddo
395
 
396
      end
397
 
398
 
399
c     ****************************************************************
400
c     * SUBROUTINE SECTION: CALCULATE SECONDARY FIELDS               *
401
c     ****************************************************************
402
 
403
c     -----------------------------------------------------------------
404
c     Calculate density
405
c     -----------------------------------------------------------------
406
 
407
      subroutine calc_rho (rho3,t3,p3,
408
     >                     nx,ny,nz,mdv)
409
 
410
c     Calculate the density <rho3> (in kg/m^3) from temperature <t3> 
411
c     (in deg C) and pressure <p3> (in hPa).
412
 
413
      implicit none
414
 
415
c     Declaration of subroutine parameters
416
      integer   nx,ny,nz
417
      real      mdv
418
      real      rho3(nx,ny,nz)
419
      real      t3  (nx,ny,nz)
420
      real      p3  (nx,ny,nz)
421
 
422
c     Declaration of physical constants
423
      real      eps 
424
      parameter (eps=0.01)
425
      real      rd
426
      parameter (rd=287.)
427
      real      tzero
428
      parameter (tzero=273.15)
429
 
430
c     Auxiliary variables
431
      integer   i,j,k
432
 
433
c     Calculation
434
      do i=1,nx
435
         do j=1,ny
436
            do k=1,nz
437
 
438
               if ((abs(t3(i,j,k)-mdv).gt.eps).and.
439
     >             (abs(p3(i,j,k)-mdv).gt.eps)) then
440
 
441
                  rho3(i,j,k)=1./rd*p3(i,j,k)/(t3(i,j,k)+tzero)
442
 
443
               else
444
                  rho3(i,j,k)=mdv
445
               endif
446
 
447
            enddo
448
         enddo
449
      enddo
450
 
451
      end
452
 
453
 
454
c     -----------------------------------------------------------------
455
c     Calculate potential temperature
456
c     -----------------------------------------------------------------
457
 
458
      subroutine calc_theta (th3,t3,p3,
459
     >                       nx,ny,nz,mdv)
460
 
461
c     Calculate potential temperature <th3> (in K) from temperature <t3> 
462
c     (in deg C) and pressure <p3> (in hPa). 
463
 
464
      implicit none
465
 
466
c     Declaration of subroutine parameters
467
      integer   nx,ny,nz
468
      real      mdv
469
      real      th3(nx,ny,nz)
470
      real      t3 (nx,ny,nz)
471
      real      p3 (nx,ny,nz)
472
 
473
c     Declaration of physical constants
474
      real      rdcp,tzero,p0
475
      parameter (rdcp=0.286)
476
      parameter (tzero=273.15)
477
      parameter (p0=100000.)
478
      real      eps
479
      parameter (eps=0.01)
480
 
481
c     Auxiliary variables
482
      integer  i,j,k
483
 
484
c     Calculation
485
      do i=1,nx
486
         do j=1,ny
487
            do k=1,nz
488
 
489
               if ((abs(t3(i,j,k)-mdv).gt.eps).and.
490
     >             (abs(p3(i,j,k)-mdv).gt.eps)) then
491
 
492
                  th3(i,j,k)=(t3(i,j,k)+tzero)*( (p0/p3(i,j,k))**rdcp )
493
 
494
               else
495
                  th3(i,j,k)=mdv
496
               endif
497
 
498
            enddo
499
         enddo
500
      enddo
501
 
502
      end
503
 
504
c     -----------------------------------------------------------------
505
c     Calculate stratification
506
c     -----------------------------------------------------------------
507
 
508
      subroutine calc_nsq (nsq3,t3,q3,
509
     >                     z3,nx,ny,nz,mdv)
510
 
511
c     Calculate stratification <nsq3> on the model grid. The grid is
512
c     specified in the horizontal by <xmin,ymin,dx,dy,nx,ny>. The number
513
c     of vertical levels is <nz>. The input field are: temperature <t3>, 
514
c     specific humidity <q3>, height <z3> and horizontal grid <x2>, <y2>. 
515
 
516
      implicit none
517
 
518
c     Declaration of subroutine parameters
519
      integer   nx,ny,nz
520
      real      nsq3(nx,ny,nz)
521
      real      t3  (nx,ny,nz)
522
      real      q3  (nx,ny,nz)
523
      real      z3  (nx,ny,nz)
524
      real      x2  (nx,ny)
525
      real      y2  (nx,ny)
526
      real      mdv
527
 
528
c     Physical and numerical parameters
529
      real      scale
530
      parameter (scale=1.)
531
      real      g
532
      parameter (g=9.80665)
533
      real      eps
534
      parameter (eps=0.01)
535
      real      tzero
536
      parameter (tzero=273.15)
537
      real      kappa
538
      parameter (kappa=0.6078)
539
      real      cp
540
      parameter (cp=1005.)
541
      real      zerodiv
542
      parameter (zerodiv=0.0000000001)
543
      real      R
544
      parameter (R=287.)
545
 
546
c     Auxiliary variables
547
      real      tv   (nx,ny,nz)
548
      real      dtvdz(nx,ny,nz)
549
      integer   i,j,k
550
      real      scaledz
551
 
552
c     Calculate 3d virtual temperature
553
      do k=1,nz
554
         do i=1,nx
555
            do j=1,ny
556
               if ((abs(t3(i,j,k)-mdv).gt.eps).and.
557
     >             (abs(q3(i,j,k)-mdv).gt.eps)) 
558
     >         then
559
                  tv(i,j,k) = (t3(i,j,k)+tzero)*(1.+kappa*q3(i,j,k))
560
               else
561
                  tv(i,j,k) = mdv
562
               endif
563
            enddo
564
         enddo
565
      enddo
566
 
567
c     Vertical derivative of virtual temperature
568
      call  deriv (dtvdz,tv,'z',z3,x2,y2,nx,ny,nz,mdv)
569
 
570
c     Stratification
571
      do i=1,nx
572
         do j=1,ny
573
            do k=1,nz
574
               if ((abs(dtvdz(i,j,k)-mdv).gt.eps).and.
575
     >             (abs(tv   (i,j,k)-mdv).gt.eps)) 
576
     >         then
577
                  nsq3(i,j,k) = g/tv(i,j,k) * (dtvdz(i,j,k) + g/cp)  
578
               else
579
                  nsq3(i,j,k) = mdv
580
               endif
581
           enddo
582
        enddo
583
      enddo
584
 
585
      end
586
 
587
c     -----------------------------------------------------------------
588
c     Calculate potential vorticity
589
c     -----------------------------------------------------------------
590
 
591
      subroutine calc_pv (pv3,t3,p3,u3,v3,
592
     >                    z3,x2,y2,f2,nx,ny,nz,mdv)
593
 
594
c     Calculate potential vorticity <pv3> on the model grid. The grid is
595
c     specified in the horizontal by <xmin,ymin,dx,dy,nx,ny>. The number
596
c     of vertical levels is <nz>. The input field are: potential temperature
597
c     <th3>, horizontal wind <u3> and <v3>, density <rho3>, height <z3>. 
598
c     The terms including the vertical velocity are neglected in the calculation 
599
c     of the PV. 
600
 
601
      implicit none
602
 
603
c     Declaration of subroutine parameters
604
      integer   nx,ny,nz
605
      real      pv3 (nx,ny,nz)
606
      real      t3  (nx,ny,nz)
607
      real      p3  (nx,ny,nz)
608
      real      u3  (nx,ny,nz)
609
      real      v3  (nx,ny,nz)
610
      real      z3  (nx,ny,nz)
611
      real      x2  (nx,ny)
612
      real      y2  (nx,ny)
613
      real      f2  (nx,ny)
614
 
615
      real      mdv
616
 
617
c     Physical and numerical parameters
618
      real      scale
619
      parameter (scale=1.E6)
620
      real      omega
621
      parameter (omega=7.292E-5)
622
      real      pi180
623
      parameter (pi180=3.141592654/180.)
624
      real      eps
625
      parameter (eps=0.01)
626
 
627
c     Auxiliary variables
628
      real      dtdz(nx,ny,nz)
629
      real      dtdx(nx,ny,nz)
630
      real      dtdy(nx,ny,nz)
631
      real      dudz(nx,ny,nz)
632
      real      dvdz(nx,ny,nz)
633
      real      dvdx(nx,ny,nz)
634
      real      dudy(nx,ny,nz)
635
      real      rho3(nx,ny,nz)
636
      real      th3 (nx,ny,nz)
637
 
638
      integer   i,j,k
639
 
640
c     Calculate density and potential temperature
641
      call calc_rho   (rho3,t3,p3,nx,ny,nz,mdv)
642
      call calc_theta (th3 ,t3,p3,nx,ny,nz,mdv)
643
 
644
c     Calculate needed derivatives
645
      call deriv (dudz, u3,'z',z3,x2,y2,nx,ny,nz,mdv)
646
      call deriv (dvdz, v3,'z',z3,x2,y2,nx,ny,nz,mdv)
647
      call deriv (dtdz,th3,'z',z3,x2,y2,nx,ny,nz,mdv)
648
      call deriv (dtdx,th3,'x',z3,x2,y2,nx,ny,nz,mdv)
649
      call deriv (dtdy,th3,'y',z3,x2,y2,nx,ny,nz,mdv)     
650
      call deriv (dvdx, v3,'x',z3,x2,y2,nx,ny,nz,mdv)
651
      call deriv (dudy, u3,'y',z3,x2,y2,nx,ny,nz,mdv)     
652
 
653
c     Calculate potential vorticity
654
      do j=1,ny
655
         do i=1,nx
656
            do k=1,nz
657
 
658
c              Evaluate PV formula with missing data check
659
               if ((abs(dtdz(i,j,k)-mdv).gt.eps).and.
660
     >             (abs(dudz(i,j,k)-mdv).gt.eps).and.
661
     >             (abs(dtdy(i,j,k)-mdv).gt.eps).and.
662
     >             (abs(dvdz(i,j,k)-mdv).gt.eps).and.
663
     >             (abs(rho3(i,j,k)-mdv).gt.eps).and.              
664
     >             (abs(dtdx(i,j,k)-mdv).gt.eps).and.
665
     >             (abs(dvdx(i,j,k)-mdv).gt.eps).and.
666
     >             (abs(dudy(i,j,k)-mdv).gt.eps)) then
667
 
668
                  pv3(i,j,k)=scale*1./rho3(i,j,k)*
669
     >                   ((dvdx(i,j,k)-dudy(i,j,k)+f2(i,j))*dtdz(i,j,k)
670
     >                                         +dudz(i,j,k)*dtdy(i,j,k)
671
     >                                         -dvdz(i,j,k)*dtdx(i,j,k))
672
               else
673
                  pv3(i,j,k)=mdv
674
               endif
675
 
676
            enddo
677
         enddo
678
      enddo
679
 
680
      end
681
 
682
 
683
c     ****************************************************************
684
c     * SUBROUTINE SECTION: GRID HANDLING                            *
685
c     ****************************************************************
686
 
687
c     -----------------------------------------------------------------
688
c     Horizontal and vertical derivatives for 3d fields
689
c     -----------------------------------------------------------------
690
 
691
      subroutine deriv (df,f,direction,
692
     >                  z3,x2,y2,nx,ny,nz,mdv)
693
 
694
c     Calculate horizontal and vertical derivatives of the 3d field <f>.
695
c     The direction of the derivative is specified in <direction> 
696
c         'x','y'          : Horizontal derivative in x and y direction
697
c         'p','z','t','m'  : Vertical derivative (pressure, height, theta, model)
698
c     The 3d field <z3> specifies the isosurfaces along which the horizontal 
699
c     derivatives are calculated or the levels for the vertical derivatives.
700
 
701
      implicit none
702
 
703
c     Input and output parameters
704
      integer    nx,ny,nz
705
      real       df (nx,ny,nz)
706
      real       f  (nx,ny,nz)
707
      real       z3 (nx,ny,nz)
708
      real       x2 (nx,ny)
709
      real       y2 (nx,ny)
710
      character  direction
711
      real       mdv
712
 
713
c     Numerical and physical parameters
714
      real       pi180
715
      parameter  (pi180=3.141592654/180.)
716
      real       zerodiv
717
      parameter  (zerodiv=0.00000001)
718
      real       eps
719
      parameter  (eps=0.01) 
720
 
721
c     Auxiliary variables
722
      integer    i,j,k
723
      real       vmin,vmax
724
      real       scale,lat
725
      real       vu,vl,vuvl,vlvu
726
      integer    o,u,w,e,n,s
727
 
728
c     Vertical derivative
729
      if ((direction.eq.'z').or.
730
     >    (direction.eq.'th').or.
731
     >    (direction.eq.'p').or.
732
     >    (direction.eq.'m').and.
733
     >    (nz.gt.1)) then
734
 
735
         do i=1,nx
736
            do j=1,ny
737
               do k=1,nz
738
 
739
                  o=k+1
740
                  if (o.gt.nz) o=nz
741
                  u=k-1
742
                  if (u.lt.1) u=1                  
743
 
744
                  if ((abs(f(i,j,o)-mdv).gt.eps).and.
745
     >                (abs(f(i,j,u)-mdv).gt.eps).and.
746
     >                (abs(f(i,j,k)-mdv).gt.eps)) then
747
 
748
                     vu = z3(i,j,k)-z3(i,j,o)
749
                     vl = z3(i,j,u)-z3(i,j,k)
750
                     vuvl = vu/(vl+zerodiv)
751
                     vlvu = 1./(vuvl+zerodiv)
752
 
753
                     df(i,j,k) = 1./(vu+vl)
754
     >                           * (vuvl*(f(i,j,u)-f(i,j,k)) 
755
     >                           +  vlvu*(f(i,j,k)-f(i,j,o))) 
756
 
757
                  else
758
                     df(i,j,k) = mdv
759
                  endif
760
 
761
               enddo
762
            enddo
763
         enddo
764
 
765
c     Horizontal derivative in the y direction: 3d
766
      elseif (direction.eq.'y') then
767
 
768
         do i=1,nx
769
            do j=1,ny
770
               do k=1,nz
771
 
772
                  s=j-1
773
                  if (s.lt.1) s=1
774
                  n=j+1
775
                  if (n.gt.ny) n=ny
776
 
777
                  if ((abs(f(i,n,k)-mdv).gt.eps).and.
778
     >                (abs(f(i,j,k)-mdv).gt.eps).and.
779
     >                (abs(f(i,s,k)-mdv).gt.eps)) then  
780
 
781
                     vu = 1000.*(y2(i,j)-y2(i,n))
782
                     vl = 1000.*(y2(i,s)-y2(i,j))
783
                     vuvl = vu/(vl+zerodiv)
784
                     vlvu = 1./(vuvl+zerodiv)
785
 
786
                     df(i,j,k) =  1./(vu+vl)
787
     >                            * (vuvl*(f(i,s,k)-f(i,j,k))
788
     >                            +  vlvu*(f(i,j,k)-f(i,n,k)))
789
 
790
                  else
791
                     df(i,j,k) = mdv
792
                  endif
793
 
794
               enddo
795
            enddo
796
         enddo
797
 
798
c     Horizontal derivative in the x direction: 3d
799
      elseif (direction.eq.'x') then
800
 
801
         do i=1,nx
802
            do j=1,ny
803
               do k=1,nz
804
 
805
                  w=i-1
806
                  if (w.lt.1) w=1
807
                  e=i+1
808
                  if (e.gt.nx) e=nx
809
 
810
                  if ((abs(f(w,j,k)-mdv).gt.eps).and.
811
     >                (abs(f(i,j,k)-mdv).gt.eps).and.
812
     >                (abs(f(e,j,k)-mdv).gt.eps)) then  
813
 
814
                     vu = 1000.*(x2(i,j)-x2(e,j))
815
                     vl = 1000.*(x2(w,j)-x2(i,j))
816
                     vuvl = vu/(vl+zerodiv)
817
                     vlvu = 1./(vuvl+zerodiv)
818
 
819
                     df(i,j,k) =  1./(vu+vl)
820
     >                            * (vuvl*(f(w,j,k)-f(i,j,k))
821
     >                            +  vlvu*(f(i,j,k)-f(e,j,k)))
822
 
823
                  else
824
                     df(i,j,k) = mdv
825
                  endif
826
 
827
               enddo
828
            enddo
829
         enddo
830
 
831
c     Undefined direction for derivative
832
      else
833
 
834
         print*,'Invalid direction of derivative... Stop'
835
         stop
836
 
837
      endif
838
 
839
      end
840
 
841
c     -----------------------------------------------------------------
842
c     Horizontal filter
843
c     -----------------------------------------------------------------
844
 
845
      subroutine filt2d (a,af,nx,ny,fil,misdat,
846
     &                   iperx,ipery,ispol,inpol)
847
 
848
c     Apply a conservative diffusion operator onto the 2d field a,
849
c     with full missing data checking.
850
c
851
c     a     real   inp  array to be filtered, dimensioned (nx,ny)
852
c     af    real   out  filtered array, dimensioned (nx,ny), can be
853
c                       equivalenced with array a in the calling routine
854
c     f1    real        workarray, dimensioned (nx+1,ny)
855
c     f2    real        workarray, dimensioned (nx,ny+1)
856
c     fil   real   inp  filter-coeff., 0<afil<=1. Maximum filtering with afil=1
857
c                       corresponds to one application of the linear filter.
858
c     misdat real  inp  missing-data value, a(i,j)=misdat indicates that
859
c                       the corresponding value is not available. The
860
c                       misdat-checking can be switched off with with misdat=0.
861
c     iperx int    inp  periodic boundaries in the x-direction (1=yes,0=no)
862
c     ipery int    inp  periodic boundaries in the y-direction (1=yes,0=no)
863
c     inpol int    inp  northpole at j=ny  (1=yes,0=no)
864
c     ispol int    inp  southpole at j=1   (1=yes,0=no)
865
c
866
c     Christoph Schaer, 1993
867
 
868
c     argument declaration
869
      integer     nx,ny
870
      real        a(nx,ny),af(nx,ny),fil,misdat
871
      integer     iperx,ipery,inpol,ispol
872
 
873
c     local variable declaration
874
      integer     i,j,is
875
      real        fh
876
      real        f1(nx+1,ny),f2(nx,ny+1)
877
 
878
c     compute constant fh
879
      fh=0.125*fil
880
 
881
c     compute fluxes in x-direction
882
      if (misdat.eq.0.) then
883
        do j=1,ny
884
        do i=2,nx
885
          f1(i,j)=a(i-1,j)-a(i,j)
886
        enddo
887
        enddo
888
      else
889
        do j=1,ny
890
        do i=2,nx
891
          if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then
892
            f1(i,j)=0.
893
          else
894
            f1(i,j)=a(i-1,j)-a(i,j)
895
          endif
896
        enddo
897
        enddo
898
      endif
899
      if (iperx.eq.1) then
900
c       do periodic boundaries in the x-direction
901
        do j=1,ny
902
          f1(1,j)=f1(nx,j)
903
          f1(nx+1,j)=f1(2,j)
904
        enddo
905
      else
906
c       set boundary-fluxes to zero
907
        do j=1,ny
908
          f1(1,j)=0.
909
          f1(nx+1,j)=0.
910
        enddo
911
      endif
912
 
913
c     compute fluxes in y-direction
914
      if (misdat.eq.0.) then
915
        do j=2,ny
916
        do i=1,nx
917
          f2(i,j)=a(i,j-1)-a(i,j)
918
        enddo
919
        enddo
920
      else
921
        do j=2,ny
922
        do i=1,nx
923
          if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then
924
            f2(i,j)=0.
925
          else
926
            f2(i,j)=a(i,j-1)-a(i,j)
927
          endif
928
        enddo
929
        enddo
930
      endif
931
c     set boundary-fluxes to zero
932
      do i=1,nx
933
        f2(i,1)=0.
934
        f2(i,ny+1)=0.
935
      enddo
936
      if (ipery.eq.1) then
937
c       do periodic boundaries in the x-direction
938
        do i=1,nx
939
          f2(i,1)=f2(i,ny)
940
          f2(i,ny+1)=f2(i,2)
941
        enddo
942
      endif
943
      if (iperx.eq.1) then
944
        if (ispol.eq.1) then
945
c         do south-pole
946
          is=(nx-1)/2
947
          do i=1,nx
948
            f2(i,1)=-f2(mod(i-1+is,nx)+1,2)
949
          enddo
950
        endif
951
        if (inpol.eq.1) then
952
c         do north-pole
953
          is=(nx-1)/2
954
          do i=1,nx
955
            f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny)
956
          enddo
957
        endif
958
      endif
959
 
960
c     compute flux-convergence -> filter
961
      if (misdat.eq.0.) then
962
        do j=1,ny
963
        do i=1,nx
964
            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
965
        enddo
966
        enddo
967
      else
968
        do j=1,ny
969
        do i=1,nx
970
          if (a(i,j).eq.misdat) then
971
            af(i,j)=misdat
972
          else
973
            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
974
          endif
975
        enddo
976
        enddo
977
      endif
978
      end