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 parameters
10
c     ---------------------------------------------------------------
11
 
12
      real      tzero
13
      parameter (tzero=273.15)
14
      real      kappa
15
      parameter (kappa=0.6078)
16
      real      rd
17
      parameter (rd=287.)
18
 
19
c     ---------------------------------------------------------------
20
c     Declaration of variables
21
c     ---------------------------------------------------------------
22
 
23
c     Variables for output Z file  : height level
24
      character*80 cfn
25
      real         varmin(4),varmax(4),stag(4)
26
      integer      vardim(4)
27
      real         mdv
28
      integer      ndim
29
      integer      nx,ny,nz
30
      real         xmin,xmax,ymin,ymax,dx,dy
31
      integer      ntimes
32
      real         aklev(1000),bklev(1000)
33
      real         aklay(1000),bklay(1000)
34
      real         time
35
      real         pollon,pollat
36
      integer      idate(5)
37
      integer      nfields
38
      character*80 field_nam(100)
39
      real,allocatable, dimension (:,:,:,:) :: field_dat
40
      real,allocatable, dimension (:,:,:)   :: z3
41
      real,allocatable, dimension (:,:)     :: x2,y2,f2,oro
42
      real,allocatable, dimension (:,:,:)   :: out1,out2
43
      real,allocatable, dimension (:,:,:)   :: inp
44
      real,allocatable,dimension (:)        :: nsqref
45
      real,allocatable,dimension (:)        :: thetaref
46
      real,allocatable,dimension (:)        :: tref
47
      real,allocatable,dimension (:)        :: rhoref
48
      real,allocatable,dimension (:)        :: pressref
49
      real,allocatable,dimension (:)        :: zref
50
      real         deltax,deltay,deltaz
51
      integer      nvars
52
      character*80 vnam(100),varname
53
      integer      isok
54
      integer      cdfid,cstid
55
      character*3  mode
56
 
57
c     Parameter file
58
      character*80 fieldname
59
      integer      nfilt
60
      character*80 ofn,gri
61
 
62
c     Auxiliary variables
63
      integer      ierr,stat
64
      integer      i,j,k,n
65
      real,allocatable, dimension (:,:)   :: tmp
66
      character*80 vnam2(100)
67
      integer      nvars2
68
 
69
c     ---------------------------------------------------------------
70
c     Preparations
71
c     ---------------------------------------------------------------
72
 
73
      print*,'*********************************************************'
74
      print*,'* qvec_analysis                                         *'
75
      print*,'*********************************************************'
76
 
77
c     Read parameter file
78
      open(10,file='fort.10')
79
       read(10,*) fieldname
80
       read(10,*) ofn
81
       read(10,*) gri
82
       read(10,*) nfilt
83
      close(10)
84
 
85
c     Decide whether to enter ANO or MOD/ORG mode
86
      mode = ofn(1:3)
87
      if ( (mode.ne.'ANO').and.
88
     >     (mode.ne.'MOD').and.
89
     >     (mode.ne.'ORG') )
90
     >then
91
         print*,'Unsupported mode ',trim(mode)
92
         stop
93
      endif
94
 
95
c     Get grid description for Z file : height level
96
      call cdfopn(ofn,cdfid,ierr)
97
      if (ierr.ne.0) goto 998
98
      call getcfn(cdfid,cfn,ierr)
99
      if (ierr.ne.0) goto 998
100
      call cdfopn(cfn,cstid,ierr)
101
      if (ierr.ne.0) goto 998
102
      call getdef(cdfid,'T',ndim,mdv,vardim,
103
     >            varmin,varmax,stag,ierr)
104
      if (ierr.ne.0) goto 998
105
      nx  =vardim(1)
106
      ny  =vardim(2)
107
      nz  =vardim(3)
108
      xmin=varmin(1)
109
      ymin=varmin(2)
110
      call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)
111
      if (ierr.ne.0) goto 998
112
      call getgrid(cstid,dx,dy,ierr)
113
      if (ierr.ne.0) goto 998
114
      xmax=xmin+real(nx-1)*dx
115
      ymax=ymin+real(ny-1)*dy
116
      call gettimes(cdfid,time,ntimes,ierr)
117
      if (ierr.ne.0) goto 998
118
      call getstart(cstid,idate,ierr)
119
      if (ierr.ne.0) goto 998
120
      call getpole(cstid,pollon,pollat,ierr)
121
      if (ierr.ne.0) goto 998
122
      call getvars(cdfid,nvars,vnam,ierr)
123
      if (ierr.ne.0) goto 998
124
      call clscdf(cstid,ierr)
125
      if (ierr.ne.0) goto 998
126
      call clscdf(cdfid,ierr)
127
      if (ierr.ne.0) goto 998
128
 
129
c     Get a list of all variables on the GRID file
130
      if (trim(gri).ne.trim(ofn)) then
131
         call cdfopn(gri,cdfid,ierr)
132
         if (ierr.ne.0) goto 998
133
         call getvars(cdfid,nvars2,vnam2,ierr)
134
         if (ierr.ne.0) goto 998
135
         do i=1,nvars2
136
            vnam(nvars+i)=vnam2(i)
137
         enddo
138
         nvars=nvars+nvars2
139
      endif
140
 
141
c     Check whether calculation of <fieldname> is supported
142
      if ( (fieldname.ne.'GEO' ).and.
143
     >     (fieldname.ne.'AGEO' ).and.
144
     >     (fieldname.ne.'DIV_UV' ).and.
145
     >     (fieldname.ne.'QVEC' ).and.
146
     >     (fieldname.ne.'DIV_Q' ) ) then
147
         print*,'Variable not supported ',trim(fieldname)
148
         stop
149
      endif
150
 
151
c     Set dependencies
152
      if (fieldname.eq.'GEO') then
153
         nfields=2
154
         field_nam(1)='T'
155
         field_nam(2)='P'
156
      elseif (fieldname.eq.'AGEO') then
157
         nfields=4
158
         field_nam(1)='U'
159
         field_nam(2)='UG'
160
         field_nam(3)='V'
161
         field_nam(4)='VG'
162
      else if (fieldname.eq.'DIV_UV') then
163
         nfields=2
164
         field_nam(1)='U'
165
         field_nam(2)='V'     
166
      else if (fieldname.eq.'QVEC') then
167
         nfields=3
168
         field_nam(1)='UG'
169
         field_nam(2)='VG'
170
         field_nam(3)='TH'
171
      else if (fieldname.eq.'DIV_Q') then
172
         nfields=2
173
         field_nam(1)='QX'
174
         field_nam(2)='QY'     
175
      endif
176
 
177
c     Allocate memory
178
      allocate(field_dat(nfields,nx,ny,nz),stat=stat)
179
      if (stat.ne.0) print*,'error allocating field_dat'
180
      allocate(out1(nx,ny,nz),stat=stat)
181
      if (stat.ne.0) print*,'error allocating out1'
182
      allocate(out2(nx,ny,nz),stat=stat)
183
      if (stat.ne.0) print*,'error allocating out2'
184
      allocate(inp(nx,ny,nz),stat=stat)
