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 hydrostatic
2
 
3
c     Calculate the geopotential and add it to the P file
4
c     Michael Sprenger / Spring 2006
5
 
6
      implicit none
7
 
8
c     ---------------------------------------------------------------
9
c     Declaration of variables
10
c     ---------------------------------------------------------------
11
 
12
c     Variables for input P file : model level
13
      character*80 ml_pfn
14
      real         ml_varmin(4),ml_varmax(4),ml_stag(4)
15
      integer      ml_vardim(4)
16
      real         ml_mdv
17
      integer      ml_ndim
18
      integer      ml_nx,ml_ny,ml_nz
19
      real         ml_xmin,ml_xmax,ml_ymin,ml_ymax,ml_dx,ml_dy
20
      integer      ml_ntimes
21
      real         ml_aklev(500),ml_bklev(500)
22
      real         ml_aklay(500),ml_bklay(500)
23
      real         ml_time
24
      real         ml_pollon,ml_pollat
25
      integer      ml_nvars
26
      character*80 ml_vnam(100)
27
      integer      ml_idate(5)
28
      real,allocatable, dimension (:,:)   :: ml_ps,ml_zb
29
      real,allocatable, dimension (:,:,:) :: ml_t3,ml_q3,ml_p3,ml_tv3
30
      real,allocatable, dimension (:,:,:) :: ml_zlay3
31
 
32
c     Variables for input Z file : pressure level
33
      character*80 pl_zfn
34
      real         pl_varmin(4),pl_varmax(4),pl_stag(4)
35
      integer      pl_vardim(4)
36
      real         pl_mdv
37
      integer      pl_ndim
38
      integer      pl_nx,pl_ny,pl_nz 
39
      real         pl_xmin,pl_xmax,pl_ymin,pl_ymax,pl_dx,pl_dy
40
      integer      pl_ntimes
41
      real         pl_aklev(500),pl_bklev(500)
42
      real         pl_aklay(500),pl_bklay(500)
43
      real         pl_time
44
      real         pl_pollon,pl_pollat
45
      integer      pl_nvars
46
      character*80 pl_vnam(100)
47
      integer      pl_idate(5)
48
      real,allocatable, dimension (:,:,:) :: pl_z3,pl_p3      
49
 
50
c     Physical and numerical parameters
51
      real      g
52
      parameter (g=9.80665)
53
      real      eps
54
      parameter (eps=0.01)
55
      real      tzero
56
      parameter (tzero=273.15)
57
      real      kappa
58
      parameter (kappa=0.6078)
59
      real      zerodiv
60
      parameter (zerodiv=0.0000000001)
61
      real      dpmin
62
      parameter (dpmin=10.)
63
      real      rdg
64
      parameter (rdg=29.271)
65
 
66
c     Auxiliary variables
67
      integer      ierr
68
      integer      cdfid,cstid
69
      character*80 cfn
70
      integer      stat
71
      real         time
72
      real         tv1(1000),z1(1000),p1(1000),f1(1000)
73
      real         spline_tv1(1000),spline_f1(1000),spline_z1(1000)
74
      real         pu,po,zu,zo,p,z,dp,p0,tvu,tvo,ff
75
      integer      i,j,k,l
76
      integer      lmin,n
77
      character*80 varname,cdfname
78
      integer      isok
79
      integer      mode
80
 
81
c     -----------------------------------------------------------------
82
c     Read input fields
83
c     -----------------------------------------------------------------
84
 
85
      print*,'*********************************************************'
86
      print*,'* hydrostatic                                           *'
87
      print*,'*********************************************************'
88
 
89
c     Read in the parameter file
90
      open(10,file='fort.10')
91
       read(10,*) ml_pfn
92
       read(10,*) pl_zfn
93
      close(10)
94
 
95
c     Decide which mode is used (1: reference from Z, 2: reference from ORO/PS)
96
      if ( trim(ml_pfn).ne.trim(pl_zfn) ) then
