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 inv_cart
2
 
3
c     ********************************************************************************
4
c     * CALCULATES FOR A PV AND THETA DISTRIBUTION OTHER PROGNOSTIC FIELDS BY        *
5
c     * MEANS OF A PV INVERSION                                                      *
6
c     *                                                                              *
7
c     * Rene Fehlmann 1994 / Code re-organization: Michael Sprenger, 2006            *
8
c     ********************************************************************************
9
 
10
c     --------------------------------------------------------------------------------
11
c     Declaration of variables, parameters, externals and common blocks            
12
c     --------------------------------------------------------------------------------
13
 
14
      implicit none
15
 
16
c     Grid parameters
17
      integer        nx,ny,nz
18
      real           xmin,ymin,zmin,xmax,ymax,zmax
19
      real           dx,dy,dz
20
      real           deltax,deltay,deltaz
21
      real           mdv
22
      real,          allocatable,dimension (:,:) :: coriol
23
 
24
c     Reference state
25
      real,   allocatable,dimension (:) :: nsqref
26
      real,   allocatable,dimension (:) :: thetaref
27
      real,   allocatable,dimension (:) :: rhoref
28
      real,   allocatable,dimension (:) :: pressref
29
      real,   allocatable,dimension (:) :: zref
30
 
31
 
32
C     Boundary conditions
33
      real,   allocatable,dimension (:,:) :: thetatop
34
      real,   allocatable,dimension (:,:) :: thetabot
35
      real,   allocatable,dimension (:,:) :: ua
36
      real,   allocatable,dimension (:,:) :: ub
37
      real,   allocatable,dimension (:,:) :: va
38
      real,   allocatable,dimension (:,:) :: vb
39
 
40
c     Potentiual vorticity and stream function
41
      real,   allocatable,dimension (:,:,:) :: psi
42
      real,   allocatable,dimension (:,:,:) :: pv
43
 
44
c     Auxiliary arrays for numerics
45
      real,   allocatable,dimension (:,:,:) :: a
46
      real,   allocatable,dimension (:,:,:) :: b
47
      real,   allocatable,dimension (:,:,:) :: c
48
 
49
c     Input parameters
50
      character*80         pvsrcfile
51
      character*80         referfile
52
 
53
c     Auxiliary variables
54
      integer             i,j,k
55
      integer             stat
56
 
57
c     --------------------------------------------------------------------------------
58
c     Preparations
59
c     --------------------------------------------------------------------------------
60
 
61
      print*,'********************************************************'
62
      print*,'* INV_CART                                             *'
63
      print*,'********************************************************'
64
 
65
c     Read parameter file
66
      open(10,file='fort.10')
67
        read(10,*) pvsrcfile
68
        read(10,*) referfile
69
      close(10)
70
      print*
71
      print*,trim(pvsrcfile)
72
 
73
 
74
c     Get lat/lon gid parameters from input file
75
      call read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv,
76
     >               pvsrcfile)
77
      print*
78
      print*,'Read_Dim: nx,ny,nz         = ',nx,ny,nz
79
      print*,'          dx,dy,dz         = ',dx,dy,dz
80
      print*,'          xmin,ymin,zmin   = ',xmin,ymin,zmin
81
      print*,'          mdv              = ',mdv
82
      print*
83
 
84
c     Count from 0, not from 1
85
      nx=nx-1
86
      ny=ny-1
87
      nz=nz-1
88
 
89
c     Allocate memory for boundary conditions
90
      allocate(thetatop(0:nx,0:ny),stat=stat)
91
      if (stat.ne.0) print*,'*** error allocating array thetatop ***'
92
      allocate(thetabot(0:nx,0:ny),stat=stat)
93
      if (stat.ne.0) print*,'*** error allocating array thetabot ***'
94
      allocate(ua(0:nx,0:nz),stat=stat)
95
      if (stat.ne.0) print*,'*** error allocating array ua ***'
96
      allocate(ub(0:nx,0:nz),stat=stat)
97
      if (stat.ne.0) print*,'*** error allocating array ub ***'
98
      allocate(va(0:ny,0:nz),stat=stat)
99
      if (stat.ne.0) print*,'*** error allocating array va ***'
100
      allocate(vb(0:ny,0:nz),stat=stat)
101
      if (stat.ne.0) print*,'*** error allocating array vb ***'
102
 
103
c     Allocate memory for 3d PV and stream function
104
      allocate(psi(0:nx,0:ny,0:nz),stat=stat)
105
      if (stat.ne.0) print*,'*** error allocating array psi ***'
106
      allocate(pv(0:nx,0:ny,0:nz),stat=stat)
107
      if (stat.ne.0) print*,'*** error allocating array pv ***'
108
 
109
c     Alllocate memory for matrix elements for inversion operator
110
      allocate(a(0:nx,0:ny,0:nz),stat=stat)
111
      if (stat.ne.0) print*,'*** error allocating array a ***'
112
      allocate(b(0:nx,0:ny,0:nz),stat=stat)
113
      if (stat.ne.0) print*,'*** error allocating array b ***'
114
      allocate(c(0:nx,0:ny,0:nz),stat=stat)
115
      if (stat.ne.0) print*,'*** error allocating array c ***'
116
 
117
c     Allocate memory for reference profile
118
      allocate(nsqref(0:2*nz),stat=stat)
119
      if (stat.ne.0) print*,'*** error allocating array nsqref ***'
120
      allocate(thetaref(0:2*nz),stat=stat)
121
      if (stat.ne.0) print*,'*** error allocating array thetaref ***'
122
      allocate(rhoref(0:2*nz),stat=stat)
123
      if (stat.ne.0) print*,'*** error allocating array rhoref ***'
124
      allocate(pressref(0:2*nz),stat=stat)
125
      if (stat.ne.0) print*,'*** error allocating array pressref ***'
126
      allocate(zref(0:2*nz),stat=stat)
127
      if (stat.ne.0) print*,'*** error allocating array zref ***'
128
 
129
c     Allocate memory for Coriolis parameter
130
      allocate(coriol(0:nx,0:ny),STAT=stat)
131
      if (stat.ne.0) print*,'error allocating coriol'
132
 
133
c     --------------------------------------------------------------------------------
134
c     Input
135
c     --------------------------------------------------------------------------------
136
 
137
c     Read reference profile and grid parameters
138
      call read_ref (nsqref,rhoref,thetaref,pressref,zref,
139
     >               nx,ny,nz,deltax,deltay,deltaz,coriol,
140
     >               referfile)
141
      deltax=1000.*deltax
142
      deltay=1000.*deltay
143
 
144
c     Read input fields from netcdf 
145
      call read_inp (pv,thetabot,thetatop,ua,ub,va,vb,
146
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,
147
     >               pvsrcfile)
148
 
149
c     --------------------------------------------------------------------------------
150
c     Perform the inversion
151
c     --------------------------------------------------------------------------------
152
 
153
C     Init matrix elements for inversion
154
      call matrixel(a,b,c,coriol,nx,ny,nz,nsqref,rhoref,
155
     >              deltax,deltay,deltaz)
156
 
157
c     Inversion
158
      call sor(psi,nsqref,thetatop,thetabot,thetaref,rhoref,coriol,
159
     >         pv,ua,ub,va,vb,a,b,c,nx,ny,nz,deltax,deltay,deltaz)
160
 
161
c     --------------------------------------------------------------------------------
162
c     Output
163
c     --------------------------------------------------------------------------------
164
 
165
c     Write output to netcdf
166
      print*
167
      call write_out (psi,thetabot,thetatop,ua,ub,va,vb,
168
     >                nx,ny,nz,deltax,deltay,deltaz,
169
     >                coriol,thetaref,rhoref,pressref,pvsrcfile)
170
 
171
 
172
      END
173
 
174
 
175
c     ********************************************************************************
176
c     * NETCDF AND ASCII INPUT AND OUTPUT                                            *
177
c     ********************************************************************************
178
 
179
c     --------------------------------------------------------------------------------
180
c     Output to netcdf file
181
c     --------------------------------------------------------------------------------
182
 
183
      SUBROUTINE write_out (psi,thetabot,thetatop,ua,ub,va,vb,nx,ny,nz,
184
     >                      dx,dy,dz,coriol,thetaref,rhoref,pref,
185
     >                      pvsrcfile)
186
 
187
c     Write the result of the inversion to output netcdf
188
c     
189
c        psi               : streamm function as calculated from inversion
190
c        thetabot,thetatop : potential temperature at lower and upper boundary
191
c        ua,ub             : Zonal wind at western and eastern boundary
192
c        va,vb             : Meridional wind at southern and northern boundary
193
c        nx,ny,nz          : Grid dimension in x, y and z direction
194
c        dx,dy,dz          : Grid resolution
195
c        coriol            : Coriolis parameter
196
c        thetaref          : Reference profile of potential temperature
197
c        rhoref            : Reference profile of density
198
c        pref              : Reference profile of pressure
199
c        pvsrcfile         : Name of the output netcdf file
200
 
201
      implicit   none
202
 