185
      if (stat.ne.0) print*,'error allocating inp'
186
      allocate(z3(nx,ny,nz),stat=stat)
187
      if (stat.ne.0) print*,'error allocating z3'
188
      allocate(x2(nx,ny),stat=stat)
189
      if (stat.ne.0) print*,'error allocating x2'
190
      allocate(y2(nx,ny),stat=stat)
191
      if (stat.ne.0) print*,'error allocating y2'
192
      allocate(f2(nx,ny),stat=stat)
193
      if (stat.ne.0) print*,'error allocating f2'
194
      allocate(oro(nx,ny),stat=stat)
195
      if (stat.ne.0) print*,'error allocating oro'
196
 
197
C     Allocate memory for reference profile
198
      allocate(rhoref  (0:2*nz),STAT=stat)
199
      if (stat.ne.0) print*,'error allocating rhoref'
200
      allocate(pressref(0:2*nz),STAT=stat)
201
      if (stat.ne.0) print*,'error allocating pressref'
202
      allocate(thetaref(0:2*nz),STAT=stat)
203
      if (stat.ne.0) print*,'error allocating thetaref'
204
      allocate(nsqref  (0:2*nz),STAT=stat)
205
      if (stat.ne.0) print*,'error allocating nsqref'
206
      allocate(zref    (0:2*nz),STAT=stat)
207
      if (stat.ne.0) print*,'error allocating zref'
208
      allocate(tref    (0:2*nz),STAT=stat)
209
      if (stat.ne.0) print*,'error allocating tref'
210
 
211
c     Allocate auxiliary fields
212
      allocate(tmp(nx,ny),stat=stat)
213
      if (stat.ne.0) print*,'error allocating tmp'
214
 
215
c     Read reference profile and grid parameters
216
      call read_ref (nsqref,rhoref,thetaref,pressref,zref,
217
     >               nx,ny,nz,deltax,deltay,deltaz,f2,oro,gri)
218
 
219
c     Read X grid
220
      varname='X'
221
      isok=0
222
      call check_varok (isok,varname,vnam,nvars)
223
      if (isok.eq.0) then
224
         print*,'Variable ',trim(varname),' missing... Stop'
225
         stop
226
      endif
227
      call cdfopn(gri,cdfid,ierr)
228
      if (ierr.ne.0) goto 998
229
      call getdat(cdfid,varname,time,0,x2,ierr)
230
      if (ierr.ne.0) goto 998
231
      print*,'R ',trim(varname),' ',trim(gri)
232
      call clscdf(cdfid,ierr)
233
      if (ierr.ne.0) goto 998
234
 
235
c     Read Y grid
236
      varname='Y'
237
      isok=0
238
      call check_varok (isok,varname,vnam,nvars)
239
      if (isok.eq.0) then
240
         print*,'Variable ',trim(varname),' missing... Stop'
241
         stop
242
      endif
243
      call cdfopn(gri,cdfid,ierr)
244
      if (ierr.ne.0) goto 998
245
      call getdat(cdfid,varname,time,0,y2,ierr)
246
      if (ierr.ne.0) goto 998
247
      print*,'R ',trim(varname),' ',trim(gri)
248
      call clscdf(cdfid,ierr)
249
      if (ierr.ne.0) goto 998
250
 
251
c     Calculate reference profile of temperature
252
      print*,'C T_ref = TH_ref * ( P_ref/1000)^kappa'
253
      do i=0,2*nz
254
         tref(i) = thetaref(i) * ( pressref(i)/1000. ) ** kappa
255
      enddo
256
 
257
c     Init height levels
258
      do i=1,nx
259
         do j=1,ny
260
            do k=1,nz
261
               z3(i,j,k)=zref(2*k-1)
262
            enddo
263
         enddo
264
      enddo
265
 
266
c     Load needed variables
267
      do n=1,nfields
268
 
269
c        Check whether variable is available on file
270
         varname=field_nam(n)
271
         isok=0
272
         call check_varok (isok,varname,vnam,nvars)
273
 
274
c        Variable is available on file
275
         if (isok.eq.1) then
276
 
277
            call cdfopn(ofn,cdfid,ierr)
278
            if (ierr.ne.0) goto 998
279
            call getdat(cdfid,varname,time,0,inp,ierr)
280
            if (ierr.ne.0) goto 998
281
            print*,'R ',trim(varname),' ',trim(ofn)
282
            call clscdf(cdfid,ierr)
283
            if (ierr.ne.0) goto 998
284
 
285
            do i=1,nx
286
               do j=1,ny
287
                  do k=1,nz
288
                     field_dat(n,i,j,k)=inp(i,j,k)
289
                  enddo
290
               enddo
291
            enddo
292
 
293
         else            
294
            print*,'Variable missing : ',trim(varname)
295
            stop
296
         endif
297
 
298
      enddo
299
 
300
c     Add reference profile if ANO file is provided
301
      if ( mode.eq.'ANO' ) then
302
 
303
          print*,'C T -> T_ano + T_ref'
304
          n=0
305
          do i=1,nfields
306
              if (field_nam(i).eq.'T') n=i
307
          enddo
308
          if (n.ne.0) then
309
           do i=1,nx
310
            do j=1,ny
311
              do k=1,nz
312
                 field_dat(n,i,j,k) = field_dat(n,i,j,k) + tref(2*k-1)
313
              enddo
314
            enddo
315
           enddo
316
          endif
317
 
318
          print*,'C TH -> TH_ano + TH_ref'
319
          n=0
320
          do i=1,nfields
321
              if (field_nam(i).eq.'TH') n=i
322
          enddo
323
          if (n.ne.0) then
324
           do i=1,nx
325
            do j=1,ny
326
              do k=1,nz
327
                 field_dat(n,i,j,k) = field_dat(n,i,j,k)+thetaref(2*k-1)
328
              enddo
329
            enddo
330
           enddo
331
          endif
332
 
333
          print*,'C P -> P_ano + P_ref'
334
          n=0
335
          do i=1,nfields
336
              if (field_nam(i).eq.'P') n=i
337
          enddo
338
          if (n.ne.0) then
339
           do i=1,nx
340
            do j=1,ny
341
              do k=1,nz
342
                 field_dat(n,i,j,k) = field_dat(n,i,j,k)+pressref(2*k-1)
343
              enddo
344
            enddo
345
           enddo
346
          endif
347
 
348
      endif
349
 
350
c     Change units: P (hPa -> Pa), T(K -> degC)
351
      n=0
352
      do i=1,nfields
353
         if (field_nam(i).eq.'P') n=i
354
      enddo
355
      if (n.ne.0) then
356
         do i=1,nx
357
            do j=1,ny
358
               do k=1,nz
359
                  field_dat(n,i,j,k)=100.*field_dat(n,i,j,k)
360
               enddo
361
            enddo
362
         enddo
363
      endif
364
 
365
      n=0
366
      do i=1,nfields
367
         if (field_nam(i).eq.'T') n=i
368
      enddo
369
      if (n.ne.0) then
370
         do i=1,nx
371
            do j=1,ny
372
               do k=1,nz