97
         mode=1
98
         print*,'Taking reference from Z ',trim(pl_zfn)
99
      else
100
         mode=2
101
         print*,'Taking reference from ORO/PS ',trim(pl_zfn)
102
      endif
103
 
104
      print*,trim(ml_pfn)
105
      print*,trim(pl_zfn)
106
 
107
c     Get grid description for P file : model level
108
      call cdfopn(ml_pfn,cdfid,ierr)
109
      if (ierr.ne.0) goto 998
110
      call getcfn(cdfid,cfn,ierr)
111
      if (ierr.ne.0) goto 998
112
      call cdfopn(cfn,cstid,ierr)   
113
      if (ierr.ne.0) goto 998
114
      call getvars(cdfid,ml_nvars,ml_vnam,ierr)
115
      varname='T'
116
      isok=0
117
      call check_varok (isok,varname,ml_vnam,ml_nvars)
118
      if (isok.ne.1) goto 998
119
      call getdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
120
     >            ml_varmin,ml_varmax,ml_stag,ierr)
121
      if (ierr.ne.0) goto 998
122
      ml_nx  =ml_vardim(1)
123
      ml_ny  =ml_vardim(2)
124
      ml_nz  =ml_vardim(3)
125
      ml_xmin=ml_varmin(1)
126
      ml_ymin=ml_varmin(2)
127
      call getlevs(cstid,ml_nz,ml_aklev,ml_bklev,ml_aklay,ml_bklay,ierr)
128
      call getgrid(cstid,ml_dx,ml_dy,ierr)
129
      ml_xmax=ml_xmin+real(ml_nx-1)*ml_dx
130
      ml_ymax=ml_ymin+real(ml_ny-1)*ml_dy
131
      call gettimes(cdfid,ml_time,ml_ntimes,ierr) 
132
      call getstart(cstid,ml_idate,ierr)
133
      call getpole(cstid,ml_pollon,ml_pollat,ierr)
134
      call clscdf(cstid,ierr)
135
      call clscdf(cdfid,ierr) 
136
 
137
c     Get grid description reference: either ORO(P) or Z(Z)
138
      call cdfopn(pl_zfn,cdfid,ierr)
139
      if (ierr.ne.0) goto 998
140
      call getcfn(cdfid,cfn,ierr)
141
      if (ierr.ne.0) goto 998
142
      call cdfopn(cfn,cstid,ierr)   
143
      if (ierr.ne.0) goto 998
144
      call getvars(cdfid,pl_nvars,pl_vnam,ierr)
145
      if (mode.eq.1) then
146
         varname='Z'
147
         isok=0
148
         call check_varok (isok,varname,pl_vnam,pl_nvars)
149
         if (isok.ne.1) goto 998
150
         call getdef(cdfid,varname,pl_ndim,pl_mdv,pl_vardim,
151
     >               pl_varmin,pl_varmax,pl_stag,ierr)
152
         if (ierr.ne.0) goto 998
153
         call getlevs(cstid,pl_nz,pl_aklev,pl_bklev,
154
     >                pl_aklay,pl_bklay,ierr)
155
         call getgrid(cstid,pl_dx,pl_dy,ierr)
156
         call gettimes(cdfid,pl_time,pl_ntimes,ierr)
157
         call getstart(cstid,pl_idate,ierr)
158
         call getpole(cstid,pl_pollon,pl_pollat,ierr)
159
 
160
      else if (mode.eq.2) then
161
         varname='ORO'
162
         isok=0
163
         call check_varok (isok,varname,pl_vnam,pl_nvars)
164
         if (isok.ne.1) goto 998
165
         call getdef(cdfid,varname,pl_ndim,pl_mdv,pl_vardim,
166
     >               pl_varmin,pl_varmax,pl_stag,ierr)
167
         if (ierr.ne.0) goto 998
168
         call getgrid(cstid,pl_dx,pl_dy,ierr)
