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 prep_iteration
2
 
3
c     ************************************************************************
4
c     * Prepare the next step for the PV inversion                           *
5
c     * Michael Sprenger / Summer, Autumn 2006                               *
6
c     ************************************************************************
7
 
8
      implicit none
9
 
10
c     ------------------------------------------------------------------------
11
c     Declaration of variables and parameters
12
c     ------------------------------------------------------------------------
13
 
14
c     Input and output file
15
      character*80   anomafile
16
      character*80   iterafile
17
 
18
c     Grid parameters
19
      integer        nx,ny,nz
20
      real           xmin,ymin,zmin,xmax,ymax,zmax
21
      real           dx,dy,dz
22
      real           mdv
23
 
24
c     Numerical epsilon and other variables
25
      real           eps
26
      parameter      (eps=0.01)
27
      real           alpha
28
 
29
c     3d arrays
30
      real,allocatable,dimension (:,:,:) :: v_iter,v_anom
31
      real,allocatable,dimension (:,:,:) :: u_iter,u_anom
32
      real,allocatable,dimension (:,:,:) :: t_iter,t_anom
33
      real,allocatable,dimension (:,:,:) :: p_iter,p_anom
34
 
35
c     Auciliary variables 
36
      integer      i,j,k
37
      integer      stat
38
      character*80 varname
39
      character*80 name
40
 
41
c     --------------------------------------------------------------------------------
42
c     Input
43
c     --------------------------------------------------------------------------------
44
 
45
      print*,'********************************************************'
46
      print*,'* PREP_ITERATION                                       *'
47
      print*,'********************************************************'
48
 
49
c     Read parameter file
50
      open(10,file='fort.10')
51
       read(10,*) iterafile
52
       read(10,*) anomafile
53
       read(10,*) name,alpha
54
       if ( trim(name).ne.'ALPHA' ) stop
55
      close(10) 
56
      print*
57
      print*,trim(anomafile)
58
      print*,trim(iterafile)
59
      print*
60
 
61
c     Get lat/lon gid parameters from input file
62
      call read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv,
63
     >               anomafile)
64
      print*,'Read_Dim: nx,ny,nz         = ',nx,ny,nz
65
      print*,'          dx,dy,dz         = ',dx,dy,dz
66
      print*,'          xmin,ymin,zmin   = ',xmin,ymin,zmin
67
      print*,'          mdv              = ',mdv
68
      print*
69
 
70
c     Count from 0, not from 1: consistent with <inv_cart.f>.
71
      nx=nx-1
72
      ny=ny-1
73
      nz=nz-1
74
 
75
c     Allocate memory for 3d arrays 
76
      allocate(u_anom (0:nx,0:ny,0:nz),STAT=stat)
77
      if (stat.ne.0) print*,'error allocating u_anom'
78
      allocate(v_anom (0:nx,0:ny,0:nz),STAT=stat)
79
      if (stat.ne.0) print*,'error allocating v_anom'
80
      allocate(t_anom  (0:nx,0:ny,0:nz),STAT=stat)
81
      if (stat.ne.0) print*,'error allocating t_anom'
82
      allocate(p_anom  (0:nx,0:ny,0:nz),STAT=stat)
83
      if (stat.ne.0) print*,'error allocating p_anom'
84
 
85
      allocate(u_iter (0:nx,0:ny,0:nz),STAT=stat)
86
      if (stat.ne.0) print*,'error allocating u_iter'
87
      allocate(v_iter (0:nx,0:ny,0:nz),STAT=stat)
88
      if (stat.ne.0) print*,'error allocating v_iter'
89
      allocate(t_iter  (0:nx,0:ny,0:nz),STAT=stat)
90
      if (stat.ne.0) print*,'error allocating t_iter'
91
      allocate(p_iter  (0:nx,0:ny,0:nz),STAT=stat)
92
      if (stat.ne.0) print*,'error allocating p_iter'
93
 
94
c     Read anomaly and iteration fields
95
      varname='U'
