Subversion Repositories pvinversion.ecmwf

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
      PROGRAM difference
2
 
3
c     **************************************************************************
4
c     * Get the difference between two runs                                    *
5
c     * Michael Sprenger / Autumn 2006                                         *
6
c     **************************************************************************
7
 
8
      implicit none
9
 
10
c     --------------------------------------------------------------------------
11
c     Declaration of subroutine parameters
12
c     --------------------------------------------------------------------------
13
 
14
c     Physical and numerical parameters
15
      real      eps
16
      parameter (eps=0.01)
17
      integer   nmax
18
      parameter (nmax=300*300*200)
19
 
20
c     Variables for input file 1
21
      character*80 i1_filename
22
      real         i1_varmin(4),i1_varmax(4),i1_stag(4)
23
      integer      i1_vardim(4)
24
      real         i1_mdv
25
      integer      i1_ndim
26
      integer      i1_nx,i1_ny,i1_nz
27
      real         i1_time
28
      integer      i1_nvars
29
      character*80 i1_vnam(100)
30
      real         i1_field(nmax)
31
 
32
c     Variables for input file 2
33
      character*80 i2_filename
34
      real         i2_varmin(4),i2_varmax(4),i2_stag(4)
35
      integer      i2_vardim(4)
36
      real         i2_mdv
37
      integer      i2_ndim
38
      integer      i2_nx,i2_ny,i2_nz
39
      real         i2_time
40
      integer      i2_nvars
41
      character*80 i2_vnam(100)
42
      real         i2_field(nmax)
43
 
44
c     Variables for output file
45
      character*80 o_filename
46
      real         o_varmin(4),o_varmax(4),o_stag(4)
47
      integer      o_vardim(4)
48
      real         o_mdv
49
      integer      o_ndim
50
      integer      o_nx,o_ny,o_nz
51
      real         o_time
52
      real         o_field(nmax)
53
 
54
c     Auxiliary variables
55
      integer      i,j,k
56
      integer      cdfid,ierr
57
      integer      isok
58
      character*80 cfn,varname
59
 
60
c     --------------------------------------------------------------------------------
61
c     Input
62
c     --------------------------------------------------------------------------------
63
 
64
      print*,'********************************************************'
65
      print*,'* DIFFERENCE                                           *'
66
      print*,'********************************************************'
67
 
68
c     Read parameter file
69
      open(10,file='fort.10')
70
       read(10,*) i1_filename
71
       read(10,*) i2_filename
72
       read(10,*) o_filename
73
      close(10)
74
      print*
75
      print*,trim(i1_filename)
76
      print*,trim(i2_filename)
77
      print*,trim(o_filename)
78
      print*
79
 
80
c     Get a list of all variables on the two input files
81
      call cdfopn(i1_filename,cdfid,ierr)
82
      if (ierr.ne.0) goto 998
83
      call getvars(cdfid,i1_nvars,i1_vnam,ierr)
84
      if (ierr.ne.0) goto 998
85
      call getdef(cdfid,i1_vnam(i1_nvars),i1_ndim,i1_mdv,i1_vardim,
86
     >            i1_varmin,i1_varmax,i1_stag,ierr)
87
      call clscdf(cdfid,ierr)