169
         call gettimes(cdfid,pl_time,pl_ntimes,ierr)
170
         call getstart(cstid,pl_idate,ierr)
171
         call getpole(cstid,pl_pollon,pl_pollat,ierr)
172
 
173
      endif
174
      pl_nx  =pl_vardim(1)
175
      pl_ny  =pl_vardim(2)
176
      pl_nz  =pl_vardim(3)
177
      pl_xmin=pl_varmin(1)
178
      pl_ymin=pl_varmin(2)
179
      pl_xmax=pl_xmin+real(pl_nx-1)*pl_dx
180
      pl_ymax=pl_ymin+real(pl_ny-1)*pl_dy
181
      call clscdf(cstid,ierr)
182
      call clscdf(cdfid,ierr) 
183
 
184
c     Consitency check for the grids
185
      if ( (ml_nx.ne.pl_nx).or.
186
     >     (ml_ny.ne.pl_ny).or.
187
     >     (abs(ml_xmin-pl_xmin    ).gt.eps).or.
188
     >     (abs(ml_ymin-pl_ymin    ).gt.eps).or.
189
     >     (abs(ml_xmax-pl_xmax    ).gt.eps).or.
190
     >     (abs(ml_ymax-pl_ymax    ).gt.eps).or.
191
     >     (abs(ml_dx  -pl_dx      ).gt.eps).or.
192
     >     (abs(ml_dy  -pl_dy      ).gt.eps).or.
193
     >     (abs(ml_time-pl_time    ).gt.eps).or.
194
     >     (abs(ml_pollon-pl_pollon).gt.eps).or.
195
     >     (abs(ml_pollat-pl_pollat).gt.eps)) then
196
         print*,'Input P and Z grids are not consistent... Stop'
197
         print*,'Xmin   ',ml_xmin    ,pl_xmin
198
         print*,'Ymin   ',ml_ymin    ,pl_ymin
199
         print*,'Xmax   ',ml_xmax    ,pl_xmax
200
         print*,'Ymax   ',ml_ymax    ,pl_ymax
201
         print*,'Dx     ',ml_dx      ,pl_dx
202
         print*,'Dy     ',ml_dy      ,pl_dy
203
         print*,'Pollon ',ml_pollon  ,pl_pollon
204
         print*,'Pollat ',ml_pollat  ,pl_pollat
205
         print*,'Time   ',ml_time    ,pl_time
206
         stop
207
      endif
208
 
209
c     Allocate memory for all fields
210
      allocate(ml_ps(ml_nx,ml_ny),stat=stat)
211
      if (stat.ne.0) print*,'*** error allocating array ml_ps ***'
212
      allocate(ml_zb(ml_nx,ml_ny),stat=stat)
213
      if (stat.ne.0) print*,'*** error allocating array ml_zb ***'
214
      allocate(ml_p3(ml_nx,ml_ny,ml_nz),stat=stat)
215
      if (stat.ne.0) print*,'*** error allocating array ml_p3 ***'
216
      allocate(ml_t3(ml_nx,ml_ny,ml_nz),stat=stat)
217
      if (stat.ne.0) print*,'*** error allocating array ml_t3 ***'
218
      allocate(ml_q3(ml_nx,ml_ny,ml_nz),stat=stat)
219
      if (stat.ne.0) print*,'*** error allocating array ml_q3 ***'
220
      allocate(ml_tv3(ml_nx,ml_ny,ml_nz),stat=stat)
221
      if (stat.ne.0) print*,'*** error allocating array ml_tv3 ***'
222
      allocate(ml_zlay3(ml_nx,ml_ny,ml_nz),stat=stat)
223
      if (stat.ne.0) print*,'*** error allocating array ml_zlay3 ***'
224
 
225
      allocate(pl_z3(pl_nx,pl_ny,pl_nz),stat=stat)
226
      if (stat.ne.0) print*,'*** error allocating array pl_z3 ***'
