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 main_calc_qgpv
2
 
3
c     ********************************************************************************
4
c     * CALCULATE QUASI-GEOSTROPHIC PV ACCORDING TO ITS DEFINITION                   *
5
c     * Michael Sprenger / Summer, Autumn 2006                                       *
6
c     ********************************************************************************
7
 
8
c     --------------------------------------------------------------------------------
9
c     Declaration of variables, parameters, externals and common blocks
10
c     --------------------------------------------------------------------------------
11
 
12
      implicit none
13
 
14
c     Input and output file
15
      character*80   diagfile
16
      character*80   referfile
17
 
18
c     Grid parameters
19
      integer        nx,ny,nz
20
      real           xmin,ymin,zmin,xmax,ymax,zmax
21
      real           dx,dy,dz
22
      real           deltax,deltay,deltaz
23
      real           mdv
24
      real,          allocatable,dimension (:,:) :: coriol
25
 
26
c     Reference state 
27
      real,   allocatable,dimension (:) :: nsqref
28
      real,   allocatable,dimension (:) :: thetaref
29
      real,   allocatable,dimension (:) :: rhoref
30
      real,   allocatable,dimension (:) :: pressref
31
      real,   allocatable,dimension (:) :: zref
32
 
33
c     3d fields for calculation of qgPV 
34
      real,allocatable,dimension (:,:,:) :: qgpv,th,uu,vv
35
 
36
c     Auxiliary variables
37
      integer      i,j,k
38
      integer      stat
39
      character*80 varname
40
 
41
c     --------------------------------------------------------------------------------
42
c     Input
43
c     --------------------------------------------------------------------------------
44
 
45
      print*,'********************************************************'
46
      print*,'* CALC_QGPV                                            *'
47
      print*,'********************************************************'
48
 
49
c     Read parameter file
50
      open(10,file='fort.10')
51
       read(10,*) diagfile
52
       read(10,*) referfile
53
      close(10) 
54
      print*
55
      print*,trim(diagfile)
56
      print*,trim(referfile)
57
      print*
58
 
59
c     Get lat/lon gid parameters from input file
60
      call read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv,
61
     >               diagfile)
62
      print*,'Read_Dim: nx,ny,nz         = ',nx,ny,nz
63
      print*,'          dx,dy,dz         = ',dx,dy,dz
64
      print*,'          xmin,ymin,zmin   = ',xmin,ymin,zmin
65
      print*,'          mdv              = ',mdv
66
      print*
67
 
68
c     Count from 0, not from 1: consistent with <inv_cart.f>.
69
      nx=nx-1
70
      ny=ny-1
71
      nz=nz-1
72
 
73
c     Allocate memory for Coriolis parameters
74
      allocate(coriol(0:nx,0:ny),STAT=stat)
75
      if (stat.ne.0) print*,'error allocating coriol'
76
 
77
c     Allocate memory for reference profile
78
      allocate(rhoref  (0:2*nz),STAT=stat)
79
      if (stat.ne.0) print*,'error allocating rhoref'
80
      allocate(pressref(0:2*nz),STAT=stat)
81
      if (stat.ne.0) print*,'error allocating pressref'
82
      allocate(thetaref(0:2*nz),STAT=stat)
83
      if (stat.ne.0) print*,'error allocating thetaref'
84
      allocate(nsqref  (0:2*nz),STAT=stat)
85
      if (stat.ne.0) print*,'error allocating nsqref'
86
      allocate(zref    (0:2*nz),STAT=stat)
87
      if (stat.ne.0) print*,'error allocating zref'
88
 
89
c     Allocate memory for 3d fields
90
      allocate(qgpv(0:nx,0:ny,0:nz),STAT=stat)
91
      if (stat.ne.0) print*,'error allocating qgpv'
92
      allocate(uu(0:nx,0:ny,0:nz),STAT=stat)
93
      if (stat.ne.0) print*,'error allocating uu'
94
      allocate(vv(0:nx,0:ny,0:nz),STAT=stat)
95
      if (stat.ne.0) print*,'error allocating vv'
96
      allocate(th(0:nx,0:ny,0:nz),STAT=stat)