88
      print*,(trim(i1_vnam(i))//'  ',i=1,i1_nvars)
89
 
90
      call cdfopn(i2_filename,cdfid,ierr)
91
      if (ierr.ne.0) goto 997
92
      call getvars(cdfid,i2_nvars,i2_vnam,ierr)
93
      if (ierr.ne.0) goto 997
94
      call clscdf(cdfid,ierr)
95
      print*,(trim(i2_vnam(i))//'  ',i=1,i2_nvars)
96
      print*
97
 
98
c     Create the output file
99
      o_ndim=i1_ndim
100
      do i=1,i1_ndim
101
         o_varmin(i)=i1_varmin(i)
102
         o_varmax(i)=i1_varmax(i)
103
      enddo
104
      cfn=trim(o_filename)//'_cst'
105
      call crecdf(trim(o_filename),cdfid,o_varmin,o_varmax,
106
     >            o_ndim,trim(cfn),ierr)
107
      if (ierr.ne.0) goto 996
108
      call clscdf(cdfid,ierr)
109
 
110
c     --------------------------------------------------------------------------------
111
c     Loop through all variables
112
c     --------------------------------------------------------------------------------
113
 
114
      do i=1,i1_nvars
115
 
116
c        Check wether the variable is available on both files
117
         varname=i1_vnam(i)
118
         if (varname.eq.'time') goto 100
119
         isok=0
120
         call check_varok (isok,varname,i2_vnam,i2_nvars)
121
         if (isok.eq.0) goto 100
122
 
123
c        Read first input file
124
         call cdfopn(i1_filename,cdfid,ierr)
125
         if (ierr.ne.0) goto 998
126
         call getdef(cdfid,varname,i1_ndim,i1_mdv,i1_vardim,
127
     >               i1_varmin,i1_varmax,i1_stag,ierr)
128
         if (ierr.ne.0) goto 998
129
         call getdat(cdfid,varname,0.,0,i1_field,ierr)
130
         if (ierr.ne.0) goto 998
131
         call clscdf(cdfid,ierr)
132
 
133
c        Read second input file
134
         call cdfopn(i2_filename,cdfid,ierr)
135
         if (ierr.ne.0) goto 998
136
         call getdef(cdfid,varname,i2_ndim,i2_mdv,i2_vardim,
137
     >               i2_varmin,i2_varmax,i2_stag,ierr)
138
         if (ierr.ne.0) goto 998
139
         call getdat(cdfid,varname,0.,0,i2_field,ierr)
140
         if (ierr.ne.0) goto 998
141
         call clscdf(cdfid,ierr)
142
 
143
c        Consistency check
144
         if (i1_ndim.ne.i2_ndim) then
145
            print*,'Inconsistent input files... Stop',i1_ndim,i2_ndim
146
            stop
147
         endif
148
         do j=1,3
149
            if ( (i1_vardim(j).ne.i2_vardim(j)).or.
150
     >           (abs(i1_varmin(j)-i2_varmin(j)).gt.eps).or.
151
     >           (abs(i1_varmax(j)-i2_varmax(j)).gt.eps).or.
152
     >           (abs(i1_stag(j)-i2_stag(j)).gt.eps)) then
153
               print*,'Inconsistent input files... Stop'
154
               print*,j,i1_varmin(j),i2_varmin(j)
155
               print*,j,i1_varmax(j),i2_varmax(j)
156
               print*,j,i1_vardim(j),i2_vardim(j)
157
               print*,j,i1_stag(j),  i2_stag(j)
158
               stop
159
            endif
160
         enddo
161
 
162
c        Get the difference
163
         do j=1,i1_vardim(1)*i1_vardim(2)*i1_vardim(3)
164
            if ( (abs(i1_field(j)-i1_mdv).gt.eps).and.
165
     >           (abs(i2_field(j)-i2_mdv).gt.eps) ) then
166
               o_field(j)=i1_field(j)-i2_field(j)
167
            else
168
               o_field(j)=i1_mdv
169
            endif
170
         enddo
171
 
172
c        Write to output file
173
         o_ndim=i1_ndim
174
         o_mdv =i1_mdv
175
         do j=1,i1_ndim
176
            o_vardim(j)=i1_vardim(j)
177
            o_varmin(j)=i1_varmin(j)
178
            o_varmax(j)=i1_varmax(j)
179
            o_stag(j)  =i1_stag(j)
180
         enddo       
181
         call cdfwopn(o_filename,cdfid,ierr)
182
         if (ierr.ne.0) goto 996
183
         call putdef(cdfid,varname,o_ndim,o_mdv,o_vardim,
184
     >               o_varmin,o_varmax,o_stag,ierr)
185
         if (ierr.ne.0) goto 996
186
         call putdat(cdfid,varname,0.,0,o_field,ierr)
187
         if (ierr.ne.0) goto 996
188
         print*,'W ',trim(varname),' ',trim(o_filename)
189
         call clscdf(cdfid,ierr)
190
         if (ierr.ne.0) goto 996
191
 
192
c        Next 
193
 100     continue
194
 
195
      enddo
196
 
197
c     -----------------------------------------------------------------
198
c     Exception handling and format specs
199
c     -----------------------------------------------------------------
200
 
201
      stop
202
 
203
 998  print*,'Problems with input file 1  ',trim(i1_filename)
204
      stop
205
 
206
 997  print*,'Problems with input file 2  ',trim(i2_filename)
207
      stop
208
 
209
 996  print*,'Problems with output file   ',trim(o_filename)
210
      stop
211
 
212
      end
213
 
214
 
215
c     --------------------------------------------------------------------------------
216
c     Check whether variable is found on netcdf file
217
c     --------------------------------------------------------------------------------
218
 
219
      subroutine check_varok (isok,varname,varlist,nvars)
220
 
221
c     Check whether the variable <varname> is in the list <varlist(nvars)>. If this is 
222
C     the case, <isok> is incremented by 1. Otherwise <isok> keeps its value.
223
 
224
      implicit none
225
 
226
c     Declaraion of subroutine parameters
227
      integer      isok
228
      integer      nvars
229
      character*80 varname
230
      character*80 varlist(nvars)
231
 
232
c     Auxiliary variables
233
      integer      i
234
 
235
c     Main
236
      do i=1,nvars
237
         if (trim(varname).eq.trim(varlist(i))) isok=isok+1
238
      enddo
239
 
240
      end
241