203
c     Declaration of subroutine parameters
204
      integer              nx,ny,nz
205
      real                 psi(0:nx,0:ny,0:nz)
206
      real                 thetatop(0:nx,0:ny)
207
      real                 thetabot(0:nx,0:ny)
208
      real                 ua(0:nx,0:nz)
209
      real                 ub(0:nx,0:nz)
210
      real                 va(0:ny,0:nz)
211
      real                 vb(0:ny,0:nz)
212
      character*(30)       pvsrcfile
213
      real                 dx,dy,dz
214
      real                 thetaref(0:2*nz)
215
      real                 rhoref(0:2*nz)
216
      real                 pref(0:2*nz)
217
      real                 coriol(0:nx,0:ny)
218
 
219
c     Numerical and physical parameters
220
      real                 eps
221
      parameter            (eps=0.01)
222
      real                 g
223
      parameter            (g=9.81)
224
      real                 preunit
225
      parameter            (preunit=0.01)
226
      real                 kappa
227
      parameter            (kappa=0.286)
228
 
229
c     Auxiliary variables
230
      integer              cdfid,stat
231
      integer              vardim(4)
232
      real                 misdat
233
      integer              ndimin
234
      real                 varmin(4),varmax(4),stag(4)
235
      integer              i,j,k,nf1,jj,kk
236
      real                 tmp(0:nx,0:ny,0:nz)
237
      real                 pr (0:nx,0:ny,0:nz)
238
      real                 th (0:nx,0:ny,0:nz)
239
      integer              ntimes
240
      real                 time(200)  
241
      character*80         vnam(100),varname
242
      integer              nvars
243
      integer              isok,ierr
244
      real                 meanpsi,meancnt
245
 
246
c     Get grid description 
247
      call cdfopn(pvsrcfile,cdfid,stat)
248
      if (stat.ne.0) goto 997
249
      call getvars(cdfid,nvars,vnam,stat)
250
      if (stat.ne.0) goto 997
251
      isok=0
252
      varname='QGPV'
253
      call check_varok(isok,varname,vnam,nvars)
254
      if (isok.eq.0) goto 997
255
      call getdef(cdfid,varname,ndimin,misdat,vardim,
256
     >            varmin,varmax,stag,stat)    
257
      if (stat.ne.0) goto 997
258
      time(1)=0.
259
      call gettimes(cdfid,time,ntimes,stat)  
260
      if (stat.ne.0) goto 997
261
      call clscdf(cdfid,stat)
262
      if (stat.ne.0) goto 997
263
 
264
c     Open output netcdf file
265
      call cdfwopn(pvsrcfile,cdfid,stat)
266
      if (stat.ne.0) goto 997
267
 
268
c     Create the variable if necessary
269
      call getvars(cdfid,nvars,vnam,stat)
270
      if (stat.ne.0) goto 997
271
      isok=0
272
      varname='PSI'
273
      call check_varok(isok,varname,vnam,nvars)
274
      if (isok.eq.0) then
275
         call putdef(cdfid,varname,ndimin,misdat,vardim,
276
     >               varmin,varmax,stag,stat)
277
         if (stat.ne.0) goto 997
278
      endif
279
      isok=0
280
      varname='U'
281
      call check_varok(isok,varname,vnam,nvars)
282
      if (isok.eq.0) then
283
         call putdef(cdfid,varname,ndimin,misdat,vardim,
284
     >               varmin,varmax,stag,stat)
285
         if (stat.ne.0) goto 997
286
      endif
287
      isok=0
288
      varname='V'
289
      call check_varok(isok,varname,vnam,nvars)
290
      if (isok.eq.0) then
291
         call putdef(cdfid,varname,ndimin,misdat,vardim,
292
     >               varmin,varmax,stag,stat)
293
         if (stat.ne.0) goto 997
294
      endif
295
      isok=0
296
      varname='TH'
297
      call check_varok(isok,varname,vnam,nvars)
298
      if (isok.eq.0) then
299
         call putdef(cdfid,varname,ndimin,misdat,vardim,
300
     >               varmin,varmax,stag,stat)
301
         if (stat.ne.0) goto 997
302
      endif
303
      isok=0
304
      varname='T'
305
      call check_varok(isok,varname,vnam,nvars)
306
      if (isok.eq.0) then
307
         call putdef(cdfid,varname,ndimin,misdat,vardim,
308
     >               varmin,varmax,stag,stat)
309
         if (stat.ne.0) goto 997
310
      endif
311
      isok=0
312
      varname='P'
313
      call check_varok(isok,varname,vnam,nvars)
314
      if (isok.eq.0) then
315
         call putdef(cdfid,varname,ndimin,misdat,vardim,
316
     >               varmin,varmax,stag,stat)
317
         if (stat.ne.0) goto 997
318
      endif
319
 
320
c     Write stream function
321
      varname='PSI'
322
      call putdat(cdfid,varname,time(1),0,psi,ierr)
323
      if (stat.ne.0) goto 997
324
      print*,'W PSI      ',trim(pvsrcfile)
325
 
326
c     Calculate and write velocity U: keep southern and northern boundary
327
      do k=0,nz
328
         do i=0,nx
329
            do j=1,ny-1
330
               tmp(i,j,k)=(psi(i,j-1,k)-psi(i,j+1,k))/
331
     >                    (2.*dy*coriol(i,j))
332
            enddo
333
            tmp(i,0,k) =ua(i,k)
334
            tmp(i,ny,k)=ub(i,k)
335
         enddo
336
      enddo
337
      varname='U'
338
      call putdat(cdfid,varname,time(1),0,tmp,ierr)
339
      if (stat.ne.0) goto 997
340
      print*,'W U        ',trim(pvsrcfile)
341
 
342
C     Calculate and write velocity V: keep western and eastern boundary
343
      do k=0,nz
344
         do j=0,ny
345
            do i=1,nx-1
346
               tmp(i,j,k)=(psi(i+1,j,k)-psi(i-1,j,k))/
347
     >                    (2.*dx*coriol(i,j))
348
            enddo
349
            tmp(0,j,k)=va(j,k)
350
            tmp(nx,j,k)=vb(j,k)
351
         enddo
352
      enddo
353
      varname='V'
354
      call putdat(cdfid,varname,time(1),0,tmp,ierr)
355
      if (stat.ne.0) goto 997
356
      print*,'W V        ',trim(pvsrcfile)
357
 
358
c     Calculate and write potential temperature: keep lower and upper boundary
359
c     Potential temperature is needed for calculation of temperature: keep it
360
      do i=0,nx
361
         do j=0,ny
362
            th(i,j,0)=thetabot(i,j)
363
            do k=1,nz-1      
364
              th(i,j,k)=thetaref(2*k)*
365
     >                  (psi(i,j,k+1)-psi(i,j,k-1))/(2.*dz*g)
366
 
367
           enddo
368
           th(i,j,nz)=thetatop(i,j)
369
         enddo
370
      enddo    
371
      varname='TH'
372
      call putdat(cdfid,varname,time(1),0,th,ierr)
373
      if (stat.ne.0) goto 997
374
      print*,'W TH       ',trim(pvsrcfile)
375
 
376
c     Calculate and write pressure: The pressure is directly proportional to the
377
c     streamfunction. But the streamfunction is determined only up to an additive
378
c     constant. Shift the streamfunction in such a way that it vanish in the mean
379
c     on the boundaries. Pressure is needed for calculation of temperature: keep it
380
      meanpsi=0.
381
      meancnt=0.
382
      do i=0,nx
383
         do j=0,ny
384
            meanpsi=meanpsi+psi(i,j,0)+psi(i,j,nz)
385
            meancnt=meancnt+2
386
         enddo
387
      enddo
388
      do i=0,nx
389
         do k=0,nz
390
            meanpsi=meanpsi+psi(i,0,k)+psi(i,ny,k)
391
            meancnt=meancnt+2 
392
         enddo
393
      enddo
394
      do j=0,ny
395
         do k=0,nz
396
            meanpsi=meanpsi+psi(0,j,k)+psi(nx,j,k)
397
            meancnt=meancnt+2 
398
         enddo
399
      enddo
400
      meanpsi=meanpsi/meancnt
401
      do i=0,nx
402
         do j=0,ny
403
            do k=0,nz
404
               kk=2*k
405
               pr(i,j,k)=preunit*rhoref(kk)*(psi(i,j,k)-meanpsi)
406
            enddo
407
         enddo
408
      enddo
409
      varname='P'
410
      call putdat(cdfid,varname,time(1),0,pr,ierr)
411
      if (stat.ne.0) goto 997
412
      print*,'W P        ',trim(pvsrcfile)
413
 
414
c     Calculate and write temperature
415
      do i=0,nx
416
         do j=0,ny
417
            do k=0,nz      
418
               kk=2*k
419
               tmp(i,j,k)=((pref(kk)/1000.)**kappa) * 
420
     >                (th(i,j,k)+kappa*thetaref(kk)*pr(i,j,k)/pref(kk))
421
           enddo
422
        enddo
423
      enddo      
424
      varname='T'