373
                  if ( field_dat(n,i,j,1).gt.100. ) then
374
                    field_dat(n,i,j,k)=field_dat(n,i,j,k) - tzero
375
                  endif
376
               enddo
377
            enddo
378
         enddo
379
      endif
380
 
381
c     ----------------------------------------------------------------
382
c     Calculations
383
c     ----------------------------------------------------------------
384
 
385
c     Geostrophic wind (GEO)
386
      if (fieldname.eq.'GEO') then
387
 
388
         print*,'C ',trim(fieldname)
389
         field_nam(1)='RHO'
390
         do i=1,nx
391
            do j=1,ny
392
               do k=1,nz
393
                  if ( mode.eq.'ANO' ) then
394
                    field_dat(1,i,j,k)=rhoref(2*k-1)
395
                  else
396
                    field_dat(1,i,j,k)=
397
     >               1./rd*field_dat(2,i,j,k)/(field_dat(1,i,j,k)+tzero)
398
                  endif
399
                enddo
400
            enddo
401
         enddo
402
         call calc_geo (out1,out2,                           ! (Ug,Vg)
403
     >                  field_dat(1,:,:,:),                  ! RHO
404
     >                  field_dat(2,:,:,:),                  ! P
405
     >                  z3,x2,y2,f2,                         ! Z,X,Y,CORIOL
406
     >                  nx,ny,nz,mdv)
407
 
408
c     Ageostrophic wind (AGEO)
409
      elseif (fieldname.eq.'AGEO') then
410
 
411
         print*,'C ',trim(fieldname)
412
         call calc_ageo (out1,out2,                           ! (Ua,Va)
413
     >                   field_dat(1,:,:,:),                  ! U
414
     >                   field_dat(2,:,:,:),                  ! UG
415
     >                   field_dat(3,:,:,:),                  ! V
416
     >                   field_dat(4,:,:,:),                  ! VG
417
     >                   nx,ny,nz,mdv)
418
 
419
c     Divergence of wind (DIV_UV)
420
      else if (fieldname.eq.'DIV_UV') then
421
 
422
         print*,'C ',trim(fieldname)
423
         call calc_div_uv (out1,                             ! Div(U,V))
424
     >                     field_dat(1,:,:,:),               ! U
425
     >                     field_dat(2,:,:,:),               ! V
426
     >                     z3,x2,y2,f2,                      ! Z,X,Y,CORIOL
427
     >                     nx,ny,nz,mdv)
428
 
429
c     Q vector (QVEC)
430
      else if (fieldname.eq.'QVEC') then
431
 
432
         print*,'C ',trim(fieldname)
433
         call calc_qvec (out1,out2,                           ! (Qx,Qy)
434
     >                   field_dat(1,:,:,:),                  ! UG
435
     >                   field_dat(2,:,:,:),                  ! VG
436
     >                   field_dat(3,:,:,:),                  ! TH
437
     >                   z3,x2,y2,f2,                         ! Z,X,Y,CORIOL
438
     >                   nx,ny,nz,mdv)
439
 
440
c     Divergence of Q vector (DIV_Q)
441
      else if (fieldname.eq.'DIV_Q') then
442
 
443
         print*,'C ',trim(fieldname)
444
         call calc_div_q  (out1,                             ! Div(U,V))
445
     >                     field_dat(1,:,:,:),               ! QX
446
     >                     field_dat(2,:,:,:),               ! QY
447
     >                     z3,x2,y2,f2,                      ! Z,X,Y,CORIOL
448
     >                     nx,ny,nz,mdv)
449
 
450
      endif
451
 
452
c     Horizontal filtering of the resulting fields
453
      print*,'F ',trim(fieldname)
454
      do k=1,nz
455
 
456
          do i=1,nx
457
             do j=1,ny
458
                tmp(i,j)=out1(i,j,k)
459
             enddo
460
          enddo
461
          do n=1,nfilt     
462
             call filt2d (tmp,tmp,nx,ny,1.,mdv,0,0,0,0)
463
          enddo
464
          do i=1,nx
465
             do j=1,ny
466
                out1(i,j,k)=tmp(i,j)
467
             enddo
468
          enddo
469
 
470
          do i=1,nx
471
             do j=1,ny
472
                tmp(i,j)=out2(i,j,k)
473
             enddo
474
          enddo
475
          do n=1,nfilt     
476
             call filt2d (tmp,tmp,nx,ny,1.,mdv,0,0,0,0)
477
          enddo
478
          do i=1,nx
479
             do j=1,ny
480
                out2(i,j,k)=tmp(i,j)
481
             enddo
482
          enddo
483
 
484
       enddo
485
 
486
c     ----------------------------------------------------------------
487
c     Save result onto netcdf file
488
c     ----------------------------------------------------------------
489
 
490
c     Open output file
491
      call cdfwopn(ofn,cdfid,ierr)
492
      if (ierr.ne.0) goto 998
493
 
494
c     Save geostrophic wind
495
      if (fieldname.eq.'GEO') then
496
         isok=0
497
         varname='UG'
498
         call check_varok(isok,varname,vnam,nvars)
499
         if (isok.eq.0) then
500
            call putdef(cdfid,varname,ndim,mdv,vardim,
501
     >                  varmin,varmax,stag,ierr)
502
            if (ierr.ne.0) goto 998
503
         endif
504
         call putdat(cdfid,varname,time,0,out1,ierr)
505
         if (ierr.ne.0) goto 998
506
         print*,'W ',trim(varname),' ',trim(ofn)
507
                  isok=0
508
         varname='VG'
509
         call check_varok(isok,varname,vnam,nvars)
510
         if (isok.eq.0) then
511
            call putdef(cdfid,varname,ndim,mdv,vardim,
512
     >                  varmin,varmax,stag,ierr)
513
            if (ierr.ne.0) goto 998
514
         endif
515
         call putdat(cdfid,varname,time,0,out2,ierr)
516
         if (ierr.ne.0) goto 998
517
         print*,'W ',trim(varname),' ',trim(ofn)
518
 
519
c     Save ageostrophic wind
520
      elseif (fieldname.eq.'AGEO') then
521
         isok=0
522
         varname='UA'
523
         call check_varok(isok,varname,vnam,nvars)
524
         if (isok.eq.0) then
525
            call putdef(cdfid,varname,ndim,mdv,vardim,
526
     >                  varmin,varmax,stag,ierr)
527
            if (ierr.ne.0) goto 998
528
         endif
529
         call putdat(cdfid,varname,time,0,out1,ierr)
530
         if (ierr.ne.0) goto 998
531
         print*,'W ',trim(varname),' ',trim(ofn)
532
                  isok=0
533
         varname='VA'
534
         call check_varok(isok,varname,vnam,nvars)
535
         if (isok.eq.0) then
536
            call putdef(cdfid,varname,ndim,mdv,vardim,
537
     >                  varmin,varmax,stag,ierr)
538
            if (ierr.ne.0) goto 998
539
         endif
540
         call putdat(cdfid,varname,time,0,out2,ierr)
541
         if (ierr.ne.0) goto 998