96
      call read_inp (u_anom,varname,anomafile,
97
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
98
      varname='V'
99
      call read_inp (v_anom,varname,anomafile,
100
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
101
      varname='T'
102
      call read_inp (t_anom,varname,anomafile,
103
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
104
      varname='P'
105
      call read_inp (p_anom,varname,anomafile,
106
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
107
      varname='U'
108
      call read_inp (u_iter,varname,iterafile,
109
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
110
      varname='V'
111
      call read_inp (v_iter,varname,iterafile,
112
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
113
      varname='T'
114
      call read_inp (t_iter,varname,iterafile,
115
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
116
      varname='P'
117
      call read_inp (p_iter,varname,iterafile,
118
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
119
 
120
c     --------------------------------------------------------------------------------
121
c     Modify the iteration fields
122
c     --------------------------------------------------------------------------------
123
 
124
      do i=0,nx
125
         do j=0,ny
126
            do k=0,nz
127
 
128
c              Update zonal velocity
129
               if ((abs(u_anom(i,j,k)-mdv).gt.eps).and.
130
     >             (abs(u_iter(i,j,k)-mdv).gt.eps)) then
131
                  u_iter(i,j,k)=u_iter(i,j,k)-alpha*u_anom(i,j,k)
132
               else
133
                  u_iter(i,j,k)=mdv
134
               endif
135
 
136
c              Update meridional velocity
137
               if ((abs(v_anom(i,j,k)-mdv).gt.eps).and.
138
     >             (abs(v_iter(i,j,k)-mdv).gt.eps)) then
139
                  v_iter(i,j,k)=v_iter(i,j,k)-alpha*v_anom(i,j,k)
140
               else
141
                  v_iter(i,j,k)=mdv
142
               endif
143
 
144
c              Update temperature
145
               if ((abs(t_anom(i,j,k)-mdv).gt.eps).and.
146
     >             (abs(t_iter(i,j,k)-mdv).gt.eps)) then
147
                  t_iter(i,j,k)=t_iter(i,j,k)-alpha*t_anom(i,j,k)
148
               else
149
                  t_iter(i,j,k)=mdv
150
               endif
151
 
152
c              Update pressure
153
               if ((abs(p_anom(i,j,k)-mdv).gt.eps).and.
154
     >             (abs(p_iter(i,j,k)-mdv).gt.eps)) then
155
                  p_iter(i,j,k)=p_iter(i,j,k)-alpha*p_anom(i,j,k)
156
               else
157
                  p_iter(i,j,k)=mdv
158
               endif
159
 
160
            enddo
161
         enddo
162
      enddo
163
 
164
 
165
c     --------------------------------------------------------------------------------
166
c     Write output
167
c     --------------------------------------------------------------------------------
168
 
169
      varname='U'
170
      call write_inp (u_iter,varname,iterafile,nx,ny,nz)
171
      varname='V'
172
      call write_inp (v_iter,varname,iterafile,nx,ny,nz)
173
      varname='T'
174
      call write_inp (t_iter,varname,iterafile,nx,ny,nz)
175
      varname='P'
176
      call write_inp (p_iter,varname,iterafile,nx,ny,nz)
177
 
178
 
179
      end
180
 
181
 
182
 
183
 
184
c     ********************************************************************************
185
c     * NETCDF INPUT AND OUTPUT                                                      *
186
c     ********************************************************************************
187
 
188
c     --------------------------------------------------------------------------------
189
c     Write input field to netcdf
190
c     --------------------------------------------------------------------------------
191
 
192
      SUBROUTINE write_inp (field,fieldname,pvsrcfile,nx,ny,nz)
193
 
194
c     Read <fieldname> from netcdf file <pvsrcfile> into <field>. The grid is specified 
195
c     by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed whether the input 
196
c     files are consitent with this grid. 
197
 
198
      implicit   none
199
 
200
c     Declaration of subroutine parameters
201
      integer              nx,ny,nz
202
      real                 field (0:nx,0:ny,0:nz)
203
      character*80         fieldname
204
      character*80         pvsrcfile
205
 
206
c     Auxiliary variables
207
      integer              cdfid,stat
208
      integer              vardim(4)
209
      real                 misdat
210
      real                 varmin(4),varmax(4),stag(4)
211
      integer              ndimin,outid,i,j,k
212
      real                 max_th
213
      real                 tmp(0:nx,0:ny,0:nz)
214
      integer              ntimes
215
      real                 time(200)       
216
      integer              nvars
217
      character*80         vnam(100),varname
218
      integer              isok
219
 
220
c     Get grid parameters
221
      call cdfopn(pvsrcfile,cdfid,stat)
222
      if (stat.ne.0) goto 998
223
      call getvars(cdfid,nvars,vnam,stat)
224
      if (stat.ne.0) goto 998
225
      isok=0
226
      varname='TH'
227
      call check_varok(isok,varname,vnam,nvars)
228
      if (isok.eq.0) goto 998
229
      call getdef(cdfid,varname,ndimin,misdat,vardim,
230
     >            varmin,varmax,stag,stat)    
231
      if (stat.ne.0) goto 998
232
      time(1)=0.
233
      call gettimes(cdfid,time,ntimes,stat)  
234
      if (stat.ne.0) goto 998
235
      call clscdf(cdfid,stat)
236
      if (stat.ne.0) goto 998
237
 
238
c     Save variables (write definition, if necessary)
239
      call cdfwopn(pvsrcfile,cdfid,stat)
240
      if (stat.ne.0) goto 998
241
      isok=0
242
      varname=fieldname
243
      call check_varok(isok,varname,vnam,nvars)
244
      if (isok.eq.0) then
245
         call putdef(cdfid,varname,ndimin,misdat,vardim,
246
     >               varmin,varmax,stag,stat)
247
         if (stat.ne.0) goto 998
248
      endif
249
      call putdat(cdfid,varname,time(1),0,field,stat)
250
      print*,'W ',trim(varname),' ',trim(pvsrcfile)
251
      if (stat.ne.0) goto 998
252
 
253
c     Close input netcdf file
254
      call clscdf(cdfid,stat)
255
      if (stat.ne.0) goto 998
256
 
257
      return
258
 
259
c     Exception handling
260
 998  print*,'Write_Inp: Problem with input netcdf file... Stop'
261
      stop
262
 
263
      end
264
 
265
 
266
c     --------------------------------------------------------------------------------
267
c     Read input fields for reference profile
268
c     --------------------------------------------------------------------------------
269
 
270
      SUBROUTINE read_inp (field,fieldname,pvsrcfile,
271
     >                     nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
272
 
273
c     Read <fieldname> from netcdf file <pvsrcfile> into <field>. The grid is specified 
274
c     by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed whether the input 
275
c     files are consitent with this grid. The missing data value is set to <mdv>.
276
 
277
      implicit   none
278
 
279
c     Declaration of subroutine parameters
280
      integer              nx,ny,nz
281
      real                 field(0:nx,0:ny,0:nz)
282
      character*80         fieldname
283
      character*80         pvsrcfile
284
      real                 dx,dy,dz,xmin,ymin,zmin
285
      real                 mdv
286
 
287
c     Numerical and physical parameters
288
      real                 eps
289
      parameter            (eps=0.01)
290
 
291
c     Auxiliary variables
292
      integer              cdfid,stat,cdfid99
293
      integer              vardim(4)
294
      real                 misdat
295
      real                 varmin(4),varmax(4),stag(4)
296
      integer              ndimin,outid,i,j,k
297
      real                 max_th
298
      real                 tmp(nx,ny,nz)
299
      integer              ntimes
300
      real                 time(200)       
301
      integer              nvars
302
      character*80         vnam(100),varname
303
      integer              isok
304
 
305
c     Open the input netcdf file
306
      call cdfopn(pvsrcfile,cdfid,stat)
307
      if (stat.ne.0) goto 998
308
 
309
c     Check whether needed variables are on file
310
      call getvars(cdfid,nvars,vnam,stat)
311
      if (stat.ne.0) goto 998
312
      isok=0
313
      varname=trim(fieldname)
314
      call check_varok(isok,varname,vnam,nvars)
315
      if (isok.eq.0) goto 998
316
 
317
c     Get the grid parameters from theta     
318
      call getdef(cdfid,varname,ndimin,misdat,vardim,
319
     >            varmin,varmax,stag,stat)    
320
      if (stat.ne.0) goto 998
321
      time(1)=0.
322
      call gettimes(cdfid,time,ntimes,stat)  
323
      if (stat.ne.0) goto 998
324
 
325
c     Check whether grid parameters are consistent
326
      if ( (vardim(1).ne.(nx+1)).or.
327
     >     (vardim(2).ne.(ny+1)).or.
328
     >     (vardim(3).ne.(nz+1)).or.
329
     >     (abs(varmin(1)-xmin).gt.eps).or.
330
     >     (abs(varmin(2)-ymin).gt.eps).or.
331
     >     (abs(varmin(3)-zmin).gt.eps).or.
332
     >     (abs((varmax(1)-varmin(1))/real(vardim(1)-1)-dx).gt.eps).or.
333
     >     (abs((varmax(2)-varmin(2))/real(vardim(2)-1)-dy).gt.eps).or.
334
     >     (abs((varmax(3)-varmin(3))/real(vardim(3)-1)-dz).gt.eps) ) 
335
     >then
336
         print*,'Input grid inconsitency...'
337
         print*,'  Nx      = ',vardim(1),nx+1
338
         print*,'  Ny      = ',vardim(2),ny+1
339
         print*,'  Nz      = ',vardim(3),nz+1
340
         print*,'  Varminx = ',varmin(1),xmin
341
         print*,'  Varminy = ',varmin(2),ymin
342
         print*,'  Varminz = ',varmin(3),zmin
343
         print*,'  Dx      = ',(varmax(1)-varmin(1))/real(nx-1),dx
344
         print*,'  Dy      = ',(varmax(2)-varmin(2))/real(ny-1),dy
345
         print*,'  Dz      = ',(varmax(3)-varmin(3))/real(nz-1),dz
346
         goto 998
347
      endif
348
 
349
c     Load variables
350
      call getdef(cdfid,varname,ndimin,misdat,vardim,
351
     >            varmin,varmax,stag,stat)
352
      if (stat.ne.0) goto 998
353
      call getdat(cdfid,varname,time(1),0,field,stat)
354
      print*, 'R ',trim(varname),' ',trim(pvsrcfile)
355
      if (stat.ne.0) goto 998
356
 
357
c     Close input netcdf file
358
      call clscdf(cdfid,stat)
359
      if (stat.ne.0) goto 998
360
 
361
c     Set missing data value to <mdv>
362
      do i=1,nx
363
         do j=1,ny
364
            do k=1,nz
365
               if (abs(field(i,j,k)-misdat).lt.eps) then
366
                  field(i,j,k)=mdv
367
               endif
368
            enddo
369
         enddo
370
      enddo
371
 
372
      return
373
 
374
c     Exception handling
375
 998  print*,'Read_Inp: Problem with input netcdf file... Stop'
376
      stop
377
 
378
      end
379
 
380
c     --------------------------------------------------------------------------------
381
c     Check whether variable is found on netcdf file
382
c     --------------------------------------------------------------------------------
383
 
384
      subroutine check_varok (isok,varname,varlist,nvars)
385
 
386
c     Check whether the variable <varname> is in the list <varlist(nvars)>. If this is 
387
C     the case, <isok> is incremented by 1. Otherwise <isok> keeps its value.
388
 
389
      implicit none
390
 
391
c     Declaraion of subroutine parameters
392
      integer      isok
393
      integer      nvars
394
      character*80 varname
395
      character*80 varlist(nvars)
396
 
397
c     Auxiliary variables
398
      integer      i
399
 
400
c     Main
401
      do i=1,nvars
402
         if (trim(varname).eq.trim(varlist(i))) isok=isok+1
403
      enddo
404
 
405
      end
406
 
407
c     --------------------------------------------------------------------------------
408
c     Get grid parameters
409
c     --------------------------------------------------------------------------------
410
 
411
      subroutine read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv,
412
     >                     pvsrcfile)
413
 
414
c     Get the grid parameters from the variable <THETA> on the input file <pvsrcfile>.
415
c     The grid parameters are
416
c        nx,ny,nz                : Number of grid points in x, y and z direction
417
c        xmin,ymin,zmin          : Minimum domain coordinates in x, y and z direction
418
c        xmax,ymax,zmax          : Maximal domain coordinates in x, y and z direction
419
c        dx,dy,dz                : Horizontal and vertical resolution
420
c     Additionally, it is checked whether the vertical grid is equally spaced. If ok,
421
c     the grid paramters are transformed from lon/lat to distance (in meters)
422
 
423
      implicit none
424
 
425
c     Declaration of subroutine parameters
426
      character*80   pvsrcfile
427
      integer        nx,ny,nz
428
      real           dx,dy,dz
429
      real           xmin,ymin,zmin,xmax,ymax,zmax
430
      real           mdv
431
 
432
c     Numerical epsilon and other physical/geoemtrical parameters
433
      real           eps
434
      parameter      (eps=0.01)
435
 
436
c     Auxiliary variables
437
      integer        cdfid,cstid
438
      integer        ierr
439
      character*80   vnam(100),varname
440
      integer        nvars
441
      integer        isok
442
      integer        vardim(4)
443
      real           misdat
444
      real           varmin(4),varmax(4),stag(4)
445
      real           aklev(1000),bklev(1000),aklay(1000),bklay(1000)
446
      real           dh
447
      character*80   csn
448
      integer        ndim
449
      integer        i
450
 
451
c     Get all grid parameters
452
      call cdfopn(pvsrcfile,cdfid,ierr)
453
      if (ierr.ne.0) goto 998
454
      call getvars(cdfid,nvars,vnam,ierr)
455
      if (ierr.ne.0) goto 998
456
      isok=0
457
      varname='TH'
458
      call check_varok(isok,varname,vnam,nvars)
459
      if (isok.eq.0) goto 998
460
      call getcfn(cdfid,csn,ierr)
461
      if (ierr.ne.0) goto 998
462
      call cdfopn(csn,cstid,ierr)
463
      if (ierr.ne.0) goto 998
464
      call getdef(cdfid,varname,ndim,misdat,vardim,varmin,varmax,
465
     >            stag,ierr)
466
      if (ierr.ne.0) goto 998
467
      nx=vardim(1)
468
      ny=vardim(2)
469
      nz=vardim(3)
470
      xmin=varmin(1)
471
      ymin=varmin(2)
472
      zmin=varmin(3)
473
      call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)
474
      if (ierr.ne.0) goto 998
475
      call getgrid(cstid,dx,dy,ierr)
476
      if (ierr.ne.0) goto 998
477
      xmax=varmax(1)
478
      ymax=varmax(2)
479
      zmax=varmax(3)
480
      dz=(zmax-zmin)/real(nz-1)
481
      call clscdf(cstid,ierr)
482
      if (ierr.ne.0) goto 998
483
      call clscdf(cdfid,ierr)
484
      if (ierr.ne.0) goto 998
485
 
486
c     Check whether the grid is equallay spaced in the vertical
487
      do i=1,nz-1
488
         dh=aklev(i+1)-aklev(i)
489
         if (abs(dh-dz).gt.eps) then
490
            print*,'Aklev: Vertical grid must be equally spaced... Stop'
491
            print*,(aklev(i),i=1,nz)
492
            stop
493
         endif
494
         dh=aklay(i+1)-aklay(i)
495
         if (abs(dh-dz).gt.eps) then
496
            print*,'Aklay: Vertical grid must be equally spaced... Stop'
497
            print*,(aklay(i),i=1,nz)
498
            stop
499
         endif
500
      enddo
501
 
502
c     Set missing data value
503
      mdv=misdat
504
 
505
      return
506
 
507
c     Exception handling
508
 998  print*,'Read_Dim: Problem with input netcdf file... Stop'
509
      stop
510
 
511
      end
512
 
513