227
      allocate(pl_p3(pl_nx,pl_ny,pl_nz),stat=stat)
228
      if (stat.ne.0) print*,'*** error allocating array pl_p3 ***'
229
 
230
c     Read T, Q, PS from P file
231
      call cdfopn(ml_pfn,cdfid,ierr)
232
      if (ierr.ne.0) goto 998
233
      isok=0
234
      varname='T'
235
      call check_varok (isok,varname, ml_vnam,ml_nvars)
236
      varname='Q'
237
      call check_varok (isok,varname, ml_vnam,ml_nvars)
238
      varname='PS'
239
      call check_varok (isok,varname,ml_vnam,ml_nvars)
240
      if (isok.ne.3) goto 998
241
      print*,'R T  ',trim(ml_pfn)
242
      call getdat(cdfid,'T',ml_time,0,ml_t3,ierr)
243
      print*,'R Q  ',trim(ml_pfn)
244
      call getdat(cdfid,'Q',ml_time,0,ml_q3,ierr)
245
      print*,'R PS ',trim(ml_pfn)
246
      call getdat(cdfid,'PS',ml_time,1,ml_ps,ierr)
247
      call clscdf(cdfid,ierr)
248
 
249
c     Read ORO from P or Z from Z file
250
      call cdfopn(pl_zfn,cdfid,ierr)
251
      if (ierr.ne.0) goto 998
252
      if (mode.eq.1) then
253
         isok=0
254
         varname='Z'
255
         call check_varok (isok,varname,pl_vnam,pl_nvars)
256
         if (isok.ne.1) goto 998
257
         print*,'R Z  ',trim(pl_zfn)
258
         call getdat(cdfid,varname,pl_time,0,pl_z3,ierr)
259
      else if (mode.eq.2) then
260
         isok=0
261
         varname='ORO'
262
         call check_varok (isok,varname,pl_vnam,pl_nvars)
263
         if (isok.ne.1) goto 998
264
         print*,'R ORO  ',trim(pl_zfn)
265
         call getdat(cdfid,varname,pl_time,1,pl_z3,ierr)
266
      endif
267
      call clscdf(cdfid,ierr)
268
 
269
c     Set the values for the pressure on the pressure levels
270
      do i=1,pl_nx
271
         do j=1,pl_ny
272
            do k=1,pl_nz
273
               if (mode.eq.1) then
274
                  pl_p3(i,j,k)=pl_aklay(k)
275
               else if (mode.eq.2) then
276
                  pl_p3(i,j,k)=ml_ps(i,j)
277
               endif
278
            enddo
279
         enddo
280
      enddo
281
 
282
c     Calculate 3d pressure field
283
      print*,'C P'
284
      do k=1,ml_nz
285
        do i=1,ml_nx
286
           do j=1,ml_ny
287
              ml_p3(i,j,k)=ml_aklay(k)+ml_bklay(k)*ml_ps(i,j)
288
            enddo
289
         enddo
290
      enddo
291
 
292
c     Calculate 3d virtual temperature
293
      print*,'C TV'
294
      do k=1,ml_nz
295
         do i=1,ml_nx
296
            do j=1,ml_ny
297
               ml_tv3(i,j,k) = (ml_t3(i,j,k)+tzero)*
298
     >                         (1.+kappa*ml_q3(i,j,k))
299
            enddo
300
         enddo
301
      enddo
302
 
303
c     -----------------------------------------------------------------
304
c     Calculate geopotential 
305
c     -----------------------------------------------------------------
306
 
307
c     Integrate hydrostatic equation towards layers
308
      print*,'C HYDROSTATIC EQUATION (LAYERS)'
309
      do i=1,ml_nx
310
         do j=1,ml_ny
311
 
312
c           Make the virtual temperature profile available 
313
            do k=1,ml_nz
314
               p1 (ml_nz-k+1)=ml_p3 (i,j,k)