425
      call putdat(cdfid,varname,time(1),0,tmp,ierr)
426
      if (stat.ne.0) goto 997
427
      print*,'W T        ',trim(pvsrcfile)
428
 
429
c     Close output netcdf file
430
      call clscdf(cdfid,stat)
431
      if (stat.ne.0) goto 997
432
 
433
      return
434
 
435
c     Exception handling
436
 997  print*,'Problem with output netcdf file... Stop'
437
      stop
438
 
439
      end
440
 
441
 
442
c     --------------------------------------------------------------------------------
443
c     Read the reference file
444
c     --------------------------------------------------------------------------------
445
 
446
      SUBROUTINE read_ref (nsqref,rhoref,thetaref,pressref,zref,
447
     >                     nx,ny,nz,deltax,deltay,deltaz,coriol,
448
     >                     pvsrcfile) 
449
 
450
c     Read the reference profile from file
451
c
452
c         thetaref             : Reference potential temperature (K)
453
c         pressref             : Reference pressure (Pa)
454
c         rhoref               : Reference density (kg/m^3)
455
c         nsqref               : Stratification (s^-1)
456
c         zref                 : Reference height (m)
457
c         nx,nny,nz            : Grid dimension in x,y,z direction
458
c         deltax,deltay,deltaz : Grid spacings used for calculations (m)
459
c         coriol               : Coriolis parameter (s^-1)
460
c         pvsrcfile            : Input file
461
 
462
      implicit   none
463
 
464
c     Declaration of subroutine parameters
465
      integer              nx,ny,nz
466
      real                 nsqref  (0:2*nz)
467
      real                 thetaref(0:2*nz)
468
      real                 rhoref  (0:2*nz)
469
      real                 pressref(0:2*nz)
470
      real                 zref    (0:2*nz)
471
      real                 deltax,deltay,deltaz
472
      real                 coriol  (0:nx,0:ny)
473
      character*80         pvsrcfile
474
 
475
c     Numerical and physical parameters
476
      real                 eps
477
      parameter            (eps=0.01)
478
 
479
c     Auxiliary variables
480
      integer              cdfid,stat
481
      integer              vardim(4)
482
      real                 misdat
483
      integer              ndimin
484
      real                 varmin(4),varmax(4),stag(4)
485
      integer              i,j,k,nf1
486
      integer              ntimes
487
      real                 time(200)  
488
      character*80         vnam(100),varname
489
      integer              nvars
490
      integer              isok,ierr
491
      real                 x(0:nx,0:ny),y(0:nx,0:ny)
492
      real                 mean,count
493
 
494
c     Get grid description from topography
495
      call cdfopn(pvsrcfile,cdfid,stat)
496
      if (stat.ne.0) goto 997
497
      call getvars(cdfid,nvars,vnam,stat)
498
      if (stat.ne.0) goto 997
499
      isok=0
500
      varname='ORO'
501
      call check_varok(isok,varname,vnam,nvars)
502
      if (isok.eq.0) goto 997
503
      call getdef(cdfid,varname,ndimin,misdat,vardim,
504
     >            varmin,varmax,stag,stat)    
505
      if (stat.ne.0) goto 997
506
      time(1)=0.
507
      call gettimes(cdfid,time,ntimes,stat)  
508
      if (stat.ne.0) goto 997
509
      call clscdf(cdfid,stat)
510
      if (stat.ne.0) goto 997
511
 
512
c     Open output netcdf file
513
      call cdfopn(pvsrcfile,cdfid,stat)
514
      if (stat.ne.0) goto 997
515
 
516
c     Create the variable if necessary
517
      call getvars(cdfid,nvars,vnam,stat)
518
      if (stat.ne.0) goto 997
519
 
520
c     Read data from netcdf file
521
      isok=0
522
      varname='NSQREF'
523
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
524
      call check_varok(isok,varname,vnam,nvars)
525
      if (isok.eq.0) goto 997
526
      call getdat(cdfid,varname,time(1),0,nsqref,stat)
527
      if (stat.ne.0) goto 997
528
 
529
      isok=0
530
      varname='RHOREF'
531
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
532
      call check_varok(isok,varname,vnam,nvars)
533
      if (isok.eq.0) goto 997
534
      call getdat(cdfid,varname,time(1),0,rhoref,stat)
535
      if (stat.ne.0) goto 997
536
 
537
      isok=0
538
      varname='THETAREF'
539
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
540
      call check_varok(isok,varname,vnam,nvars)
541
      if (isok.eq.0) goto 997
542
      call getdat(cdfid,varname,time(1),0,thetaref,stat)
543
      if (stat.ne.0) goto 997
544
 
545
      isok=0
546
      varname='PREREF'
547
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
548
      call check_varok(isok,varname,vnam,nvars)
549
      if (isok.eq.0) goto 997
550
      call getdat(cdfid,varname,time(1),0,pressref,stat)
551
      if (stat.ne.0) goto 997
552
 
553
      isok=0
554
      varname='ZREF'
555
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
556
      call check_varok(isok,varname,vnam,nvars)
557
      if (isok.eq.0) goto 997
558
      call getdat(cdfid,varname,time(1),0,zref,stat)
559
      if (stat.ne.0) goto 997
560
 
561
      isok=0
562
      varname='CORIOL'
563
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
564
      call check_varok(isok,varname,vnam,nvars)
565
      if (isok.eq.0) goto 997
566
      call getdat(cdfid,varname,time(1),0,coriol,stat)
567
      if (stat.ne.0) goto 997
568
 
569
      isok=0
570
      varname='X'
571
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
572
      call check_varok(isok,varname,vnam,nvars)
573
      if (isok.eq.0) goto 997
574
      call getdat(cdfid,varname,time(1),0,x,stat)
575
      if (stat.ne.0) goto 997
576
 
577
      isok=0
578
      varname='Y'
579
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
580
      call check_varok(isok,varname,vnam,nvars)
581
      if (isok.eq.0) goto 997
582
      call getdat(cdfid,varname,time(1),0,y,stat)
583
      if (stat.ne.0) goto 997
584
 
585
c     Close  netcdf file
586
      call clscdf(cdfid,stat)
587
      if (stat.ne.0) goto 997
588
 
589
c     Determine the grid spacings <deltax, deltay, deltaz>
590
      mean=0.
591
      count=0.
592
      do i=1,nx
593
         do j=0,ny
594
            mean=mean+abs(x(i,j)-x(i-1,j))
595
            count=count+1.
596
         enddo
597
      enddo
598
      deltax=mean/count
599
 
600
      mean=0.
601
      count=0.
602
      do j=1,ny
603
         do i=0,nx
604
            mean=mean+abs(y(i,j)-y(i,j-1))
605
            count=count+1.
606
         enddo
607
      enddo
608
      deltay=mean/count
609
 
610
      mean=0.
611
      count=0.
612
      do k=1,nz-1
613
         mean=mean+abs(zref(k+1)-zref(k-1))
614
         count=count+1.
615
      enddo
616
      deltaz=mean/count
617
 
618
      return
619
 
620
c     Exception handling
621
 997  print*,'Read_Ref: Problem with input netcdf file... Stop'
622
      stop
623
 
624
      end
625
 
626
 
627
c     --------------------------------------------------------------------------------
628
c     Get grid parameters
629
c     --------------------------------------------------------------------------------
630
 
631
      subroutine read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv,
632
     >                     pvsrcfile)
633
 
634
c     Get the grid parameters from the variable <THETA> on the input file <pvsrcfile>.
635
c     The grid parameters are
636
c        nx,ny,nz                : Number of grid points in x, y and z direction
637
c        xmin,ymin,zmin          : Minimum domain coordinates in x, y and z direction
638
c        xmax,ymax,zmax          : Maximal domain coordinates in x, y and z direction
639
c        dx,dy,dz                : Horizontal and vertical resolution
640
c     Additionally, it is checked whether the vertical grid is equally spaced. If ok,
641
c     the grid paramters are transformed from lon/lat to distance (in meters)
642
 
643
      implicit none
644
 
645
c     Declaration of subroutine parameters
646
      character*80   pvsrcfile
647
      integer        nx,ny,nz
648
      real           dx,dy,dz
649
      real           xmin,ymin,zmin,xmax,ymax,zmax
650
      real           mdv
651
 
652
c     Numerical epsilon and other physical/geoemtrical parameters
653
      real           eps
654
      parameter      (eps=0.01)
655
 
656
c     Auxiliary variables
657
      integer        cdfid,cstid
658
      integer        ierr
659
      character*80   vnam(100),varname
660
      integer        nvars
661
      integer        isok
662
      integer        vardim(4)
663
      real           misdat
664
      real           varmin(4),varmax(4),stag(4)
665
      real           aklev(1000),bklev(1000),aklay(1000),bklay(1000)
666
      real           dh
667
      character*80   csn
668
      integer        ndim
669
      integer        i
670
 
671
c     Get all grid parameters
672
      call cdfopn(pvsrcfile,cdfid,ierr)
673
      if (ierr.ne.0) goto 998
674
      call getvars(cdfid,nvars,vnam,ierr)