542
         print*,'W ',trim(varname),' ',trim(ofn)
543
 
544
c     Save divergence of wind field
545
      else if (fieldname.eq.'DIV_UV') then
546
         isok=0
547
         varname='DIV_UV'
548
         call check_varok(isok,varname,vnam,nvars)
549
         if (isok.eq.0) then
550
            call putdef(cdfid,varname,ndim,mdv,vardim,
551
     >                  varmin,varmax,stag,ierr)
552
            if (ierr.ne.0) goto 998
553
         endif
554
         call putdat(cdfid,varname,time,0,out1,ierr)
555
         if (ierr.ne.0) goto 998
556
         print*,'W ',trim(varname),' ',trim(ofn)
557
 
558
c     Save components of Q vector
559
      else if (fieldname.eq.'QVEC') then
560
         isok=0
561
         varname='QX'
562
         call check_varok(isok,varname,vnam,nvars)
563
         if (isok.eq.0) then
564
            call putdef(cdfid,varname,ndim,mdv,vardim,
565
     >                  varmin,varmax,stag,ierr)
566
            if (ierr.ne.0) goto 998
567
         endif
568
         call putdat(cdfid,varname,time,0,out1,ierr)
569
         if (ierr.ne.0) goto 998
570
         print*,'W ',trim(varname),' ',trim(ofn)
571
                  isok=0
572
         varname='QY'
573
         call check_varok(isok,varname,vnam,nvars)
574
         if (isok.eq.0) then
575
            call putdef(cdfid,varname,ndim,mdv,vardim,
576
     >                  varmin,varmax,stag,ierr)
577
            if (ierr.ne.0) goto 998
578
         endif
579
         call putdat(cdfid,varname,time,0,out2,ierr)
580
         if (ierr.ne.0) goto 998
581
         print*,'W ',trim(varname),' ',trim(ofn)
582
 
583
c     Save divergence of wind field
584
      else if (fieldname.eq.'DIV_Q') then
585
         isok=0
586
         varname='DIV_Q'
587
         call check_varok(isok,varname,vnam,nvars)
588
         if (isok.eq.0) then
589
            call putdef(cdfid,varname,ndim,mdv,vardim,
590
     >                  varmin,varmax,stag,ierr)
591
            if (ierr.ne.0) goto 998
592
         endif
593
         call putdat(cdfid,varname,time,0,out1,ierr)
594
         if (ierr.ne.0) goto 998
595
         print*,'W ',trim(varname),' ',trim(ofn)
596
 
597
      endif
598
 
599
c     Close output file
600
      call clscdf(cdfid,ierr)
601
      if (ierr.ne.0) goto 998
602
 
603
 
604
c     ----------------------------------------------------------------
605
c     Exception handling
606
c     ----------------------------------------------------------------
607
 
608
      stop
609
 
610
 998  print*,'Problem with netcdf file'
611
      stop
612
 
613
      end
614
 
615
 
616
 
617
c     ****************************************************************
618
c     * SUBROUTINE SECTION: AUXILIARY ROUTINES                       *
619
c     ****************************************************************     
620
 
621
c     ----------------------------------------------------------------
622
c     Check whether variable is found on netcdf file
623
c     ----------------------------------------------------------------
624
 
625
      subroutine check_varok (isok,varname,varlist,nvars)
626
 
627
c     Check whether the variable <varname> is in the list <varlist(nvars)>. 
628
c     If this is the case, <isok> is incremented by 1. Otherwise <isok> 
629
c     keeps its value.
630
 
631
      implicit none
632
 
633
c     Declaraion of subroutine parameters
634
      integer      isok
635
      integer      nvars
636
      character*80 varname
637
      character*80 varlist(nvars)
638
 
639
c     Auxiliary variables
640
      integer      i
641
 
642
c     Main
643
      do i=1,nvars
644
         if (trim(varname).eq.trim(varlist(i))) isok=isok+1
645
      enddo
646
 
647
      end
648
 
649
 
650
c     ****************************************************************
651
c     * SUBROUTINE SECTION: CALCULATE SECONDARY FIELDS               *
652
c     ****************************************************************
653
 
654
c     ----------------------------------------------------------------
655
c     Calculate geostrophic wind components
656
c     ----------------------------------------------------------------
657
 
658
      subroutine calc_geo (ug,vg,rho,p,
659
     >                     z3,x2,y2,f2,nx,ny,nz,mdv) 
660
 
661
c     Calculate the geostrophic wind components (ug,vg) if the temperature 
662
c     (t) and the pressure (p) are given. The grid and the missing data
663
c     value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv) 
664
 
665
      implicit none
666
 
667
c     Declaration of subroutine parameters
668
      integer   nx,ny,nz
669
      real      ug (nx,ny,nz)
670
      real      vg (nx,ny,nz)
671
      real      rho(nx,ny,nz)
672
      real      p  (nx,ny,nz)
673
      real      z3 (nx,ny,nz)
674
      real      x2 (nx,ny)
675
      real      y2 (nx,ny)
676
      real      f2 (nx,ny)
677
      real      mdv
678
 
679
c     Physical parameters and numerical constants
680
      real      eps
681
      parameter (eps=0.01)
682
      real       g
683
      parameter  (g=9.80616)
684
 
685
 
686
c     Auxiliray variables
687
      integer   i,j,k
688
      real      dpdx(nx,ny,nz)
689
      real      dpdy(nx,ny,nz)
690
 
691
c     Calculate horizontal derivatives of pressure
692
      call deriv(dpdx,p,'x',z3,x2,y2,nx,ny,nz,mdv)
693
      call deriv(dpdy,p,'y',z3,x2,y2,nx,ny,nz,mdv)
694
 
695
c     Calculation
696
      do i=1,nx
697
         do j=1,ny
698
            do k=1,nz
699
 
700
               if ((abs(rho (i,j,k)-mdv).gt.eps).and.
701
     >             (abs(dpdx(i,j,k)-mdv).gt.eps).and.
702
     >             (abs(dpdy(i,j,k)-mdv).gt.eps)) then
703
 
704
                  ug(i,j,k)=-1./(rho(i,j,k)*f2(i,j))*dpdy(i,j,k)
705
                  vg(i,j,k)= 1./(rho(i,j,k)*f2(i,j))*dpdx(i,j,k)
706
 
707
               endif
708
 
709
            enddo
710
         enddo
711
      enddo
712
 
713
      end
714
 
715
c     ----------------------------------------------------------------
716
c     Calculate ageostrophic wind components
717
c     ----------------------------------------------------------------
718
 
719
      subroutine calc_ageo (ua,va,u,ug,v,vg,nx,ny,nz,mdv)
720
 
721
c     Calculate the geostrophic wind components (ug,vg) if the temperature
722
c     (t) and the pressure (p) are given. The grid and the missing data
723
c     value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv)
724
 
725
      implicit none
726
 
727
c     Declaration of subroutine parameters
728
      integer   nx,ny,nz
729
      real      ua (nx,ny,nz)
730
      real      va (nx,ny,nz)
731
      real      u  (nx,ny,nz)
732
      real      ug (nx,ny,nz)