315
               tv1(ml_nz-k+1)=ml_tv3(i,j,k)
316
            enddo
317
            call spline (p1,tv1,ml_nz,1.e30,1.e30,spline_tv1)
318
 
319
c           Loop over all model levels
320
            do k=1,ml_nz
321
 
322
c              Get pressure at the grid point
323
               p = ml_aklay(k)+ml_bklay(k)*ml_ps(i,j)
324
 
325
c              Find nearest pressure level which is above topography
326
               if (mode.eq.1) then
327
                  lmin=pl_nz
328
                  do l=1,pl_nz
329
                     if ((abs(p-pl_p3(i,j,l))).lt.
330
     >                   (abs(p-pl_p3(i,j,lmin)))
331
     >                    .and.
332
     >                    (pl_p3(i,j,l).lt.ml_ps(i,j)) ) then
333
                        lmin=l
334
                     endif
335
                  enddo
336
               else if (mode.eq.2) then
337
                  lmin=1
338
               endif
339
 
340
c              Integrate hydrostatic equation from this level to the grid point
341
               p0 = pl_p3(i,j,lmin)
342
               n  = nint(abs(p-p0)/dpmin)
343
               if (n.lt.1) n=1
344
               dp = (p-p0)/real(n)
345
 
346
               pu  = p0
347
               z   = pl_z3(i,j,lmin)
348
               call splint(p1,tv1,spline_tv1,ml_nz,pu,tvu)
349
               do l=1,n
350
                  po  = pu+dp
351
                  call splint(p1,tv1,spline_tv1,ml_nz,po,tvo)
352
                  z   = z + rdg*0.5*(tvu+tvo)*alog(pu/po)
353
                  tvu = tvo
354
                  pu  = po
355
               enddo
356
 
357
c              Set the geopotential at the grid point
358
               ml_zlay3(i,j,k) = z
359
 
360
            enddo
361
 
362
         enddo
363
 
364
      enddo
365
 
366
c     -----------------------------------------------------------------
367
c     Calculate height of topography 
368
c     -----------------------------------------------------------------
369
 
370
      if (mode.eq.1) then 
371
         print*,'C TOPOGRAPHY'
372
         do i=1,ml_nx
373
            do j=1,ml_ny
374
 
375
c              Make the z/p profile available 
376
               do k=1,ml_nz
377
                  p1(ml_nz-k+1)=ml_p3(i,j,k)
378
                  z1(ml_nz-k+1)=ml_zlay3(i,j,k)
379
               enddo
380
 
381
c              Cubic spline interpolation
382
               call spline (p1,z1,ml_nz,1.e30,1.e30,spline_z1) 
383
               call splint (p1,z1,spline_z1,ml_nz,ml_ps(i,j),ml_zb(i,j))
384
 
385
            enddo
386
         enddo
387
      endif
388
 
389
c     -----------------------------------------------------------------
390
c     Write the topography and geopotential to P file
391
c     -----------------------------------------------------------------
392
 
393
c     Open P file
394
      call cdfwopn(trim(ml_pfn),cdfid,ierr)   
395
 
396
c     Write orography (levels)
397
      if (mode.eq.1) then
398
         varname='ORO'
399
         print*,'W ',trim(varname),' ',trim(ml_pfn)
400
         isok=0
401
         call check_varok (isok,varname,ml_vnam,ml_nvars)
402
         if (isok.eq.0) then
403
            ml_vardim(3)=1
404
            call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
405
     >                  ml_varmin,ml_varmax,ml_stag,ierr)         
406
            ml_vardim(3)=ml_nz
407
            if (ierr.ne.0) goto 997
408
         endif
409
         call putdat(cdfid,varname,ml_time,1,ml_zb,ierr)
410
         if (ierr.ne.0) goto 997
411
      endif
412
 
413
c     Write geopotential on layers
414
      varname='Z_DIA'
415
      print*,'W ',trim(varname),' ',trim(ml_pfn)