675
      if (ierr.ne.0) goto 998
676
      isok=0
677
      varname='QGPV'
678
      call check_varok(isok,varname,vnam,nvars)
679
      if (isok.eq.0) goto 998
680
      call getcfn(cdfid,csn,ierr)
681
      if (ierr.ne.0) goto 998
682
      call cdfopn(csn,cstid,ierr)
683
      if (ierr.ne.0) goto 998
684
      call getdef(cdfid,varname,ndim,misdat,vardim,varmin,varmax,
685
     >            stag,ierr)
686
      if (ierr.ne.0) goto 998
687
      nx=vardim(1)
688
      ny=vardim(2)
689
      nz=vardim(3)
690
      xmin=varmin(1)
691
      ymin=varmin(2)
692
      zmin=varmin(3)
693
      call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)
694
      if (ierr.ne.0) goto 998
695
      call getgrid(cstid,dx,dy,ierr)
696
      if (ierr.ne.0) goto 998
697
      xmax=varmax(1)
698
      ymax=varmax(2)
699
      zmax=varmax(3)
700
      dz=(zmax-zmin)/real(nz-1)
701
      call clscdf(cstid,ierr)
702
      if (ierr.ne.0) goto 998
703
      call clscdf(cdfid,ierr)
704
      if (ierr.ne.0) goto 998
705
 
706
c     Check whether the grid is equallay spaced in the vertical
707
      do i=1,nz-1
708
         dh=aklev(i+1)-aklev(i)
709
         if (abs(dh-dz).gt.eps) then
710
            print*,'Aklev: Vertical grid must be equally spaced... Stop'
711
            print*,(aklev(i),i=1,nz)
712
            stop
713
         endif
714
         dh=aklay(i+1)-aklay(i)
715
         if (abs(dh-dz).gt.eps) then
716
            print*,'Aklay: Vertical grid must be equally spaced... Stop'
717
            print*,(aklay(i),i=1,nz)
718
            stop
719
         endif
720
      enddo
721
 
722
c     Set missing data value
723
      mdv=misdat
724
 
725
      return
726
 
727
c     Exception handling
728
 998  print*,'Read_Dim: Problem with input netcdf file... Stop'
729
      stop
730
 
731
      end
732
 
733
 
734
c     -------------------------------------------------------------------------------
735
c     Read the input netcdf file
736
c     --------------------------------------------------------------------------------
737
 
738
      SUBROUTINE read_inp (pv,thetabot,thetatop,ua,ub,va,vb,
739
     >                     nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,
740
     >                     pvsrcfile)
741
 
742
c     Read all needed field from netcdf file <pvsrcfile>. The input fields are:
743
c        pv                : quasigeostrophic potential vorticity
744
c        thetabot,thetatop : potential temperature at lower and upper boundary
745
c        ua,ub             : Zonal wind at western and eastern boundary
746
c        va,vb             : Meridional wind at southern and northern boundary
747
c     The grid is specified by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed
748
c     whether the input files are consitent with this grid. The input netcdf file must
749
c     contain the variables <QGPV,THETA,U,V>. If the netcdf file also contains the fields
750
c     <DQGPV,DTHETA,DU,DV>, these increments are added to <QGPV,THETA,U,V>.
751
 
752
      implicit   none
753
 
754
c     Declaration of subroutine parameters
755
      integer              nx,ny,nz
756
      real                 pv(0:nx,0:ny,0:nz)
757
      real                 thetatop(0:nx,0:ny)
758
      real                 thetabot(0:nx,0:ny)
759
      real                 ua(0:nx,0:nz)
760
      real                 ub(0:nx,0:nz)
761
      real                 va(0:ny,0:nz)
762
      real                 vb(0:ny,0:nz)
763
      character*(30)       pvsrcfile
764
      real                 dx,dy,dz,xmin,ymin,zmin
765
 
766
c     Numerical and physical parameters
767
      real                 eps
768
      parameter            (eps=0.01)
769
 
770
c     Auxiliary variables
771
      integer              cdfid,stat
772
      integer              vardim(4)
773
      real                 misdat
774
      real                 varmin(4),varmax(4),stag(4)
775
      integer              ndimin,outid,i,j,k
776
      real                 max_th
777
      real                 tmp(0:nx,0:ny,0:nz)
778
      integer              ntimes
779
      real                 time(200)       
780
      integer              nvars
781
      character*80         vnam(100),varname
782
      integer              isok
783
 
784
c     Open the input netcdf file
785
      call cdfopn(pvsrcfile,cdfid,stat)
786
      if (stat.ne.0) goto 998
787
 
788
c     Check whether needed variables are on file
789
      call getvars(cdfid,nvars,vnam,stat)
790
      if (stat.ne.0) goto 998
791
      isok=0
792
      varname='TH'
793
      call check_varok(isok,varname,vnam,nvars)
794
      varname='U'
795
      call check_varok(isok,varname,vnam,nvars)
796
      varname='V'
797
      call check_varok(isok,varname,vnam,nvars)
798
      varname='QGPV'
799
      call check_varok(isok,varname,vnam,nvars)
800
      if (isok.ne.4) goto 998
801
 
802
c     Get the grid parameters 
803
      varname='QGPV'
804
      call getdef(cdfid,varname,ndimin,misdat,vardim,
805
     >            varmin,varmax,stag,stat)    
806
      if (stat.ne.0) goto 998
807
      time(1)=0.
808
      call gettimes(cdfid,time,ntimes,stat)  
809
      if (stat.ne.0) goto 998
810
 
811
c     Check whether grid parameters are consitent
812
      if ( (vardim(1).ne.nx+1).or.
813
     >     (vardim(2).ne.ny+1).or.
814
     >     (vardim(3).ne.nz+1).or.
815
     >     (abs(varmin(1)-xmin).gt.eps).or.
816
     >     (abs(varmin(2)-ymin).gt.eps).or.
817
     >     (abs(varmin(3)-zmin).gt.eps).or.
818
     >     (abs((varmax(1)-varmin(1))/real(nx)-dx).gt.eps).or.
819
     >     (abs((varmax(2)-varmin(2))/real(ny)-dy).gt.eps).or.
820
     >     (abs((varmax(3)-varmin(3))/real(nz)-dz).gt.eps) ) then
821
         print*,'Input grid inconsitency...'
822
         print*,'xmin : ',xmin,varmin(1)
823
         print*,'ymin : ',ymin,varmin(2)
824
         print*,'zmin : ',zmin,varmin(3)
825
         print*,'dx   : ',dx,(varmax(1)-varmin(1))/real(nx)
826
         print*,'dy   : ',dy,(varmax(2)-varmin(2))/real(ny)
827
         print*,'dz   : ',dz,(varmax(3)-varmin(3))/real(nz)
828
         print*,'nx   : ',nx
829
         print*,'ny   : ',ny
830
         print*,'nz   : ',nz
831
         goto 998
832
      endif
833
 
834
c     THETA: Load upper and lower boundary values
835
      varname='TH'
836
      call getdat(cdfid,varname,time(1),0,tmp,stat)  
837
      print*,'R TH       ',trim(pvsrcfile)
838
      if (stat.ne.0) goto 998
839
      do i=0,nx
840
         do j=0,ny
841
            if ( abs(tmp(i,j,0)-misdat).lt.eps ) then
842
               thetabot(i,j)=0.
843
            else
844
               thetabot(i,j)=tmp(i,j,0)     
845
            endif
846
            if ( abs(tmp(i,j,nz)-misdat).lt.eps ) then
847
               thetatop(i,j)=0.
848
            else
849
               thetatop(i,j)=tmp(i,j,nz)
850
            endif
851
         enddo
852
      enddo
853
 
854
c     U: Load zonal velocity at southern and northern boundary
855
      varname='U'
856
      call getdef(cdfid,varname,ndimin,misdat,vardim,
857
     >            varmin,varmax,stag,stat)
858
      if (stat.ne.0) goto 998
859
      call getdat(cdfid,varname,time(1),0,tmp,stat)
860
      print*,'R U        ',trim(pvsrcfile)
861
      if (stat.ne.0) goto 998
862
      do i=0,nx
863
         do k=0,nz
864
            if ( abs(tmp(i,0,k)-misdat).lt.eps ) then
865
               ua(i,k)=0.
866
            else
867
               ua(i,k)=tmp(i,0,k)
868
            endif
869
            if ( abs(tmp(i,ny,k)-misdat).lt.eps ) then
870
               ub(i,k)=0.
871
            else
872
               ub(i,k)=tmp(i,ny,k)
873
           endif
874
         enddo
875
      enddo
876
 
877
c     Load meridional velocity at western and eastern boundary
878
      varname='V'
879
      call getdef(cdfid,varname,ndimin,misdat,vardim,
880
     >            varmin,varmax,stag,stat)
881
      if (stat.ne.0) goto 998
882
      call getdat(cdfid,varname,time(1),0,tmp,stat)
883
      print*,'R V        ',trim(pvsrcfile)
884
      if (stat.ne.0) goto 998
885
      do j=0,ny