733
      real      v  (nx,ny,nz)
734
      real      vg (nx,ny,nz)
735
      real      mdv
736
 
737
c     Parameters
738
      real      eps
739
      parameter (eps=0.01)
740
 
741
c     Auxiliray variables
742
      integer   i,j,k
743
 
744
c     Calculation
745
      do i=1,nx
746
         do j=1,ny
747
            do k=1,nz
748
 
749
               if ((abs(u (i,j,k)-mdv).gt.eps).and.
750
     >             (abs(ug(i,j,k)-mdv).gt.eps).and.
751
     >             (abs(v (i,j,k)-mdv).gt.eps).and.
752
     >             (abs(vg(i,j,k)-mdv).gt.eps)) then
753
 
754
                  ua(i,j,k)=u(i,j,k) - ug(i,j,k)
755
                  va(i,j,k)=v(i,j,k) - vg(i,j,k)
756
 
757
               endif
758
 
759
            enddo
760
         enddo
761
      enddo
762
 
763
      end
764
 
765
c     ----------------------------------------------------------------
766
c     Calculate divergence of wind field
767
c     ----------------------------------------------------------------
768
 
769
      subroutine calc_div_uv (div_uv,u,v,
770
     >                        z3,x2,y2,f2,nx,ny,nz,mdv) 
771
 
772
c     Calculate the divergence (div_uv) of the horizontal wind field if+
773
c     the wind components (u,v) are given. The grid and the missing data
774
c     value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv) 
775
 
776
      implicit none
777
 
778
c     Declaration of subroutine parameters
779
      integer   nx,ny,nz
780
      real      div_uv(nx,ny,nz)
781
      real      u     (nx,ny,nz)
782
      real      v     (nx,ny,nz)
783
      real      z3    (nx,ny,nz)
784
      real      x2    (nx,ny)
785
      real      y2    (nx,ny)
786
      real      f2    (nx,ny)
787
      real      mdv
788
 
789
c     Physical parameters and numerical constants
790
      real      eps
791
      parameter (eps=0.01)
792
 
793
c     Auxiliray variables
794
      integer   i,j,k
795
      real      rho
796
      real      dudx(nx,ny,nz)
797
      real      dvdy(nx,ny,nz)
798
 
799
c     Calculate horizontal derivatives of pressure
800
      call deriv(dudx,u,'x',z3,x2,y2,nx,ny,nz,mdv)
801
      call deriv(dvdy,v,'y',z3,x2,y2,nx,ny,nz,mdv)
802
 
803
c     Calculation
804
      do i=1,nx
805
         do j=1,ny
806
            do k=1,nz
807
 
808
               if ((abs(dudx(i,j,k)-mdv).gt.eps).and.
809
     >             (abs(dvdy(i,j,k)-mdv).gt.eps)) then
810
 
811
                  div_uv(i,j,k)=dudx(i,j,k)+dvdy(i,j,k)
812
 
813
               endif
814
 
815
            enddo
816
         enddo
817
      enddo
818
 
819
      end
820
 
821
c     ----------------------------------------------------------------
822
c     Calculate Q vector
823
c     ----------------------------------------------------------------
824
 
825
      subroutine calc_qvec (qx3,qy3,
826
     >                     th3,u3,v3,
827
     >                     z3,x2,y2,f2,nx,ny,nz,mdv) 
828
 
829
c     Calculate teh Q vector components <qx3> and <qy3>, as well as the divergence
830
c     of the Q vector <divq3> on the model grid. The grid is specified in the horizontal
831
c     by <xmin,ymin,dx,dy,nx,ny>. The number of vertical levels is <nz>. The input field
832
c     are: potential temperature <th3>, horizontal wind <u3> and <v3>, pressure <p3>.
833
c     The calculation follows the one described in "Weather Analysis, Dusan Djuric"
834
 
835
      implicit none
836
 
837
c     Declaration of subroutine parameters
838
      integer   nx,ny,nz
839
      real      qx3(nx,ny,nz)
840
      real      qy3(nx,ny,nz)
841
      real      th3(nx,ny,nz)
842
      real      u3 (nx,ny,nz)
843
      real      v3 (nx,ny,nz)
844
      real      z3 (nx,ny,nz)
845
      real      x2 (nx,ny)
846
      real      y2 (nx,ny)
847
      real      f2 (nx,ny)
848
      real      mdv
849
 
850
c     Physical and numerical parameters
851
      real      scale1,scale2
852
      parameter (scale1=1.E10,scale2=1.E14)
853
      real      eps
854
      parameter (eps=0.01)
855
      real      g
856
      parameter (g=9.80616)
857
      real      tzero
858
      parameter (tzero=273.16)
859
 
860
c     Auxiliary variables
861
      real      dudx(nx,ny,nz)
862
      real      dudy(nx,ny,nz)
863
      real      dvdx(nx,ny,nz)
864
      real      dvdy(nx,ny,nz)
865
      real      dtdx(nx,ny,nz)
866
      real      dtdy(nx,ny,nz)
867
      integer   i,j,k
868
 
869
c     Needed derivatives
870
      call deriv (dudx, u3,'x',z3,x2,y2,nx,ny,nz,mdv)
871
      call deriv (dudy, u3,'y',z3,x2,y2,nx,ny,nz,mdv)
872
      call deriv (dvdx, v3,'x',z3,x2,y2,nx,ny,nz,mdv)
873
      call deriv (dvdy, v3,'y',z3,x2,y2,nx,ny,nz,mdv)
874
      call deriv (dtdx,th3,'x',z3,x2,y2,nx,ny,nz,mdv)
875
      call deriv (dtdy,th3,'y',z3,x2,y2,nx,ny,nz,mdv)
876
 
877
c     Calculate vector components of the Q vector
878
      do i=1,nx
879
         do j=1,ny
880
            do k=1,nz
881
 
882
c              Evaluate Q vector formula with missing data check
883
               if ((abs(dudx(i,j,k)-mdv).gt.eps).and.
884
     >             (abs(dudy(i,j,k)-mdv).gt.eps).and.
885
     >             (abs(dvdx(i,j,k)-mdv).gt.eps).and.
886
     >             (abs(dvdy(i,j,k)-mdv).gt.eps).and.
887
     >             (abs(dtdx(i,j,k)-mdv).gt.eps).and.
888
     >             (abs(dtdy(i,j,k)-mdv).gt.eps)) then
889
 
890
                  qx3(i,j,k) = -g/tzero * (dudx(i,j,k)*dtdx(i,j,k)
891
     >                                     +dvdx(i,j,k)*dtdy(i,j,k))
892
                  qy3(i,j,k) = -g/tzero * (dudy(i,j,k)*dtdx(i,j,k)
893
     >                                     +dvdy(i,j,k)*dtdy(i,j,k))
894
 
895
               else
896
                  qx3(i,j,k)=mdv
897
                  qy3(i,j,k)=mdv
898
               endif
899
 
900
            enddo
901
         enddo
902
      enddo
903
 
904
c     Scale the output
905
      do i=1,nx
906
         do j=1,ny
907
            do k=1,nz