97
      if (stat.ne.0) print*,'error allocating th'
98
 
99
 
100
c     --------------------------------------------------------------------------------
101
c     Calculate the qgPV from definition and put it onto file
102
c     --------------------------------------------------------------------------------
103
 
104
c     Read data from file
105
      varname='U'
106
      call read_inp (uu,varname,diagfile,
107
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
108
      varname='V'
109
      call read_inp (vv,varname,diagfile,
110
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
111
      varname='TH'
112
      call read_inp (th,varname,diagfile,
113
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
114
 
115
c     Read reference profile and grid parameters
116
      call read_ref (nsqref,rhoref,thetaref,pressref,zref,
117
     >               nx,ny,nz,deltax,deltay,deltaz,coriol,
118
     >               referfile)
119
      deltax=1000.*deltax
120
      deltay=1000.*deltay
121
      print*,'Deltax,deltay,deltaz =',deltax,deltay,deltaz
122
 
123
c     Calculate qgPV
124
      print*,'C qgPV'
125
      call calc_qgpv (qgpv,uu,vv,th,
126
     >                rhoref,pressref,nsqref,thetaref,coriol,
127
     >                nx,ny,nz,deltax,deltay,deltaz,mdv)
128
 
129
c     Write result to netcdf file
130
      varname='QGPV_DIA'
131
      call write_inp (qgpv,varname,diagfile,nx,ny,nz)
132
 
133
      end
134
 
135
c     ********************************************************************************
136
c     * NETCDF OUTPUT                                                                *
137
c     ********************************************************************************
138
 
139
c     --------------------------------------------------------------------------------
140
c     Write input field to netcdf
141
c     --------------------------------------------------------------------------------
142
 
143
      SUBROUTINE write_inp (field,fieldname,pvsrcfile,nx,ny,nz)
144
 
145
c     Read <fieldname> from netcdf file <pvsrcfile> into <field>. The grid is specified 
146
c     by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed whether the input 
147
c     files are consitent with this grid. 
148
 
149
      implicit   none
150
 
151
c     Declaration of subroutine parameters
152
      integer              nx,ny,nz
153
      real                 field (0:nx,0:ny,0:nz)
154
      character*80         fieldname
155
      character*80         pvsrcfile
156
 
157
c     Auxiliary variables
158
      integer              cdfid,stat
159
      integer              vardim(4)
160
      real                 misdat
161
      real                 varmin(4),varmax(4),stag(4)
162
      integer              ndimin,outid,i,j,k
163
      real                 max_th
164
      real                 tmp(0:nx,0:ny,0:nz)
165
      integer              ntimes
166
      real                 time(200)       
167
      integer              nvars
168
      character*80         vnam(100),varname
169
      integer              isok
170
 
171
c     Get grid parameters from THETA
172
      call cdfopn(pvsrcfile,cdfid,stat)
173
      if (stat.ne.0) goto 998
174
      call getvars(cdfid,nvars,vnam,stat)
175
      if (stat.ne.0) goto 998
176
      isok=0
177
      varname='TH'
178
      call check_varok(isok,varname,vnam,nvars)
179
      if (isok.eq.0) goto 998
180
      call getdef(cdfid,varname,ndimin,misdat,vardim,
181
     >            varmin,varmax,stag,stat)    
182
      if (stat.ne.0) goto 998
183
      time(1)=0.
184
      call gettimes(cdfid,time,ntimes,stat)  
185
      if (stat.ne.0) goto 998
186
      call clscdf(cdfid,stat)
187
      if (stat.ne.0) goto 998
188
 
189
c     Save variables (write definition, if necessary)
190
      call cdfwopn(pvsrcfile,cdfid,stat)
191
      if (stat.ne.0) goto 998
192
      isok=0
193
      varname=fieldname
194
      call check_varok(isok,varname,vnam,nvars)
195
      if (isok.eq.0) then
196
         call putdef(cdfid,varname,ndimin,misdat,vardim,
197
     >               varmin,varmax,stag,stat)
198
         if (stat.ne.0) goto 998
199
      endif
200
      call putdat(cdfid,varname,time(1),0,field,stat)
201
      print*,'W ',trim(varname),' ',trim(pvsrcfile)
202
      if (stat.ne.0) goto 998
203
 
204
c     Close input netcdf file
205
      call clscdf(cdfid,stat)
206
      if (stat.ne.0) goto 998
207
 
208
      return
209
 
210
c     Exception handling
211
 998  print*,'Write_Inp: Problem with input netcdf file... Stop'
212
      stop
213
 
214
      end
215
 
216
c     --------------------------------------------------------------------------------
217
c     Read input fields 
218
c     --------------------------------------------------------------------------------
219
 
220
      SUBROUTINE read_inp (field,fieldname,pvsrcfile,
221
     >                     nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
222
 
223
c     Read <fieldname> from netcdf file <pvsrcfile> into <field>. The grid is specified 
224
c     by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed whether the input 
225
c     files are consitent with this grid. The missing data value is set to <mdv>.
226
 
227
      implicit   none
228
 
229
c     Declaration of subroutine parameters
230
      integer              nx,ny,nz
231
      real                 field(0:nx,0:ny,0:nz)
232
      character*80         fieldname
233
      character*80         pvsrcfile
234
      real                 dx,dy,dz,xmin,ymin,zmin
235
      real                 mdv
236
 
237
c     Numerical and physical parameters
238
      real                 eps
239
      parameter            (eps=0.01)
240
 
241
c     Auxiliary variables
242
      integer              cdfid,stat,cdfid99
243
      integer              vardim(4)
244
      real                 misdat
245
      real                 varmin(4),varmax(4),stag(4)
246
      integer              ndimin,outid,i,j,k
247
      real                 max_th
248
      real                 tmp(nx,ny,nz)
249
      integer              ntimes
250
      real                 time(200)       
251
      integer              nvars
252
      character*80         vnam(100),varname
253
      integer              isok
254
 
255
c     Open the input netcdf file
256
      call cdfopn(pvsrcfile,cdfid,stat)
257
      if (stat.ne.0) goto 998
258
 
259
c     Check whether needed variables are on file
260
      call getvars(cdfid,nvars,vnam,stat)
261
      if (stat.ne.0) goto 998
262
      isok=0
263
      varname=trim(fieldname)
264
      call check_varok(isok,varname,vnam,nvars)
265
      if (isok.eq.0) goto 998
266
 
267
c     Get the grid parameters from theta     
268
      call getdef(cdfid,varname,ndimin,misdat,vardim,
269
     >            varmin,varmax,stag,stat)    
270
      if (stat.ne.0) goto 998
271
      time(1)=0.
272
      call gettimes(cdfid,time,ntimes,stat)  
273
      if (stat.ne.0) goto 998
274
 
275
c     Check whether grid parameters are consistent
276
      if ( (vardim(1).ne.(nx+1)).or.
277
     >     (vardim(2).ne.(ny+1)).or.
278
     >     (vardim(3).ne.(nz+1)).or.
279
     >     (abs(varmin(1)-xmin).gt.eps).or.
280
     >     (abs(varmin(2)-ymin).gt.eps).or.
281
     >     (abs(varmin(3)-zmin).gt.eps).or.
282
     >     (abs((varmax(1)-varmin(1))/real(vardim(1)-1)-dx).gt.eps).or.
283
     >     (abs((varmax(2)-varmin(2))/real(vardim(2)-1)-dy).gt.eps).or.
284
     >     (abs((varmax(3)-varmin(3))/real(vardim(3)-1)-dz).gt.eps) ) 
285
     >then
286
         print*,'Input grid inconsitency...'
287
         print*,'  Nx      = ',vardim(1),nx+1
288
         print*,'  Ny      = ',vardim(2),ny+1
289
         print*,'  Nz      = ',vardim(3),nz+1
290
         print*,'  Varminx = ',varmin(1),xmin
291
         print*,'  Varminy = ',varmin(2),ymin
292
         print*,'  Varminz = ',varmin(3),zmin
293
         print*,'  Dx      = ',(varmax(1)-varmin(1))/real(nx-1),dx
294
         print*,'  Dy      = ',(varmax(2)-varmin(2))/real(ny-1),dy
295
         print*,'  Dz      = ',(varmax(3)-varmin(3))/real(nz-1),dz
296
         goto 998
297
      endif
298
 
299
c     Load variables
300
      call getdef(cdfid,varname,ndimin,misdat,vardim,
301
     >            varmin,varmax,stag,stat)
302
      if (stat.ne.0) goto 998
303
      call getdat(cdfid,varname,time(1),0,field,stat)
304
      print*, 'R ',trim(varname),' ',trim(pvsrcfile)
305
      if (stat.ne.0) goto 998
306
 
307
c     Close input netcdf file
308
      call clscdf(cdfid,stat)
309
      if (stat.ne.0) goto 998
310
 
311
c     Set missing data value to <mdv>
312
      do i=1,nx
313
         do j=1,ny
314
            do k=1,nz
315
               if (abs(field(i,j,k)-misdat).lt.eps) then
316
                  field(i,j,k)=mdv
317
               endif
318
            enddo
319
         enddo
320
      enddo
321
 
322
      return
323
 
324
c     Exception handling
325
 998  print*,'Read_Inp: Problem with input netcdf file... Stop'
326
      stop
327
 
328
      end
329
 
330
c     --------------------------------------------------------------------------------
331
c     Read refernece profile from netcdf
332
c     --------------------------------------------------------------------------------
333
 
334
      SUBROUTINE read_ref (nsqref,rhoref,thetaref,pressref,zref,
335
     >                     nx,ny,nz,deltax,deltay,deltaz,coriol,
336
     >                     pvsrcfile) 
337
 
338
c     Read the reference profile from file
339
c
340
c         thetaref             : Reference potential temperature (K)
341
c         pressref             : Reference pressure (Pa)
342
c         rhoref               : Reference density (kg/m^3)
343
c         nsqref               : Stratification (s^-1)
344
c         zref                 : Reference height (m)
345
c         nx,nny,nz            : Grid dimension in x,y,z direction
346
c         deltax,deltay,deltaz : Grid spacings used for calculations (m)
347
c         coriol               : Coriolis parameter (s^-1)
348
c         pvsrcfile            : Input file
349
 
350
      implicit   none
351
 
352
c     Declaration of subroutine parameters
353
      integer              nx,ny,nz
354
      real                 nsqref  (0:2*nz)
355
      real                 thetaref(0:2*nz)
356
      real                 rhoref  (0:2*nz)
357
      real                 pressref(0:2*nz)
358
      real                 zref    (0:2*nz)
359
      real                 deltax,deltay,deltaz
360
      real                 coriol  (0:nx,0:ny)
361
      character*80         pvsrcfile
362
 
363
c     Numerical and physical parameters
364
      real                 eps
365
      parameter            (eps=0.01)
366
 
367
c     Auxiliary variables
368
      integer              cdfid,stat
369
      integer              vardim(4)
370
      real                 misdat
371
      integer              ndimin
372
      real                 varmin(4),varmax(4),stag(4)
373
      integer              i,j,k,nf1
374
      integer              ntimes
375
      real                 time(200)  
376
      character*80         vnam(100),varname
377
      integer              nvars
378
      integer              isok,ierr
379
      real                 x(0:nx,0:ny),y(0:nx,0:ny)
380
      real                 mean,count
381
 
382
c     Get grid description from topography
383
      call cdfopn(pvsrcfile,cdfid,stat)
384
      if (stat.ne.0) goto 997
385
      call getvars(cdfid,nvars,vnam,stat)
386
      if (stat.ne.0) goto 997
387
      isok=0
388
      varname='ORO'
389
      call check_varok(isok,varname,vnam,nvars)
390
      if (isok.eq.0) goto 997
391
      call getdef(cdfid,varname,ndimin,misdat,vardim,
392
     >            varmin,varmax,stag,stat)    
393
      if (stat.ne.0) goto 997
394
      time(1)=0.
395
      call gettimes(cdfid,time,ntimes,stat)  
396
      if (stat.ne.0) goto 997
397
      call clscdf(cdfid,stat)
398
      if (stat.ne.0) goto 997
399
 
400
c     Open output netcdf file
401
      call cdfopn(pvsrcfile,cdfid,stat)
402
      if (stat.ne.0) goto 997
403
 
404
c     Create the variable if necessary
405
      call getvars(cdfid,nvars,vnam,stat)
406
      if (stat.ne.0) goto 997
407
 
408
c     Read data from netcdf file
409
      isok=0
410
      varname='NSQREF'
411
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
412
      call check_varok(isok,varname,vnam,nvars)
413
      if (isok.eq.0) goto 997
414
      call getdat(cdfid,varname,time(1),0,nsqref,stat)
415
      if (stat.ne.0) goto 997
416
 
417
      isok=0
418
      varname='RHOREF'
419
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
420
      call check_varok(isok,varname,vnam,nvars)
421
      if (isok.eq.0) goto 997
422
      call getdat(cdfid,varname,time(1),0,rhoref,stat)
423
      if (stat.ne.0) goto 997
424
 
425
      isok=0
426
      varname='THETAREF'
427
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
428
      call check_varok(isok,varname,vnam,nvars)
429
      if (isok.eq.0) goto 997
430
      call getdat(cdfid,varname,time(1),0,thetaref,stat)
431
      if (stat.ne.0) goto 997
432
 
433
      isok=0
434
      varname='PREREF'
435
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
436
      call check_varok(isok,varname,vnam,nvars)
437
      if (isok.eq.0) goto 997
438
      call getdat(cdfid,varname,time(1),0,pressref,stat)
439
      if (stat.ne.0) goto 997
440
 
441
      isok=0
442
      varname='ZREF'
443
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
444
      call check_varok(isok,varname,vnam,nvars)
445
      if (isok.eq.0) goto 997
446
      call getdat(cdfid,varname,time(1),0,zref,stat)
447
      if (stat.ne.0) goto 997
448
 
449
      isok=0
450
      varname='CORIOL'
451
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
452
      call check_varok(isok,varname,vnam,nvars)
453
      if (isok.eq.0) goto 997
454
      call getdat(cdfid,varname,time(1),0,coriol,stat)
455
      if (stat.ne.0) goto 997
456
 
457
      isok=0
458
      varname='X'
459
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
460
      call check_varok(isok,varname,vnam,nvars)
461
      if (isok.eq.0) goto 997
462
      call getdat(cdfid,varname,time(1),0,x,stat)
463
      if (stat.ne.0) goto 997
464
 
465
      isok=0
466
      varname='Y'
467
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
468
      call check_varok(isok,varname,vnam,nvars)
469
      if (isok.eq.0) goto 997
470
      call getdat(cdfid,varname,time(1),0,y,stat)
471
      if (stat.ne.0) goto 997
472
 
473
c     Close  netcdf file
474
      call clscdf(cdfid,stat)
475
      if (stat.ne.0) goto 997
476
 
477
c     Determine the grid spacings <deltax, deltay, deltaz>
478
      mean=0.
479
      count=0.
480
      do i=1,nx
481
         do j=0,ny
482
            mean=mean+abs(x(i,j)-x(i-1,j))
483
            count=count+1.
484
         enddo
485
      enddo
486
      deltax=mean/count
487
 
488
      mean=0.
489
      count=0.
490
      do j=1,ny
491
         do i=0,nx
492
            mean=mean+abs(y(i,j)-y(i,j-1))
493
            count=count+1.
494
         enddo
495
      enddo
496
      deltay=mean/count
497
 
498
      mean=0.
499
      count=0.
500
      do k=1,nz-1
501
         mean=mean+abs(zref(k+1)-zref(k-1))
502
         count=count+1.
503
      enddo
504
      deltaz=mean/count
505
 
506
      return
507
 
508
c     Exception handling
509
 997  print*,'Read_Ref: Problem with input netcdf file... Stop'
510
      stop
511
 
512
      end
513
 
514
 
515
c     --------------------------------------------------------------------------------
516
c     Check whether variable is found on netcdf file
517
c     --------------------------------------------------------------------------------
518
 
519
      subroutine check_varok (isok,varname,varlist,nvars)
520
 
521
c     Check whether the variable <varname> is in the list <varlist(nvars)>. If this is 
522
C     the case, <isok> is incremented by 1. Otherwise <isok> keeps its value.
523
 
524
      implicit none
525
 
526
c     Declaraion of subroutine parameters
527
      integer      isok
528
      integer      nvars
529
      character*80 varname
530
      character*80 varlist(nvars)
531
 
532
c     Auxiliary variables
533
      integer      i
534
 
535
c     Main
536
      do i=1,nvars
537
         if (trim(varname).eq.trim(varlist(i))) isok=isok+1
538
      enddo
539
 
540
      end
541
 
542
c     --------------------------------------------------------------------------------
543
c     Get grid parameters
544
c     --------------------------------------------------------------------------------
545
 
546
      subroutine read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv,
547
     >                     pvsrcfile)
548
 
549
c     Get the grid parameters from the variable <THETA> on the input file <pvsrcfile>.
550
c     The grid parameters are
551
c        nx,ny,nz                : Number of grid points in x, y and z direction
552
c        xmin,ymin,zmin          : Minimum domain coordinates in x, y and z direction
553
c        xmax,ymax,zmax          : Maximal domain coordinates in x, y and z direction
554
c        dx,dy,dz                : Horizontal and vertical resolution
555
c     Additionally, it is checked whether the vertical grid is equally spaced. If ok,
556
c     the grid paramters are transformed from lon/lat to distance (in meters)
557
 
558
      implicit none
559
 
560
c     Declaration of subroutine parameters
561
      character*80   pvsrcfile
562
      integer        nx,ny,nz
563
      real           dx,dy,dz
564
      real           xmin,ymin,zmin,xmax,ymax,zmax
565
      real           mdv
566
 
567
c     Numerical epsilon and other physical/geoemtrical parameters
568
      real           eps
569
      parameter      (eps=0.01)
570
 
571
c     Auxiliary variables
572
      integer        cdfid,cstid
573
      integer        ierr
574
      character*80   vnam(100),varname
575
      integer        nvars
576
      integer        isok
577
      integer        vardim(4)
578
      real           misdat
579
      real           varmin(4),varmax(4),stag(4)
580
      real           aklev(1000),bklev(1000),aklay(1000),bklay(1000)
581
      real           dh
582
      character*80   csn
583
      integer        ndim
584
      integer        i
585
 
586
c     Get all grid parameters
587
      call cdfopn(pvsrcfile,cdfid,ierr)
588
      if (ierr.ne.0) goto 998
589
      call getvars(cdfid,nvars,vnam,ierr)
590
      if (ierr.ne.0) goto 998
591
      isok=0
592
      varname='TH'
593
      call check_varok(isok,varname,vnam,nvars)
594
      if (isok.eq.0) goto 998
595
      call getcfn(cdfid,csn,ierr)
596
      if (ierr.ne.0) goto 998
597
      call cdfopn(csn,cstid,ierr)
598
      if (ierr.ne.0) goto 998
599
      call getdef(cdfid,varname,ndim,misdat,vardim,varmin,varmax,
600
     >            stag,ierr)
601
      if (ierr.ne.0) goto 998
602
      nx=vardim(1)
603
      ny=vardim(2)
604
      nz=vardim(3)
605
      xmin=varmin(1)
606
      ymin=varmin(2)
607
      zmin=varmin(3)
608
      call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)
609
      if (ierr.ne.0) goto 998
610
      call getgrid(cstid,dx,dy,ierr)
611
      if (ierr.ne.0) goto 998
612
      xmax=varmax(1)
613
      ymax=varmax(2)
614
      zmax=varmax(3)
615
      dz=(zmax-zmin)/real(nz-1)
616
      call clscdf(cstid,ierr)
617
      if (ierr.ne.0) goto 998
618
      call clscdf(cdfid,ierr)
619
      if (ierr.ne.0) goto 998
620
 
621
c     Check whether the grid is equallay spaced in the vertical
622
      do i=1,nz-1
623
         dh=aklev(i+1)-aklev(i)
624
         if (abs(dh-dz).gt.eps) then
625
            print*,'Aklev: Vertical grid must be equally spaced... Stop'
626
            print*,(aklev(i),i=1,nz)
627
            stop
628
         endif
629
         dh=aklay(i+1)-aklay(i)
630
         if (abs(dh-dz).gt.eps) then
631
            print*,'Aklay: Vertical grid must be equally spaced... Stop'
632
            print*,(aklay(i),i=1,nz)
633
            stop
634
         endif
635
      enddo
636
 
637
c     Set missing data value
638
      mdv=misdat
639
 
640
      return
641
 
642
c     Exception handling
643
 998  print*,'Read_Dim: Problem with input netcdf file... Stop'
644
      stop
645
 
646
      end
647
 
648
 
649
 
650
c     --------------------------------------------------------------------------------
651
c     Calculate qgPV  from wind and theta 
652
c     --------------------------------------------------------------------------------       
653
 
654
      subroutine calc_qgpv (qgpv,uu,vv,th,
655
     >                      rhoref,pressref,nsqref,thetaref,coriol,
656
     >                      nx,ny,nz,deltax,deltay,deltaz,mdv)
657
 
658
c     Calculate qgPV from  wind and potential temperature according to
659
c     equation 2.9 p 16 Thesis Rene Fehlmann. Note a cartesian grid with
660
c     equidistant grid spacings <deltax,deltay,deltaz> is assumend. No
661
c     'correction' is made for spherical geoemtry.
662
 
663
      implicit none
664
 
665
c     Declaration of subroutine parameters
666
      integer   nx,ny,nz
667
      real      qgpv(0:nx,0:ny,0:nz)
668
      real      uu(0:nx,0:ny,0:nz)
669
      real      vv(0:nx,0:ny,0:nz)
670
      real      th(0:nx,0:ny,0:nz)
671
      real      rhoref(0:2*nz)
672
      real      nsqref(0:2*nz)
673
      real      thetaref(0:2*nz)
674
      real      pressref(0:2*nz)
675
      real      deltax,deltay,deltaz
676
      real      mdv
677
      real      coriol(0:nx,0:ny)
678
 
679
c     Numerical epsilon and physical constants
680
      real       g
681
      parameter  (g=9.81)
682
      real       eps
683
      parameter  (eps=0.01)
684
      real       scale
685
      parameter  (scale=1e6)
686
 
687
c     Auxiliary variables
688
      integer  i,j,k
689
      integer  kk,jj
690
      real     dvdx(0:nx,0:ny,0:nz)
691
      real     dudy(0:nx,0:ny,0:nz)
692
      real     dtdz(0:nx,0:ny,0:nz)
693
      real     t1,t2
694
 
695
c     Calculate horizontal derivatives dudy and dvdx of velocity
696
      do k=0,nz
697
         do j=0,ny
698
            do i=0,nx
699
 
700
c              Calculate dudy
701
               if (j.eq.0) then
702
                  if ( (abs(uu(i,1,k)-mdv).gt.eps).and.
703
     >                 (abs(uu(i,0,k)-mdv).gt.eps)) then
704
                     dudy(i,0,k)=(uu(i,1,k)-uu(i,0,k))/deltay
705
                  else
706
                     dudy(i,0,k)=mdv
707
                  endif
708
               elseif (j.eq.ny) then
709
                  if ( (abs(uu(i,ny,  k)-mdv).gt.eps).and.
710
     >                 (abs(uu(i,ny-1,k)-mdv).gt.eps)) then
711
                     dudy(i,ny,k)=(uu(i,ny,k)-uu(i,ny-1,k))/deltay
712
                  else
713
                     dudy(i,ny,k)=mdv
714
                  endif
715
               else
716
                  if ( (abs(uu(i,j+1,k)-mdv).gt.eps).and.
717
     >                 (abs(uu(i,j-1,k)-mdv).gt.eps)) then
718
                     dudy(i,j,k)=(uu(i,j+1,k)-uu(i,j-1,k))/(2.*deltay)
719
                  else
720
                     dudy(i,j,k)=mdv
721
                  endif
722
               endif
723
 
724
c              Calculate dvdx
725
               if (i.eq.0) then
726
                  if ((abs(vv(1,j,k)-mdv).gt.eps).and.
727
     >                (abs(vv(0,j,k)-mdv).gt.eps)) then
728
                     dvdx(0,j,k)=(vv(1,j,k)-vv(0,j,k))/deltax
729
                  else
730
                     dvdx(0,j,k)=mdv
731
                  endif
732
               elseif (i.eq.nx) then
733
                  if ((abs(vv(nx,  j,k)-mdv).gt.eps).and.
734
     >                (abs(vv(nx-1,j,k)-mdv).gt.eps)) then
735
                     dvdx(nx,j,k)=(vv(nx,j,k)-vv(nx-1,j,k))/deltax
736
                  else
737
                     dvdx(nx,j,k)=mdv
738
                  endif
739
               else
740
                  if ((abs(vv(i+1,j,k)-mdv).gt.eps).and.
741
     >                (abs(vv(i-1,j,k)-mdv).gt.eps)) then
742
                     dvdx(i,j,k)=(vv(i+1,j,k)-vv(i-1,j,k))/(2.*deltax)
743
                  else
744
                     dvdx(i,j,k)=mdv
745
                  endif
746
               endif
747
 
748
            enddo               
749
         enddo
750
      enddo
751
 
752
c     Calculate vertical derivative of potential temperature
753
      do i=0,nx
754
         do j=0,ny
755
            do k=0,nz
756
 
757
               if (k.eq.0) then
758
                  if ((abs(th(i,j,2)-mdv).gt.eps).and.
759
     >                (abs(th(i,j,1)-mdv).gt.eps)) then
760
                     t1=rhoref(2)*th(i,j,1)/thetaref(2)/nsqref(2)*g
761
                     t2=rhoref(0)*th(i,j,0)/thetaref(0)/nsqref(0)*g
762
                     dtdz(i,j,0)=(t1-t2)/deltaz
763
                  else
764
                     dtdz(i,j,0)=mdv
765
                  endif
766
              else if (k.eq.nz) then
767
                 if ((abs(th(i,j,nz  )-mdv).gt.eps).and.
768
     >               (abs(th(i,j,nz-1)-mdv).gt.eps)) then
769
                    kk=2*nz
770
                    t1=rhoref(kk  )*th(i,j,nz  )/
771
     >                     thetaref(kk)/nsqref(kk)*g
772
                    t2=rhoref(kk-2)*th(i,j,nz-1)/
773
     >                     thetaref(kk-2)/nsqref(kk-2)*9.8
774
                    dtdz(i,j,nz)=(t1-t2)/deltaz
775
                 else
776
                    dtdz(i,j,nz)=mdv
777
                 endif
778
              else
779
                 if ((abs(th(i,j,k+1)-mdv).gt.eps).and. 
780
     >               (abs(th(i,j,k  )-mdv).gt.eps).and.
781
     >               (abs(th(i,j,k-1)-mdv).gt.eps)) then
782
                    kk=2*k
783
                    t1=rhoref(kk+1)*(th(i,j,k)+th(i,j,k+1))/2.
784
     >                 /thetaref(kk+1)/nsqref(kk+1)*g
785
                    t2=rhoref(kk-1)*(th(i,j,k)+th(i,j,k-1))/2.
786
     >                 /thetaref(kk-1)/nsqref(kk-1)*g
787
                    dtdz(i,j,k)=(t1-t2)/deltaz
788
                 else
789
                    dtdz(i,j,k)=mdv
790
                 endif
791
              endif
792
 
793
           enddo
794
        enddo
795
      enddo
796
 
797
c     Calculate qgPV
798
      do i=0,nx
799
         do j=0,ny
800
            do k=0,nz
801
 
802
               kk=2*k
803
 
804
               if ((abs(dudy(i,j,k)-mdv).gt.eps).and. 
805
     >             (abs(dvdx(i,j,k)-mdv).gt.eps).and.
806
     >             (abs(dtdz(i,j,k)-mdv).gt.eps)) then
807
 
808
                  qgpv(i,j,k)=-dudy(i,j,k)+dvdx(i,j,k)+ 
809
     >                        coriol(i,j)*dtdz(i,j,k)/rhoref(kk)
810
               else
811
                  qgpv(i,j,k)=mdv
812
               endif
813
 
814
            enddo
815
         enddo
816
      enddo
817
 
818
      end
819