886
         do k=0,nz
887
            if ( abs(tmp(0,j,k)-misdat).lt.eps ) then
888
               va(j,k)=0.
889
            else
890
               va(j,k)=tmp(0,j,k)
891
            endif
892
            if ( abs(tmp(nx,j,k)-misdat).lt.eps ) then
893
               vb(j,k)=0.
894
            else
895
               vb(j,k)=tmp(nx,j,k)
896
            endif
897
         enddo
898
      enddo      
899
 
900
c     Load qgPV
901
      varname='QGPV'
902
      call getdef(cdfid,varname,ndimin,misdat,vardim,
903
     >            varmin,varmax,stag,stat)
904
      if (stat.ne.0) goto 998
905
      call getdat(cdfid,varname,time(1),0,tmp,stat)
906
      print*,'R QGPV     ',trim(pvsrcfile)
907
      if (stat.ne.0) goto 998
908
      do i=0,nx
909
         do j=0,ny
910
            do k=0,nz
911
               if ( abs(tmp(i,j,k)-misdat).lt.eps ) then
912
                  pv(i,j,k)=0.
913
               else
914
                  pv(i,j,k)=tmp(i,j,k)
915
               endif
916
            enddo
917
         enddo
918
      enddo
919
 
920
c     Close input netcdf file
921
      call clscdf(cdfid,stat)
922
      if (stat.ne.0) goto 998
923
 
924
      return
925
 
926
c     Exception handling
927
 998  print*,'Problem with input netcdf file... Stop'
928
      stop
929
 
930
      end
931
 
932
c     --------------------------------------------------------------------------------
933
c     Check whether variable is found on netcdf file
934
c     --------------------------------------------------------------------------------
935
 
936
      subroutine check_varok (isok,varname,varlist,nvars)
937
 
938
c     Check whether the variable <varname> is in the list <varlist(nvars)>. If this is 
939
C     the case, <isok> is incremented by 1. Otherwise <isok> keeps its value.
940
 
941
      implicit none
942
 
943
c     Declaraion of subroutine parameters
944
      integer      isok
945
      integer      nvars
946
      character*80 varname
947
      character*80 varlist(nvars)
948
 
949
c     Auxiliary variables
950
      integer      i
951
 
952
c     Main
953
      do i=1,nvars
954
         if (trim(varname).eq.trim(varlist(i))) isok=isok+1
955
      enddo
956
 
957
      end
958
 
959
 
960
c     ********************************************************************************
961
c     * INVERSION ROUTINES                                                           *
962
c     ********************************************************************************
963
 
964
c     --------------------------------------------------------------------------------
965
c     SOR algorithm (successive over relaxation)
966
c     --------------------------------------------------------------------------------
967
 
968
      SUBROUTINE sor(psi,nsq,thetatop,thetabot,thetaref,rhoref,
969
     >               coriol,pv,ua,ub,va,vb,a,b,c,nx,ny,nz,dx,dy,dz)
970
 
971
c     Solve the qgPV equation by succesive over relaxation (SOR). The subroutine
972
c     parameters are:
973
c
974
c        psi                 : Streamfunction, i.e. result of the PV inversion
975
c        nsq,rhoref,thetaref : Reference profile
976
c        thetatop,thetabot   : Upper and lower boundary condition
977
c        pv                  : quasigeostrophic potential vorticity (qgPV)
978
c        ua,ub,va,vb         : lateral boundary condition for wind
979
c        a,b,c               : Matrices for the inversion operator
980
c        nx,ny,nz,dx,dy,dz   : Grid specification
981
c        coriol              : Coriolis parameter
982
 
983
      implicit none
984
 
985
c     Declaration of subroutine parameters
986
      integer              nx,ny,nz
987
      real                 dx,dy,dz
988
      real                 psi     (0:nx,0:ny,0:nz)
989
      real                 nsq     (0:2*nz)
990
      real                 thetatop(0:nx,0:ny)
991
      real                 thetabot(0:nx,0:ny)
992
      real                 thetaref(0:2*nz)
993
      real                 rhoref  (0:2*nz)
994
      real                 pv      (0:nx,0:ny,0:nz)
995
      real                 ua      (0:nx,0:nz)
996
      real                 ub      (0:nx,0:nz)
997
      real                 va      (0:ny,0:nz)
998
      real                 vb      (0:ny,0:nz)
999
      real                 a       (0:nx,0:ny,0:nz)
1000
      real                 b       (0:nx,0:ny,0:nz)
1001
      real                 c       (0:nx,0:ny,0:nz)
1002
      real                 coriol  (0:nx,0:ny)
1003
 
1004
c     Numerical and physical parameters
1005
      real                 maxspec
1006
      parameter            (maxspec=2.0)
1007
      integer              nofiter
1008
      parameter            (nofiter=500)
1009
      real                 omega
1010
      parameter            (omega=1.81)
1011
 
1012
c     Auxiliary variables
1013
      integer              counter
1014
      integer              i,j,k
1015
      real                 deltasq,psigauge
1016
      real                 specx,specy,specz
1017
      real                 helpx,helpy,helpz
1018
 
1019
c     Init the output array
1020
      do i=0,nx
1021
         do j=0,ny
1022
            do k=0,nz
1023
               psi(i,j,k)=0.
1024
            enddo
1025
         enddo
1026
      enddo
1027
 
1028
c     Calculate the spectrum of the matrix 
1029
      i=nx/2
1030
      specx=4.*a(i,0,0)/
1031
     >      (2.*a(i,0,0)+b(i,0,0)+b(i,1,0)+c(i,0,0)+c(i,0,1))
1032
      specy=2.*(b(i,0,0)+b(i,1,0))/
1033
     >      (2.*a(i,0,0)+b(i,0,0)+b(i,1,0)+c(i,0,0)+c(i,0,1))
1034
      specz=2.*(c(i,0,0)+c(i,0,1))/
1035
     >      (2.*a(i,0,0)+b(i,0,0)+b(i,1,0)+c(i,0,0)+c(i,0,1))
1036
      do k=1,nz-2
1037
         do j=1,ny-2
1038
 
1039
            helpx=4.*a(i,j,k)/
1040
     >            (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k)+c(i,j,k-1))
1041
            if (helpx.gt.specx) specx=helpx
1042
 
1043
            helpy=2.*(b(i,j,k)+b(i,j+1,k))/
1044
     >            (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k)+c(i,j,k-1))
1045
            if (helpy.gt.specy) specy=helpy
1046
 
1047
            helpz=2.*(c(i,j,k)+c(i,j,k+1))/
1048
     >            (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k)+c(i,j,k-1))
1049
            if (helpz.gt.specz) specz=helpz
1050
 
1051
         enddo
1052
      enddo
1053
 
1054
c     Check whether the dimensions of the grid are sufficient
1055
      print *
1056
      print *, 'Spectrum of the matrix in each direction '
1057
      print *, 'Spectrum = ', specx, specy, specz
1058
      print *
1059
      if ((maxspec*specx.lt.specy).or.(maxspec*specx.lt.specz)) then
1060
         print*,' Nx too small... Stop'
1061
         stop
1062
      endif
1063
      if ((maxspec*specy.lt.specx).or.(maxspec*specy.lt.specz)) then
1064
         print*,'Ny too small... Stop'
1065
         stop
1066
      endif
1067
      if ((maxspec*specz.lt.specx).or.(maxspec*specz.lt.specy)) then
1068
         print*,'Nz too small... Stop'
1069
         stop
1070
      endif
1071
 
1072
c     Calculate error: control variable for the iteration
1073
      psigauge=0.
1074
      deltasq=0.
1075
      do k=1,nz-1
1076
         do i=1,nx-1
1077
            do j=1,ny-1
1078
               deltasq=deltasq+(-pv(i,j,k)+(
1079
     >              a(i,j,k)*(psi(i+1,j,k)+psi(i-1,j,k)-
1080
     >              2.*psi(i,j,k)) +
1081
     >              b(i,j,k)*(psi(i,j+1,k)-psi(i,j,k))-
1082
     >              b(i,j-1,k)*(psi(i,j,k)-psi(i,j-1,k))+ 
1083
     >              c(i,j,k)*(psi(i,j,k+1)-psi(i,j,k))-
1084
     >              c(i,j,k-1)*(psi(i,j,k)-psi(i,j,k-1))
1085
     >              )/(dx*dy*dz*rhoref(2*k)))**2.
1086
 
1087
            enddo
1088
         enddo
1089
      enddo
1090
      print 102, 'psigauge', psigauge, 'deltasq',
1091
     >           deltasq/(real(nx)*real(ny)*real(nz))
1092
 
1093
c     Iterations
1094
      do counter=1,nofiter
1095
 
1096
C        Perform one iteration step
1097
         call psiappsor(omega,pv,psi,nsq,rhoref,thetatop,
1098
     >                  thetabot,thetaref,coriol,ua,ub,va,vb,
1099
     >                  a,b,c,nx,ny,nz,dx,dy,dz)
1100
 
1101
 
1102
c        Adjustment
1103
         if (mod(counter,100).eq.0) then