908
               if (abs(qx3(i,j,k)-mdv).gt.eps) then
909
                  qx3(i,j,k)=scale1*qx3(i,j,k)
910
               endif
911
               if (abs(qy3(i,j,k)-mdv).gt.eps) then
912
                  qy3(i,j,k)=scale1*qy3(i,j,k)
913
               endif
914
            enddo
915
         enddo
916
      enddo
917
 
918
      end
919
 
920
c     ----------------------------------------------------------------
921
c     Calculate divergence of wind field
922
c     ----------------------------------------------------------------
923
 
924
      subroutine calc_div_q (div_q,qx,qy,
925
     >                       z3,x2,y2,f2,nx,ny,nz,mdv) 
926
 
927
c     Calculate the divergence (div_q) of the Q vector field if
928
c     the components (qx,qy) are given. The grid and the missing data
929
c     value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv) 
930
 
931
      implicit none
932
 
933
c     Declaration of subroutine parameters
934
      integer   nx,ny,nz
935
      real      div_q(nx,ny,nz)
936
      real      qx   (nx,ny,nz)
937
      real      qy   (nx,ny,nz)
938
      real      z3   (nx,ny,nz)
939
      real      x2   (nx,ny)
940
      real      y2   (nx,ny)
941
      real      f2   (nx,ny)
942
      real      mdv
943
 
944
c     Physical parameters and numerical constants
945
      real      eps
946
      parameter (eps=0.01)
947
 
948
c     Auxiliray variables
949
      integer   i,j,k
950
      real      rho
951
      real      dqxdx(nx,ny,nz)
952
      real      dqydy(nx,ny,nz)
953
 
954
c     Calculate horizontal derivatives of pressure
955
      call deriv(dqxdx,qx,'x',z3,x2,y2,nx,ny,nz,mdv)
956
      call deriv(dqydy,qy,'y',z3,x2,y2,nx,ny,nz,mdv)
957
 
958
c     Calculation
959
      do i=1,nx
960
         do j=1,ny
961
            do k=1,nz
962
 
963
               if ((abs(dqxdx(i,j,k)-mdv).gt.eps).and.
964
     >             (abs(dqydy(i,j,k)-mdv).gt.eps)) then
965
 
966
                  div_q(i,j,k)=dqxdx(i,j,k)+dqydy(i,j,k)
967
 
968
               endif
969
 
970
            enddo
971
         enddo
972
      enddo
973
 
974
      end
975
 
976
 
977
c     ****************************************************************
978
c     * SUBROUTINE SECTION: GRID HANDLING                            *
979
c     ****************************************************************
980
 
981
c     -----------------------------------------------------------------
982
c     Horizontal and vertical derivatives for 3d fields
983
c     -----------------------------------------------------------------
984
 
985
      subroutine deriv (df,f,direction,
986
     >                  z3,x2,y2,nx,ny,nz,mdv)
987
 
988
c     Calculate horizontal and vertical derivatives of the 3d field <f>.
989
c     The direction of the derivative is specified in <direction> 
990
c         'x','y'          : Horizontal derivative in x and y direction
991
c         'p','z','t','m'  : Vertical derivative (pressure, height, theta, model)
992
c     The 3d field <z3> specifies the isosurfaces along which the horizontal 
993
c     derivatives are calculated or the levels for the vertical derivatives.
994
 
995
      implicit none
996
 
997
c     Input and output parameters
998
      integer    nx,ny,nz
999
      real       df (nx,ny,nz)
1000
      real       f  (nx,ny,nz)
1001
      real       z3 (nx,ny,nz)
1002
      real       x2 (nx,ny)
1003
      real       y2 (nx,ny)
1004
      character  direction
1005
      real       mdv
1006
 
1007
c     Numerical and physical parameters
1008
      real       pi180
1009
      parameter  (pi180=3.141592654/180.)
1010
      real       deltay
1011
      parameter  (deltay=111.1775E3)
1012
      real       zerodiv
1013
      parameter  (zerodiv=0.00000001)
1014
      real       eps
1015
      parameter  (eps=0.01) 
1016
 
1017
c     Auxiliary variables
1018
      integer    i,j,k
1019
      real       vmin,vmax
1020
      real       scale,lat
1021
      real       vu,vl,vuvl,vlvu
1022
      integer    o,u,w,e,n,s
1023
 
1024
c     Vertical derivative
1025
      if ((direction.eq.'z').or.
1026
     >    (direction.eq.'th').or.
1027
     >    (direction.eq.'p').or.
1028
     >    (direction.eq.'m').and.
1029
     >    (nz.gt.1)) then
1030
 
1031
         do i=1,nx
1032
            do j=1,ny
1033
               do k=1,nz
1034
 
1035
                  o=k+1
1036
                  if (o.gt.nz) o=nz
1037
                  u=k-1
1038
                  if (u.lt.1) u=1                  
1039
 
1040
                  if ((abs(f(i,j,o)-mdv).gt.eps).and.
1041
     >                (abs(f(i,j,u)-mdv).gt.eps).and.
1042
     >                (abs(f(i,j,k)-mdv).gt.eps)) then
1043
 
1044
                     vu = z3(i,j,k)-z3(i,j,o)
1045
                     vl = z3(i,j,u)-z3(i,j,k)
1046
                     vuvl = vu/(vl+zerodiv)
1047
                     vlvu = 1./(vuvl+zerodiv)
1048
 
1049
                     df(i,j,k) = 1./(vu+vl)
1050
     >                           * (vuvl*(f(i,j,u)-f(i,j,k)) 
1051
     >                           +  vlvu*(f(i,j,k)-f(i,j,o))) 
1052
 
1053
                  else
1054
                     df(i,j,k) = mdv
1055
                  endif
1056
 
1057
               enddo
1058
            enddo
1059
         enddo
1060
 
1061
c     Horizontal derivative in the y direction: 3d
1062
      elseif (direction.eq.'y') then
1063
 
1064
         do i=1,nx
1065
            do j=1,ny
1066
               do k=1,nz
1067
 
1068
                  s=j-1
1069
                  if (s.lt.1) s=1
1070
                  n=j+1
1071
                  if (n.gt.ny) n=ny
1072
 
1073
                  if ((abs(f(i,n,k)-mdv).gt.eps).and.
1074
     >                (abs(f(i,j,k)-mdv).gt.eps).and.
1075
     >                (abs(f(i,s,k)-mdv).gt.eps)) then  
1076
 
1077
                     vu = 1000.*(y2(i,j)-y2(i,n))
1078
                     vl = 1000.*(y2(i,s)-y2(i,j))
1079
                     vuvl = vu/(vl+zerodiv)
1080
                     vlvu = 1./(vuvl+zerodiv)
1081
 
1082
                     df(i,j,k) =  1./(vu+vl)
1083
     >                            * (vuvl*(f(i,s,k)-f(i,j,k))
1084
     >                            +  vlvu*(f(i,j,k)-f(i,n,k)))
1085
 
1086
                  else
1087
                     df(i,j,k) = mdv
1088
                  endif
1089
 
1090
               enddo
