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