1104
            psigauge=0.
1105
            do i=0,nx
1106
               do j=0,ny
1107
                  if (psi(i,j,0).lt.psigauge) then 
1108
                     psigauge=psi(i,j,0)
1109
                  endif
1110
               enddo
1111
            enddo
1112
            do k=0,nz
1113
               do i=0,nx
1114
                  do j=0,ny
1115
                     psi(i,j,k)=psi(i,j,k)-psigauge
1116
                  enddo
1117
               enddo
1118
            enddo
1119
         endif
1120
 
1121
c        Calculate error: control variable for the iteration        
1122
         if (mod(counter,nofiter/10).eq.0) then
1123
            deltasq=0.
1124
            do k=1,nz-1
1125
               do i=1,nx-1
1126
                  do j=1,ny-1
1127
                     deltasq=deltasq+(-pv(i,j,k)+(
1128
     >                 a(i,j,k)*(psi(i+1,j,k)+psi(i-1,j,k)-
1129
     >                    2.*psi(i,j,k)) +
1130
     >                 b(i,j,k)*(psi(i,j+1,k)-psi(i,j,k))-
1131
     >                 b(i,j-1,k)*(psi(i,j,k)-psi(i,j-1,k))+ 
1132
     >                 c(i,j,k)*(psi(i,j,k+1)-psi(i,j,k))-
1133
     >                 c(i,j,k-1)*(psi(i,j,k)-psi(i,j,k-1))
1134
     >                    )/(dx*dy*dz*rhoref(2*k)))**2.
1135
                  enddo
1136
               enddo
1137
            enddo
1138
            print 102, 'psigauge', psigauge, 'deltasq',
1139
     >           deltasq/(real(nx)*real(ny)*real(nz))
1140
         endif
1141
 
1142
      enddo
1143
 
1144
      return
1145
 
1146
c     Format specifications
1147
 102  format (a11, ' = ',e10.3,a11, ' = ',e10.3)
1148
 
1149
      end
1150
 
1151
c     --------------------------------------------------------------------------------
1152
c     SOR algorithm (successive over relaxation)
1153
c     --------------------------------------------------------------------------------
1154
 
1155
      subroutine psiappsor(omega,pv,psi,nsq,rhoref,thetatop,
1156
     >                     thetabot,thetaref,coriol,ua,ub,va,vb,
1157
     >                     a,b,c,nx,ny,nz,dx,dy,dz)
1158
 
1159
c     Perform one relaxation step
1160
c
1161
c        psi                 : Streamfunction, i.e. result of the PV inversion
1162
c        nsq,rhoref,thetaref : Reference profile
1163
c        thetatop,thetabot   : Upper and lower boundary condition
1164
c        pv                  : quasigeostrophic potential vorticity (qgPV)
1165
c        ua,ub,va,vb         : lateral boundary condition for wind
1166
c        a,b,c               : Matrices for the inversion operator
1167
c        nx,ny,nz,dx,dy,dz   : Grid specification
1168
c        nofiter             : Number of iterations
1169
c        omega               : Relaxation parameter
1170
c        coriol              : Coriolis parameter
1171
 
1172
      implicit none
1173
 
1174
c     Declaration of subroutine parameters
1175
      integer              nx,ny,nz
1176
      real                 pv(0:nx,0:ny,0:nz)
1177
      real                 psi(0:nx,0:ny,0:nz)
1178
      real                 nsq(0:2*nz)
1179
      real                 rhoref(0:2*nz)
1180
      real                 thetatop(0:nx,0:ny)
1181
      real                 thetabot(0:nx,0:ny)
1182
      real                 thetaref(0:2*nz)
1183
      real                 ua(0:nx,0:nz)
1184
      real                 ub(0:nx,0:nz)
1185
      real                 va(0:ny,0:nz)
1186
      real                 vb(0:ny,0:nz)
1187
      real                 a(0:nx,0:ny,0:nz)
1188
      real                 b(0:nx,0:ny,0:nz)
1189
      real                 c(0:nx,0:ny,0:nz)
1190
      real                 coriol(0:nx,0:ny)
1191
      real                 dx,dy,dz
1192
      real                 omega
1193
 
1194
c     Numerical and physical parameters
1195
      real                 g
1196
      parameter            (g=9.81)
1197
 
1198
c     Auxiliary variables
1199
      integer              i,j,k
1200
      real                 dxy,dxz,dyz,dxyz
1201
 
1202
c     Set the area and volume infinitesimals for integration
1203
      dxy=dx*dy
1204
      dxz=dx*dz
1205
      dyz=dy*dz
1206
      dxyz=dx*dy*dz
1207
 
1208
c     Inner
1209
      do k=1,nz-1
1210
         do i=1,nx-1
1211
            do j=1,ny-1
1212
               psi(i,j,k)=omega*(-dxyz*
1213
     >                    rhoref(2*k)*pv(i,j,k)+a(i,j,k)*
1214
     >                    (psi(i+1,j,k)+psi(i-1,j,k))+b(i,j,k)*
1215
     >                    psi(i,j+1,k)+b(i,j-1,k)*psi(i,j-1,k)+c(i,j,k)*
1216
     >                    psi(i,j,k+1)+c(i,j,k-1)*psi(i,j,k-1))/
1217
     >                    (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k-1)+
1218
     >                    c(i,j,k))+(1.-omega)*psi(i,j,k)
1219
            enddo
1220
         enddo
1221
      enddo
1222
 
1223
c     ZY plane
1224
      do k=1,nz-1
1225
         do j=1,ny-1
1226
            psi(0,j,k)=omega*(-dyz*
1227
     >           rhoref(2*k)*(dx*pv(0,j,k)+va(j,k))+
1228
     >           a(0,j,k)*psi(1,j,k)+
1229
     >           b(0,j,k)*psi(0,j+1,k)+b(0,j-1,k)*psi(0,j-1,k)+
1230
     >           c(0,j,k)*psi(0,j,k+1)+c(0,j,k-1)*psi(0,j,k-1))/
1231
     >           (a(0,j,k)+b(0,j,k)+b(0,j-1,k)+c(0,j,k-1)+c(0,j,k))
1232
     >           +(1.-omega)*psi(0,j,k)
1233
c     
1234
            psi(nx,j,k)=omega*(-dyz*
1235
     >           rhoref(2*k)*(dx*pv(nx,j,k)-vb(j,k))+
1236
     >           a(nx,j,k)*psi(nx-1,j,k)+
1237
     >           b(nx,j,k)*psi(nx,j+1,k)+b(nx,j-1,k)*psi(nx,j-1,k)+
1238
     >           c(nx,j,k)*psi(nx,j,k+1)+c(nx,j,k-1)*psi(nx,j,k-1))/
1239
     >           (a(nx,j,k)+b(nx,j-1,k)+b(nx,j,k)+c(nx,j,k-1)+c(nx,j,k))
1240
     >           +(1.-omega)*psi(nx,j,k)
1241
         enddo
1242
      enddo
1243
 
1244
c     ZX plane
1245
      do k=1,nz-1
1246
         do i=1,nx-1
1247
            psi(i,0,k)=omega*(-dxz*
1248
     >           rhoref(2*k)*(dy*pv(i,0,k)-ua(i,k))+
1249
     >           a(i,0,k)*(psi(i+1,0,k)+psi(i-1,0,k))+
1250
     >           b(i,0,k)*psi(i,1,k)+
1251
     >           c(i,0,k)*psi(i,0,k+1)+c(i,0,k-1)*psi(i,0,k-1))/
1252
     >           (2.*a(i,0,k)+b(i,0,k)+c(i,0,k-1)+c(i,0,k))
1253
     >           +(1.-omega)*psi(i,0,k)
1254
c     
1255
            psi(i,ny,k)=omega*(-dxz*
1256
     >           rhoref(2*k)*(dy*pv(i,ny,k)+ub(i,k))+
1257
     >           a(i,ny-1,k)*(psi(i+1,ny,k)+psi(i-1,ny,k))+
1258
     >           b(i,ny-1,k)*psi(i,ny-1,k)+
1259
     >           c(i,ny-1,k)*psi(i,ny,k+1)+c(i,ny-1,k-1)*
1260
     >           psi(i,ny,k-1))/(2.*a(i,ny-1,k)+b(i,ny-1,k)+
1261
     >           c(i,ny-1,k-1)+c(i,ny-1,k))
1262
     >           +(1.-omega)*psi(i,ny,k)
1263
         enddo
1264
      enddo
1265
 
1266
c     XY plane
1267
      do i=1,nx-1
1268
         do j=1,ny-1
1269
            psi(i,j,0)=omega*(-dxy*rhoref(0)*(
1270
     >           dz*pv(i,j,0)+g*coriol(i,j)*thetabot(i,j)/
1271
     >           (nsq(0)*thetaref(0)))+
1272
     >           a(i,j,0)*(psi(i+1,j,0)+psi(i-1,j,0))+
1273
     >           b(i,j,0)*psi(i,j+1,0)+b(i,j-1,0)*psi(i,j-1,0)+
1274
     >           c(i,j,0)*psi(i,j,1))/