1091
            enddo
1092
         enddo
1093
 
1094
c     Horizontal derivative in the x direction: 3d
1095
      elseif (direction.eq.'x') then
1096
 
1097
         do i=1,nx
1098
            do j=1,ny
1099
               do k=1,nz
1100
 
1101
                  w=i-1
1102
                  if (w.lt.1) w=1
1103
                  e=i+1
1104
                  if (e.gt.nx) e=nx
1105
 
1106
                  if ((abs(f(w,j,k)-mdv).gt.eps).and.
1107
     >                (abs(f(i,j,k)-mdv).gt.eps).and.
1108
     >                (abs(f(e,j,k)-mdv).gt.eps)) then  
1109
 
1110
                     vu = 1000.*(x2(i,j)-x2(e,j))
1111
                     vl = 1000.*(x2(w,j)-x2(i,j))
1112
                     vuvl = vu/(vl+zerodiv)
1113
                     vlvu = 1./(vuvl+zerodiv)
1114
 
1115
                     df(i,j,k) =  1./(vu+vl)
1116
     >                            * (vuvl*(f(w,j,k)-f(i,j,k))
1117
     >                            +  vlvu*(f(i,j,k)-f(e,j,k)))
1118
 
1119
                  else
1120
                     df(i,j,k) = mdv
1121
                  endif
1122
 
1123
               enddo
1124
            enddo
1125
         enddo
1126
 
1127
c     Undefined direction for derivative
1128
      else
1129
 
1130
         print*,'Invalid direction of derivative... Stop'
1131
         stop
1132
 
1133
      endif
1134
 
1135
      end
1136
 
1137
c     -----------------------------------------------------------------
1138
c     Horizontal filter
1139
c     -----------------------------------------------------------------
1140
 
1141
      subroutine filt2d (a,af,nx,ny,fil,misdat,
1142
     &                   iperx,ipery,ispol,inpol)
1143
 
1144
c     Apply a conservative diffusion operator onto the 2d field a,
1145
c     with full missing data checking.
1146
c
1147
c     a     real   inp  array to be filtered, dimensioned (nx,ny)
1148
c     af    real   out  filtered array, dimensioned (nx,ny), can be
1149
c                       equivalenced with array a in the calling routine
1150
c     f1    real        workarray, dimensioned (nx+1,ny)
1151
c     f2    real        workarray, dimensioned (nx,ny+1)
1152
c     fil   real   inp  filter-coeff., 0<afil<=1. Maximum filtering with afil=1
1153
c                       corresponds to one application of the linear filter.
1154
c     misdat real  inp  missing-data value, a(i,j)=misdat indicates that
1155
c                       the corresponding value is not available. The
1156
c                       misdat-checking can be switched off with with misdat=0.
1157
c     iperx int    inp  periodic boundaries in the x-direction (1=yes,0=no)
1158
c     ipery int    inp  periodic boundaries in the y-direction (1=yes,0=no)
1159
c     inpol int    inp  northpole at j=ny  (1=yes,0=no)
1160
c     ispol int    inp  southpole at j=1   (1=yes,0=no)
1161
c
1162
c     Christoph Schaer, 1993
1163
 
1164
c     argument declaration
1165
      integer     nx,ny
1166
      real        a(nx,ny),af(nx,ny),fil,misdat
1167
      integer     iperx,ipery,inpol,ispol
1168
 
1169
c     local variable declaration
1170
      integer     i,j,is
1171
      real        fh
1172
      real        f1(nx+1,ny),f2(nx,ny+1)
1173
 
1174
c     compute constant fh
1175
      fh=0.125*fil
1176
 
1177
c     compute fluxes in x-direction
1178
      if (misdat.eq.0.) then
1179
        do j=1,ny
1180
        do i=2,nx
1181
          f1(i,j)=a(i-1,j)-a(i,j)
1182
        enddo
1183
        enddo
1184
      else
1185
        do j=1,ny
1186
        do i=2,nx
1187
          if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then
1188
            f1(i,j)=0.
1189
          else
1190
            f1(i,j)=a(i-1,j)-a(i,j)
1191
          endif
1192
        enddo
1193
        enddo
1194
      endif
1195
      if (iperx.eq.1) then
1196
c       do periodic boundaries in the x-direction
1197
        do j=1,ny
1198
          f1(1,j)=f1(nx,j)
1199
          f1(nx+1,j)=f1(2,j)
1200
        enddo
1201
      else
1202
c       set boundary-fluxes to zero
1203
        do j=1,ny
1204
          f1(1,j)=0.
1205
          f1(nx+1,j)=0.
1206
        enddo
1207
      endif
1208
 
1209
c     compute fluxes in y-direction
1210
      if (misdat.eq.0.) then
1211
        do j=2,ny
1212
        do i=1,nx
1213
          f2(i,j)=a(i,j-1)-a(i,j)
1214
        enddo
1215
        enddo
1216
      else
1217
        do j=2,ny
1218
        do i=1,nx
1219
          if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then
1220
            f2(i,j)=0.
1221
          else
1222
            f2(i,j)=a(i,j-1)-a(i,j)
1223
          endif
1224
        enddo
1225
        enddo
1226
      endif
1227
c     set boundary-fluxes to zero
1228
      do i=1,nx
1229
        f2(i,1)=0.
1230
        f2(i,ny+1)=0.
1231
      enddo
1232
      if (ipery.eq.1) then
1233
c       do periodic boundaries in the x-direction
1234
        do i=1,nx
1235
          f2(i,1)=f2(i,ny)
1236
          f2(i,ny+1)=f2(i,2)
1237
        enddo
1238
      endif
1239
      if (iperx.eq.1) then
1240
        if (ispol.eq.1) then
1241
c         do south-pole
1242
          is=(nx-1)/2
1243
          do i=1,nx
1244
            f2(i,1)=-f2(mod(i-1+is,nx)+1,2)
1245
          enddo
1246
        endif
1247
        if (inpol.eq.1) then
1248
c         do north-pole
1249
          is=(nx-1)/2
1250
          do i=1,nx
1251
            f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny)
1252
          enddo
1253
        endif
1254
      endif
1255
 
1256
c     compute flux-convergence -> filter
1257
      if (misdat.eq.0.) then
1258
        do j=1,ny
1259
        do i=1,nx
1260
            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
1261
        enddo
1262
        enddo
1263
      else
1264
        do j=1,ny
1265
        do i=1,nx
1266
          if (a(i,j).eq.misdat) then
1267
            af(i,j)=misdat
1268
          else
1269
            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
1270
          endif
1271
        enddo
1272
        enddo
1273
      endif
1274
      end
1275
 
1276
c     --------------------------------------------------------------------------------
1277
c     Read refernece profile from netcdf
1278
c     --------------------------------------------------------------------------------
1279
 
1280
      SUBROUTINE read_ref (nsqref,rhoref,thetaref,pressref,zref,
1281
     >                     nx,ny,nz,deltax,deltay,deltaz,coriol,oro,
1282
     >                     pvsrcfile)
1283
 