416
      isok=0
417
      call check_varok (isok,varname,ml_vnam,ml_nvars)
418
      if (isok.eq.0) then
419
         ml_stag(3)=-0.5
420
         call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
421
     >               ml_varmin,ml_varmax,ml_stag,ierr)         
422
         if (ierr.ne.0) goto 997
423
      endif
424
      call putdat(cdfid,varname,ml_time,0,ml_zlay3,ierr)
425
      if (ierr.ne.0) goto 997
426
 
427
c     Close P file
428
      call clscdf(cdfid,ierr)
429
 
430
c     -----------------------------------------------------------------
431
c     Exception handling and format specs
432
c     -----------------------------------------------------------------
433
 
434
      stop
435
 
436
 998  print*,'Z: Problems with input from m level'
437
      stop
438
 
439
 997  print*,'Z: Problems with output on m level'
440
      stop
441
 
442
 996  print*,'F: Problems with input from m level'
443
      stop
444
 
445
 995  print*,'F: Problems with output on z level'
446
      stop
447
 
448
 
449
      end
450
 
451
c     -------------------------------------------------------------
452
c     Natural cubic spline
453
c     -------------------------------------------------------------
454
 
455
      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
456
      INTEGER n,NMAX
457
      REAL yp1,ypn,x(n),y(n),y2(n)
458
      PARAMETER (NMAX=500)
459
      INTEGER i,k
460
      REAL p,qn,sig,un,u(NMAX)
461
      if (yp1.gt..99e30) then
462
        y2(1)=0.
463
        u(1)=0.
464
      else
465
        y2(1)=-0.5
466
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
467
      endif
468
      do 11 i=2,n-1
469
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
470
        p=sig*y2(i-1)+2.
471
        y2(i)=(sig-1.)/p
472
        u(i)=(6.*((y(i+1)-y(i))/(x(i+
473
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
474
     *u(i-1))/p
475
11    continue
476
      if (ypn.gt..99e30) then
477
        qn=0.
478
        un=0.
479
      else
480
        qn=0.5
481
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
482
      endif
483
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
484
      do 12 k=n-1,1,-1
485
        y2(k)=y2(k)*y2(k+1)+u(k)
486
12    continue
487
      return
488
      END
489
 
490
 
491
      SUBROUTINE splint(xa,ya,y2a,n,x,y)
492
      INTEGER n
493
      REAL x,y,xa(n),y2a(n),ya(n)
494
      INTEGER k,khi,klo
495
      REAL a,b,h
496
      klo=1
497
      khi=n
498
1     if (khi-klo.gt.1) then
499
        k=(khi+klo)/2
500
        if(xa(k).gt.x)then
501
          khi=k
502
        else
503
          klo=k
504
        endif
505
      goto 1
506
      endif
507
      h=xa(khi)-xa(klo)
508
      if (h.eq.0.) pause 'bad xa input in splint'
509
      a=(xa(khi)-x)/h
510
      b=(x-xa(klo))/h
511
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
512
     *2)/6.
513
      return
514
      END
515
 
516
c     ----------------------------------------------------------------
517
c     Check whether variable is found on netcdf file
518
c     ----------------------------------------------------------------
519
 
520
      subroutine check_varok (isok,varname,varlist,nvars)
521
 
522
c     Check whether the variable <varname> is in the list <varlist(nvars)>. 
523
c     If this is the case, <isok> is incremented by 1. Otherwise <isok> 
524
c     keeps its value.
525
 
526
      implicit none
527
 
528
c     Declaraion of subroutine parameters
529
      integer      isok
530
      integer      nvars
531
      character*80 varname
532
      character*80 varlist(nvars)
533
 
534
c     Auxiliary variables
535
      integer      i
536
 
537
c     Main
538
      do i=1,nvars
539
         if (trim(varname).eq.trim(varlist(i))) isok=isok+1
540
      enddo
541
 
542
      end