1275
     >           (2.*a(i,j,0)+b(i,j-1,0)+b(i,j,0)+c(i,j,0))
1276
     >           +(1.-omega)*psi(i,j,0)
1277
c     
1278
            psi(i,j,nz)=omega*(-dxy*rhoref(2*nz)*(
1279
     >           dz*pv(i,j,nz)-g*coriol(i,j)*thetatop(i,j)/
1280
     >           (nsq(2*nz)*thetaref(2*nz)))+
1281
     >           a(i,j,nz)*(psi(i+1,j,nz)+psi(i-1,j,nz))+
1282
     >           b(i,j,nz)*psi(i,j+1,nz)+b(i,j-1,nz)*psi(i,j-1,nz)+
1283
     >           c(i,j,nz-1)*psi(i,j,nz-1))/
1284
     >           (2.*a(i,j,nz)+b(i,j-1,nz)+b(i,j,nz)+c(i,j,nz-1))
1285
     >           +(1.-omega)*psi(i,j,nz)
1286
         enddo
1287
      enddo
1288
 
1289
c     Y edges
1290
      do j=1,ny-1
1291
         psi(0,j,0)=omega*(-dy*rhoref(0)*(dxz*pv(0,j,0)+
1292
     >        dz*va(j,0)+dx*g*coriol(0,j)*thetabot(0,j)/
1293
     >        (nsq(0)*thetaref(0)))+
1294
     >        a(0,j,0)*psi(1,j,0)+
1295
     >        b(0,j,0)*psi(0,j+1,0)+b(0,j-1,0)*psi(0,j-1,0)+
1296
     >        c(0,j,0)*psi(0,j,1))/
1297
     >        (a(0,j,0)+b(0,j-1,0)+b(0,j,0)+c(0,j,0))
1298
     >        +(1.-omega)*psi(0,j,0)
1299
c     
1300
         psi(nx,j,0)=omega*(-dy*rhoref(0)*(dxz*pv(nx,j,0)-
1301
     >        dz*vb(j,0)+dx*g*coriol(nx,j)*thetabot(nx,j)/
1302
     >        (nsq(0)*thetaref(0)))+
1303
     >        a(nx,j,0)*psi(nx-1,j,0)+
1304
     >        b(nx,j,0)*psi(nx,j+1,0)+b(nx,j-1,0)*psi(nx,j-1,0)+
1305
     >        c(nx,j,0)*psi(nx,j,1))/
1306
     >        (a(nx,j,0)+b(nx,j-1,0)+b(nx,j,0)+c(nx,j,0))
1307
     >        +(1.-omega)*psi(nx,j,0)
1308
c
1309
         psi(0,j,nz)=omega*(-dy*rhoref(2*nz)*(dxz*pv(0,j,nz)+
1310
     >        dz*va(j,nz)-dx*g*coriol(0,j)*thetatop(0,j)/
1311
     >        (nsq(2*nz)*thetaref(2*nz)))+
1312
     >        a(0,j,nz)*psi(1,j,nz)+
1313
     >        b(0,j,nz)*psi(0,j+1,nz)+b(0,j-1,nz)*psi(0,j-1,nz)+
1314
     >        c(0,j,nz-1)*psi(0,j,nz-1))/
1315
     >        (a(0,j,nz)+b(0,j-1,nz)+b(0,j,nz)+c(0,j,nz-1))
1316
     >        +(1.-omega)*psi(0,j,nz)
1317
c     
1318
         psi(nx,j,nz)=omega*(-dy*rhoref(2*nz)*(dxz*pv(nx,j,nz)-
1319
     >        dz*vb(j,nz)-dx*g*coriol(nx,j)*thetatop(nx,j)/
1320
     >        (nsq(2*nz)*thetaref(2*nz)))+
1321
     >        a(nx,j,nz)*psi(nx-1,j,nz)+
1322
     >        b(nx,j,nz)*psi(nx,j+1,nz)+b(nx,j-1,nz)*psi(nx,j-1,nz)+
1323
     >        c(nx,j,nz-1)*psi(nx,j,nz-1))/
1324
     >        (a(nx,j,nz)+b(nx,j-1,nz)+b(nx,j,nz)+c(nx,j,nz-1))
1325
     >        +(1.-omega)*psi(nx,j,nz)
1326
      enddo
1327
 
1328
c     X edges
1329
      do i=1,nx-1
1330
         psi(i,0,0)=omega*(-dx*rhoref(0)*(dyz*pv(i,0,0)-
1331
     >        dz*ua(i,0)+dy*g*coriol(i,0)*thetabot(i,0)/
1332
     >        (nsq(0)*thetaref(0)))+
1333
     >        a(i,0,0)*(psi(i+1,0,0)+psi(i-1,0,0))+
1334
     >        b(i,0,0)*psi(i,1,0)+
1335
     >        c(i,0,0)*psi(i,0,1))/
1336
     >        (2.*a(i,0,0)+b(i,0,0)+c(i,0,0))
1337
     >        +(1.-omega)*psi(i,0,0)
1338
c     
1339
         psi(i,ny,0)=omega*(-dx*rhoref(0)*(dyz*pv(i,ny,0)+
1340
     >        dz*ub(i,0)+dy*g*coriol(i,ny)*thetabot(i,ny)/
1341
     >        (nsq(0)*thetaref(0)))+
1342
     >        a(i,ny,0)*(psi(i+1,ny,0)+psi(i-1,ny,0))+
1343
     >        b(i,ny-1,0)*psi(i,ny-1,0)+
1344
     >        c(i,ny,0)*psi(i,ny,1))/
1345
     >        (2.*a(i,ny,0)+b(i,ny-1,0)+c(i,ny,0))
1346
     >        +(1.-omega)*psi(i,ny,0)
1347
c     
1348
         psi(i,0,nz)=omega*(-dx*rhoref(2*nz)*(dyz*pv(i,0,nz)-
1349
     >        dz*ua(i,nz)-dy*g*coriol(i,0)*thetatop(i,0)/
1350
     >        (nsq(2*nz)*thetaref(2*nz)))+
1351
     >        a(i,0,nz)*(psi(i+1,0,nz)+psi(i-1,0,nz))+
1352
     >        b(i,0,nz)*psi(i,1,nz)+
1353
     >        c(i,0,nz-1)*psi(i,0,nz-1))/
1354
     >        (2.*a(i,0,nz)+b(i,0,nz)+c(i,0,nz-1))
1355
     >        +(1.-omega)*psi(i,0,nz)
1356
c     
1357
         psi(i,ny,nz)=omega*(-dx*rhoref(2*nz)*(dyz*pv(i,ny,nz)+
1358
     >        dz*ub(i,nz)-dy*g*coriol(i,ny)*thetatop(i,ny)/
1359
     >        (nsq(2*nz)*thetaref(2*nz)))+
1360
     >        a(i,ny,nz)*(psi(i+1,ny,nz)+psi(i-1,ny,nz))+
1361
     >        b(i,ny-1,nz)*psi(i,ny-1,nz)+
1362
     >        c(i,ny,nz-1)*psi(i,ny,nz-1))/
1363
     >        (2.*a(i,ny,nz)+b(i,ny-1,nz)+c(i,ny,nz-1))
1364
     >        +(1.-omega)*psi(i,ny,nz)
1365
      enddo
1366
 
1367
c     Z edges
1368
      do k=1,nz-1
1369
         psi(0,0,k)=omega*(-dz*rhoref(2*k)*(dxy*pv(0,0,k)+
1370
     >        dy*va(0,k)-dx*ua(0,k))+
1371
     >        a(0,0,k)*psi(1,0,k)+
1372
     >        b(0,0,k)*psi(0,1,k)+
1373
     >        c(0,0,k)*psi(0,0,k+1)+c(0,0,k-1)*psi(0,0,k-1))/
1374
     >        (a(0,0,k)+b(0,0,k)+c(0,0,k-1)+c(0,0,k))
1375
     >        +(1.-omega)*psi(0,0,k)
1376
c
1377
         psi(nx,0,k)=omega*(-dz*rhoref(2*k)*(dxy*pv(nx,0,k)-
1378
     >        dy*vb(0,k)-dx*ua(nx,k))+
1379
     >        a(nx,0,k)*psi(nx-1,0,k)+
1380
     >        b(nx,0,k)*psi(nx,1,k)+
1381
     >        c(nx,0,k)*psi(nx,0,k+1)+c(nx,0,k-1)*psi(nx,0,k-1))/
1382
     >        (a(nx,0,k)+b(nx,0,k)+c(nx,0,k-1)+c(nx,0,k))
1383
     >        +(1.-omega)*psi(nx,0,k)
1384
c
1385
         psi(0,ny,k)=omega*(-dz*rhoref(2*k)*(dxy*pv(0,ny,k)+
1386
     >        dy*va(ny,k)+dx*ub(0,k))+
1387
     >        a(0,ny,k)*psi(1,ny,k)+
1388
     >        b(0,ny-1,k)*psi(0,ny-1,k)+
1389
     >        c(0,ny,k)*psi(0,ny,k+1)+c(0,ny,k-1)*psi(0,ny,k-1))/