1284
c     Read the reference profile from file
1285
c
1286
c         thetaref             : Reference potential temperature (K)
1287
c         pressref             : Reference pressure (Pa)
1288
c         rhoref               : Reference density (kg/m^3)
1289
c         nsqref               : Stratification (s^-1)
1290
c         zref                 : Reference height (m)
1291
c         nx,nny,nz            : Grid dimension in x,y,z direction
1292
c         deltax,deltay,deltaz : Grid spacings used for calculations (m)
1293
c         coriol               : Coriolis parameter (s^-1)
1294
c         oro                  : Height of orography (m)
1295
c         pvsrcfile            : Input file
1296
 
1297
      implicit   none
1298
 
1299
c     Declaration of subroutine parameters
1300
      integer              nx,ny,nz
1301
      real                 nsqref  (0:2*nz)
1302
      real                 thetaref(0:2*nz)
1303
      real                 rhoref  (0:2*nz)
1304
      real                 pressref(0:2*nz)
1305
      real                 zref    (0:2*nz)
1306
      real                 deltax,deltay,deltaz
1307
      real                 coriol  (0:nx,0:ny)
1308
      real                 oro     (0:nx,0:ny)
1309
      character*80         pvsrcfile
1310
 
1311
c     Numerical and physical parameters
1312
      real                 eps
1313
      parameter            (eps=0.01)
1314
 
1315
c     Auxiliary variables
1316
      integer              cdfid,stat
1317
      integer              vardim(4)
1318
      real                 misdat
1319
      integer              ndimin
1320
      real                 varmin(4),varmax(4),stag(4)
1321
      integer              i,j,k,nf1
1322
      integer              ntimes
1323
      real                 time(200)
1324
      character*80         vnam(100),varname
1325
      integer              nvars
1326
      integer              isok,ierr
1327
      real                 x(0:nx,0:ny),y(0:nx,0:ny)
1328
      real                 mean,count
1329
 
1330
c     Get grid description from topography
1331
      call cdfopn(pvsrcfile,cdfid,stat)
1332
      if (stat.ne.0) goto 997
1333
      call getvars(cdfid,nvars,vnam,stat)
1334
      if (stat.ne.0) goto 997
1335
      isok=0
1336
      varname='ORO'
1337
      call check_varok(isok,varname,vnam,nvars)
1338
      if (isok.eq.0) goto 997
1339
      call getdef(cdfid,varname,ndimin,misdat,vardim,
1340
     >            varmin,varmax,stag,stat)
1341
      if (stat.ne.0) goto 997
1342
      time(1)=0.
1343
      call gettimes(cdfid,time,ntimes,stat)
1344
      if (stat.ne.0) goto 997
1345
      call clscdf(cdfid,stat)
1346
      if (stat.ne.0) goto 997
1347
 
1348
c     Open output netcdf file
1349
      call cdfopn(pvsrcfile,cdfid,stat)
1350
      if (stat.ne.0) goto 997
1351
 
1352
c     Create the variable if necessary
1353
      call getvars(cdfid,nvars,vnam,stat)
1354
      if (stat.ne.0) goto 997
1355
 
1356
c     Read data from netcdf file
1357
      isok=0
1358
      varname='NSQREF'
1359
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
1360
      call check_varok(isok,varname,vnam,nvars)
1361
      if (isok.eq.0) goto 997
1362
      call getdat(cdfid,varname,time(1),0,nsqref,stat)
1363
      if (stat.ne.0) goto 997
1364
 
1365
      isok=0
1366
      varname='RHOREF'
1367
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
1368
      call check_varok(isok,varname,vnam,nvars)
1369
      if (isok.eq.0) goto 997
1370
      call getdat(cdfid,varname,time(1),0,rhoref,stat)
1371
      if (stat.ne.0) goto 997
1372
 
1373
      isok=0
1374
      varname='THETAREF'
1375
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
1376
      call check_varok(isok,varname,vnam,nvars)
1377
      if (isok.eq.0) goto 997
1378
      call getdat(cdfid,varname,time(1),0,thetaref,stat)
1379
      if (stat.ne.0) goto 997
1380
 
1381
      isok=0
1382
      varname='PREREF'
1383
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
1384
      call check_varok(isok,varname,vnam,nvars)
1385
      if (isok.eq.0) goto 997
1386
      call getdat(cdfid,varname,time(1),0,pressref,stat)
1387
      if (stat.ne.0) goto 997
1388
 
1389
      isok=0
1390
      varname='ZREF'
1391
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
1392
      call check_varok(isok,varname,vnam,nvars)
1393
      if (isok.eq.0) goto 997
1394
      call getdat(cdfid,varname,time(1),0,zref,stat)
1395
      if (stat.ne.0) goto 997
1396
 
1397
      isok=0
1398
      varname='CORIOL'
1399
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
1400
      call check_varok(isok,varname,vnam,nvars)
1401
      if (isok.eq.0) goto 997
1402
      call getdat(cdfid,varname,time(1),0,coriol,stat)
1403
      if (stat.ne.0) goto 997
1404
 
1405
      isok=0
1406
      varname='ORO'
1407
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
1408
      call check_varok(isok,varname,vnam,nvars)
1409
      if (isok.eq.0) goto 997
1410
      call getdat(cdfid,varname,time(1),0,oro,stat)
1411
      if (stat.ne.0) goto 997
1412
 
1413
      isok=0
1414
      varname='X'
1415
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
1416
      call check_varok(isok,varname,vnam,nvars)
1417
      if (isok.eq.0) goto 997
1418
      call getdat(cdfid,varname,time(1),0,x,stat)
1419
      if (stat.ne.0) goto 997
1420
 
1421
      isok=0
1422
      varname='Y'
1423
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
1424
      call check_varok(isok,varname,vnam,nvars)
1425
      if (isok.eq.0) goto 997
1426
      call getdat(cdfid,varname,time(1),0,y,stat)
1427
      if (stat.ne.0) goto 997
1428
 
1429
c     Close  netcdf file
1430
      call clscdf(cdfid,stat)
1431
      if (stat.ne.0) goto 997
1432
 
1433
c     Determine the grid spacings <deltax, deltay, deltaz>
1434
      mean=0.
1435
      count=0.
1436
      do i=1,nx
1437
         do j=0,ny
1438
            mean=mean+abs(x(i)-x(i-1))
1439
            count=count+1.
1440
         enddo
1441
      enddo
1442
      deltax=mean/count
1443
 
1444
      mean=0.
1445
      count=0.
1446
      do j=1,ny
1447
         do i=0,nx
1448
            mean=mean+abs(y(j)-y(j-1))
1449
            count=count+1.
1450
         enddo
1451
      enddo
1452
      deltay=mean/count
1453
 
1454
      mean=0.
1455
      count=0.
1456
      do k=1,nz-1
1457
         mean=mean+abs(zref(k+1)-zref(k-1))
1458
         count=count+1.
1459
      enddo
1460
      deltaz=mean/count
1461
 
1462
      return
1463
 
1464
c     Exception handling
1465
 997  print*,'Read_Ref: Problem with input netcdf file... Stop'
1466
      stop
1467
 
1468
      end