1390
     >        (a(0,ny,k)+b(0,ny-1,k)+c(0,ny,k-1)+c(0,ny,k))
1391
     >        +(1.-omega)*psi(0,ny,k)
1392
c
1393
         psi(nx,ny,k)=omega*(-dz*rhoref(2*k)*(dxy*pv(nx,ny,k)-
1394
     >        dy*vb(ny,k)+dx*ub(nx,k))+
1395
     >        a(nx,ny,k)*psi(nx-1,ny,k)+
1396
     >        b(nx,ny-1,k)*psi(nx,ny-1,k)+
1397
     >        c(nx,ny,k)*psi(nx,ny,k+1)+c(nx,ny,k-1)*psi(nx,ny,k-1))/
1398
     >        (a(nx,ny,k)+b(nx,ny-1,k)+c(nx,ny,k-1)+c(nx,ny,k))
1399
     >        +(1.-omega)*psi(nx,ny,k)
1400
      enddo
1401
 
1402
c     Points
1403
      psi(0,0,0)=omega*(-rhoref(0)*(dxyz*pv(0,0,0)+dyz*va(0,0)-
1404
     >        dxz*ua(0,0)+dxy*g*coriol(0,0)*thetabot(0,0)/
1405
     >        (nsq(0)*thetaref(0)))+
1406
     >        a(0,0,0)*psi(1,0,0)+
1407
     >        b(0,0,0)*psi(0,1,0)+
1408
     >        c(0,0,0)*psi(0,0,1))/
1409
     >        (a(0,0,0)+b(0,0,0)+c(0,0,0))+
1410
     >        (1.-omega)*psi(0,0,0)
1411
c
1412
      psi(nx,0,0)=omega*(-rhoref(0)*(dxyz*pv(nx,0,0)-dyz*vb(0,0)-
1413
     >        dxz*ua(nx,0)+dxy*g*coriol(nx,0)*thetabot(nx,0)/
1414
     >        (nsq(0)*thetaref(0)))+
1415
     >        a(nx,0,0)*psi(nx-1,0,0)+
1416
     >        b(nx,0,0)*psi(nx,1,0)+
1417
     >        c(nx,0,0)*psi(nx,0,1))/
1418
     >        (a(nx,0,0)+b(nx,0,0)+c(nx,0,0))+
1419
     >        (1.-omega)*psi(nx,0,0)
1420
c
1421
      psi(0,ny,0)=omega*(-rhoref(0)*(dxyz*pv(0,ny,0)+dyz*va(ny,0)+
1422
     >        dxz*ub(0,0)+dxy*g*coriol(0,ny)*thetabot(0,ny)/
1423
     >        (nsq(0)*thetaref(0)))+
1424
     >        a(0,ny,0)*psi(1,ny,0)+
1425
     >        b(0,ny-1,0)*psi(0,ny-1,0)+
1426
     >        c(0,ny,0)*psi(0,ny,1))/
1427
     >        (a(0,ny,0)+b(0,ny-1,0)+c(0,ny,0))+
1428
     >        (1.-omega)*psi(0,ny,0)
1429
c
1430
      psi(nx,ny,0)=omega*(-rhoref(0)*(dxyz*pv(nx,ny,0)-dyz*vb(ny,0)+
1431
     >        dxz*ub(nx,0)+dxy*g*coriol(nx,ny)*thetabot(nx,ny)/
1432
     >        (nsq(0)*thetaref(0)))+
1433
     >        a(nx,ny,0)*psi(nx-1,ny,0)+
1434
     >        b(nx,ny-1,0)*psi(nx,ny-1,0)+
1435
     >        c(nx,ny,0)*psi(nx,ny,1))/
1436
     >        (a(nx,ny,0)+b(nx,ny-1,0)+c(nx,ny,0))+
1437
     >        (1.-omega)*psi(nx,ny,0)
1438
c
1439
      psi(0,0,nz)=omega*(-rhoref(2*nz)*(dxyz*pv(0,0,nz)+dyz*va(0,nz)-
1440
     >        dxz*ua(0,nz)-dxy*g*coriol(0,0)*thetatop(0,0)/
1441
     >        (nsq(2*nz)*thetaref(2*nz)))+
1442
     >        a(0,0,nz)*psi(1,0,nz)+
1443
     >        b(0,0,nz)*psi(0,1,nz)+
1444
     >        c(0,0,nz-1)*psi(0,0,nz-1))/
1445
     >       (a(0,0,nz)+b(0,0,nz)+c(0,0,nz-1))+
1446
     >       (1.-omega)*psi(0,0,nz)
1447
c
1448
      psi(nx,0,nz)=omega*(-rhoref(2*nz)*(dxyz*pv(nx,0,nz)-dyz*vb(0,nz)
1449
     >        -dxz*ua(nx,nz)-dxy*g*coriol(nx,0)*thetatop(nx,0)/
1450
     >        (nsq(2*nz)*thetaref(2*nz)))+
1451
     >        a(nx,0,nz)*psi(nx-1,0,nz)+
1452
     >        b(nx,0,nz)*psi(nx,1,nz)+
1453
     >        c(nx,0,nz-1)*psi(nx,0,nz-1))/
1454
     >        (a(nx,0,nz)+b(nx,0,nz)+c(nx,0,nz-1))+
1455
     >        (1.-omega)*psi(nx,0,nz)
1456
c
1457
      psi(0,ny,nz)=omega*(-rhoref(2*nz)*(dxyz*pv(0,ny,nz)+
1458
     >        dyz*va(ny,nz)+dxz*ub(0,nz)-
1459
     >        dxy*g*coriol(0,ny)*thetatop(0,ny)/
1460
     >        (nsq(2*nz)*thetaref(2*nz)))+
1461
     >        a(0,ny,nz)*psi(1,ny,nz)+
1462
     >        b(0,ny-1,nz)*psi(0,ny-1,nz)+
1463
     >        c(0,ny,nz-1)*psi(0,ny,nz-1))/
1464
     >        (a(0,ny,nz)+b(0,ny-1,nz)+c(0,ny,nz-1))+
1465
     >        (1.-omega)*psi(0,ny,nz)
1466
c
1467
      psi(nx,ny,nz)=omega*(-rhoref(2*nz)*(dxyz*pv(nx,ny,nz)-
1468
     >        dyz*vb(ny,nz)+dxz*ub(nx,nz)-
1469
     >        dxy*g*coriol(nx,ny)*thetatop(nx,ny)/
1470
     >        (nsq(2*nz)*thetaref(2*nz)))+
1471
     >        a(nx,ny,nz)*psi(nx-1,ny,nz)+
1472
     >        b(nx,ny-1,nz)*psi(nx,ny-1,nz)+
1473
     >        c(nx,ny,nz-1)*psi(nx,ny,nz-1))/
1474
     >        (a(nx,ny,nz)+b(nx,ny-1,nz)+c(nx,ny,nz-1))+
1475
     >        (1.-omega)*psi(nx,ny,nz)
1476
c
1477
 
1478
      end
1479
 
1480
c     --------------------------------------------------------------------------------
1481
c     Init matrix elements for the inversion
1482
c     --------------------------------------------------------------------------------
1483
 
1484
      subroutine matrixel(a,b,c,coriol,nx,ny,nz,nsq,rhoref,dx,dy,dz)
1485
 
1486
c     Define the coefficients for the inversion problem (see page 119ff in Rene's 
1487
c     dissertation).
1488
 
1489
      implicit none
1490
 
1491
c     Declaration of subroutine parameters
1492
      integer   nx,nz,ny
1493
      real      a     (0:nx,0:ny,0:nz)
1494
      real      b     (0:nx,0:ny,0:nz)
1495
      real      c     (0:nx,0:ny,0:nz)
1496
      real      coriol(0:nx,0:ny)
1497
      real      nsq   (0:2*nz)
1498
      real      rhoref(0:2*nz)
1499
      real      dx,dy,dz
1500
 
1501
c     Auxiliary variables
1502
      integer   i,j,k
1503
 
1504
c     Calculate coefficients
1505
      do i=0,nx
1506
         do j=0,ny
1507
            do k=0,nz
1508
               a(i,j,k)=dy*dz*rhoref(2*k)/(dx*coriol(i,j))
1509
               if (j.lt.ny) then
1510
                  b(i,j,k)=dx*dz*rhoref(2*k)/
1511
     >                     (dy*0.5*(coriol(i,j)+coriol(i,j+1)))
1512
               else
1513
                  b(i,j,k)=0.
1514
               endif
1515
               if (k.lt.nz) then
1516
                  c(i,j,k)=dx*dy*rhoref(2*k+1)*coriol(i,j)/
1517
     >                     (dz*nsq(2*k+1))
1518
               else
1519
                  c(i,j,k)=0.
1520
               endif
1521
            enddo
1522
         enddo
1523
      enddo
1524
 
1525
      end
1526
 
1527
 
1528
 
1529