| 3 |
michaesp |
1 |
PROGRAM p2z
|
|
|
2 |
|
|
|
3 |
c Calculate the geopotential and add it to the P file
|
|
|
4 |
c Michael Sprenger / Spring 2006
|
|
|
5 |
|
|
|
6 |
implicit none
|
|
|
7 |
|
|
|
8 |
c ---------------------------------------------------------------
|
|
|
9 |
c Declaration of variables
|
|
|
10 |
c ---------------------------------------------------------------
|
|
|
11 |
|
|
|
12 |
c Variables for input P file : model level
|
|
|
13 |
character*80 ml_pfn
|
|
|
14 |
real ml_varmin(4),ml_varmax(4),ml_stag(4)
|
|
|
15 |
integer ml_vardim(4)
|
|
|
16 |
real ml_mdv
|
|
|
17 |
integer ml_ndim
|
|
|
18 |
integer ml_nx,ml_ny,ml_nz
|
|
|
19 |
real ml_xmin,ml_xmax,ml_ymin,ml_ymax,ml_dx,ml_dy
|
|
|
20 |
integer ml_ntimes
|
|
|
21 |
real ml_aklev(500),ml_bklev(500)
|
|
|
22 |
real ml_aklay(500),ml_bklay(500)
|
|
|
23 |
real ml_time
|
|
|
24 |
real ml_pollon,ml_pollat
|
|
|
25 |
integer ml_nvars
|
|
|
26 |
character*80 ml_vnam(100)
|
|
|
27 |
integer ml_idate(5)
|
|
|
28 |
real,allocatable, dimension (:,:) :: ml_ps,ml_zb
|
|
|
29 |
real,allocatable, dimension (:,:,:) :: ml_t3,ml_q3,ml_p3,ml_tv3
|
|
|
30 |
real,allocatable, dimension (:,:,:) :: ml_z3
|
|
|
31 |
|
|
|
32 |
c Variables for input Z file : pressure level
|
|
|
33 |
character*80 pl_zfn
|
|
|
34 |
real pl_varmin(4),pl_varmax(4),pl_stag(4)
|
|
|
35 |
integer pl_vardim(4)
|
|
|
36 |
real pl_mdv
|
|
|
37 |
integer pl_ndim
|
|
|
38 |
integer pl_nx,pl_ny,pl_nz
|
|
|
39 |
real pl_xmin,pl_xmax,pl_ymin,pl_ymax,pl_dx,pl_dy
|
|
|
40 |
integer pl_ntimes
|
|
|
41 |
real pl_aklev(500),pl_bklev(500)
|
|
|
42 |
real pl_aklay(500),pl_bklay(500)
|
|
|
43 |
real pl_time
|
|
|
44 |
real pl_pollon,pl_pollat
|
|
|
45 |
integer pl_nvars
|
|
|
46 |
character*80 pl_vnam(100)
|
|
|
47 |
integer pl_idate(5)
|
|
|
48 |
real,allocatable, dimension (:,:,:) :: pl_z3,pl_p3
|
|
|
49 |
|
|
|
50 |
c Variables for output Z file : height level
|
|
|
51 |
character*80 zl_ofn
|
|
|
52 |
real zl_varmin(4),zl_varmax(4),zl_stag(4)
|
|
|
53 |
integer zl_vardim(4)
|
|
|
54 |
real zl_mdv
|
|
|
55 |
integer zl_ndim
|
|
|
56 |
integer zl_nx,zl_ny,zl_nz
|
|
|
57 |
real zl_xmin,zl_xmax,zl_ymin,zl_ymax,zl_dx,zl_dy
|
|
|
58 |
integer zl_ntimes
|
|
|
59 |
real zl_aklev(500),zl_bklev(500)
|
|
|
60 |
real zl_aklay(500),zl_bklay(500)
|
|
|
61 |
real zl_time
|
|
|
62 |
real zl_pollon,zl_pollat
|
|
|
63 |
integer zl_idate(5)
|
|
|
64 |
real,allocatable, dimension (:,:,:) :: zl_field
|
|
|
65 |
real,allocatable, dimension (:,:,:) :: zl_p
|
|
|
66 |
|
|
|
67 |
|
|
|
68 |
c Variables for input P,S file : model level
|
|
|
69 |
character*80 in_ofn
|
|
|
70 |
real in_varmin(4),in_varmax(4),in_stag(4)
|
|
|
71 |
integer in_vardim(4)
|
|
|
72 |
real in_mdv
|
|
|
73 |
integer in_ndim
|
|
|
74 |
integer in_nx,in_ny,in_nz
|
|
|
75 |
real in_xmin,in_xmax,in_ymin,in_ymax,in_dx,in_dy
|
|
|
76 |
integer in_ntimes
|
|
|
77 |
real in_aklev(500),in_bklev(500)
|
|
|
78 |
real in_aklay(500),in_bklay(500)
|
|
|
79 |
real in_time
|
|
|
80 |
real in_pollon,in_pollat
|
|
|
81 |
integer in_nvars
|
|
|
82 |
character*80 in_vnam(100)
|
|
|
83 |
integer in_idate(5)
|
|
|
84 |
real,allocatable, dimension (:,:,:) :: in_field
|
|
|
85 |
|
|
|
86 |
c Physical and numerical parameters
|
|
|
87 |
real g
|
|
|
88 |
parameter (g=9.80665)
|
|
|
89 |
real eps
|
|
|
90 |
parameter (eps=0.01)
|
|
|
91 |
real tzero
|
|
|
92 |
parameter (tzero=273.15)
|
|
|
93 |
real kappa
|
|
|
94 |
parameter (kappa=0.6078)
|
|
|
95 |
real zerodiv
|
|
|
96 |
parameter (zerodiv=0.0000000001)
|
|
|
97 |
real dpmin
|
|
|
98 |
parameter (dpmin=10.)
|
|
|
99 |
real rdg
|
|
|
100 |
parameter (rdg=29.271)
|
|
|
101 |
|
|
|
102 |
c Flag for test mode
|
|
|
103 |
integer test
|
|
|
104 |
parameter (test=0)
|
|
|
105 |
character*80 testfile
|
|
|
106 |
parameter (testfile='TEST')
|
|
|
107 |
|
|
|
108 |
c Variables and levels for interpolation onto z levels
|
|
|
109 |
character*80 levelfile,varfile
|
|
|
110 |
integer nvar,nlev
|
|
|
111 |
character*80 varinp(100),varout(100),varsrc(100)
|
|
|
112 |
real zlev(500)
|
|
|
113 |
|
|
|
114 |
c Auxiliary variables
|
|
|
115 |
integer ierr
|
|
|
116 |
integer cdfid,cstid
|
|
|
117 |
character*80 cfn
|
|
|
118 |
integer stat
|
|
|
119 |
real time
|
|
|
120 |
real tv1(1000),z1(1000),p1(1000),f1(1000)
|
|
|
121 |
real spline_tv1(1000),spline_f1(1000),spline_z1(1000)
|
|
|
122 |
real pu,po,zu,zo,p,z,dp,p0,tvu,tvo,ff
|
|
|
123 |
integer i,j,k,l
|
|
|
124 |
integer lmin,n
|
|
|
125 |
character*80 varname,cdfname
|
|
|
126 |
integer idate(5),stdate(5),datar(14)
|
|
|
127 |
integer isok
|
|
|
128 |
character*80 name
|
|
|
129 |
real zmin,dz
|
|
|
130 |
|
|
|
131 |
c -----------------------------------------------------------------
|
|
|
132 |
c Read input fields
|
|
|
133 |
c -----------------------------------------------------------------
|
|
|
134 |
|
|
|
135 |
print*,'*********************************************************'
|
|
|
136 |
print*,'* p2z *'
|
|
|
137 |
print*,'*********************************************************'
|
|
|
138 |
|
|
|
139 |
c Read in the parameter file
|
|
|
140 |
open(10,file='fort.10')
|
|
|
141 |
|
|
|
142 |
read(10,*) ml_pfn
|
|
|
143 |
read(10,*) pl_zfn
|
|
|
144 |
read(10,*) zl_ofn
|
|
|
145 |
|
|
|
146 |
read(10,*) name,zmin
|
|
|
147 |
if ( trim(name).ne.'GEO_ZMIN') stop
|
|
|
148 |
read(10,*) name,nlev
|
|
|
149 |
if ( trim(name).ne.'GEO_NZ' ) stop
|
|
|
150 |
read(10,*) name,dz
|
|
|
151 |
if ( trim(name).ne.'GEO_DZ' ) stop
|
|
|
152 |
do i=1,nlev
|
|
|
153 |
zlev(i)=zmin+real(i-1)*dz
|
|
|
154 |
enddo
|
|
|
155 |
|
|
|
156 |
nvar=1
|
|
|
157 |
102 continue
|
|
|
158 |
read(10,*,end=103) varinp(nvar),varout(nvar),varsrc(nvar)
|
|
|
159 |
nvar=nvar+1
|
|
|
160 |
goto 102
|
|
|
161 |
103 continue
|
|
|
162 |
nvar=nvar-1
|
|
|
163 |
print*
|
|
|
164 |
do i=1,nvar
|
|
|
165 |
write(*,'(a10,a10,a30)') trim(varinp(i)),trim(varout(i)),
|
|
|
166 |
> trim(varsrc(i))
|
|
|
167 |
enddo
|
|
|
168 |
print*
|
|
|
169 |
|
|
|
170 |
close(10)
|
|
|
171 |
|
|
|
172 |
|
|
|
173 |
c Get grid description for P file : model level
|
|
|
174 |
call cdfopn(ml_pfn,cdfid,ierr)
|
|
|
175 |
if (ierr.ne.0) goto 998
|
|
|
176 |
call getcfn(cdfid,cfn,ierr)
|
|
|
177 |
if (ierr.ne.0) goto 998
|
|
|
178 |
call cdfopn(cfn,cstid,ierr)
|
|
|
179 |
if (ierr.ne.0) goto 998
|
|
|
180 |
call getvars(cdfid,ml_nvars,ml_vnam,ierr)
|
|
|
181 |
varname='T'
|
|
|
182 |
isok=0
|
|
|
183 |
call check_varok (isok,varname,ml_vnam,ml_nvars)
|
|
|
184 |
if (isok.ne.1) goto 998
|
|
|
185 |
call getdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
|
|
|
186 |
> ml_varmin,ml_varmax,ml_stag,ierr)
|
|
|
187 |
if (ierr.ne.0) goto 998
|
|
|
188 |
ml_nx =ml_vardim(1)
|
|
|
189 |
ml_ny =ml_vardim(2)
|
|
|
190 |
ml_nz =ml_vardim(3)
|
|
|
191 |
ml_xmin=ml_varmin(1)
|
|
|
192 |
ml_ymin=ml_varmin(2)
|
|
|
193 |
call getlevs(cstid,ml_nz,ml_aklev,ml_bklev,ml_aklay,ml_bklay,ierr)
|
|
|
194 |
call getgrid(cstid,ml_dx,ml_dy,ierr)
|
|
|
195 |
ml_xmax=ml_xmin+real(ml_nx-1)*ml_dx
|
|
|
196 |
ml_ymax=ml_ymin+real(ml_ny-1)*ml_dy
|
|
|
197 |
call gettimes(cdfid,ml_time,ml_ntimes,ierr)
|
|
|
198 |
call getstart(cstid,ml_idate,ierr)
|
|
|
199 |
call getpole(cstid,ml_pollon,ml_pollat,ierr)
|
|
|
200 |
call clscdf(cstid,ierr)
|
|
|
201 |
call clscdf(cdfid,ierr)
|
|
|
202 |
|
|
|
203 |
c Get grid description for Z file : pressure level
|
|
|
204 |
call cdfopn(pl_zfn,cdfid,ierr)
|
|
|
205 |
if (ierr.ne.0) goto 998
|
|
|
206 |
call getcfn(cdfid,cfn,ierr)
|
|
|
207 |
if (ierr.ne.0) goto 998
|
|
|
208 |
call cdfopn(cfn,cstid,ierr)
|
|
|
209 |
if (ierr.ne.0) goto 998
|
|
|
210 |
call getvars(cdfid,pl_nvars,pl_vnam,ierr)
|
|
|
211 |
varname='Z'
|
|
|
212 |
isok=0
|
|
|
213 |
call check_varok (isok,varname,pl_vnam,pl_nvars)
|
|
|
214 |
if (isok.ne.1) goto 998
|
|
|
215 |
call getdef(cdfid,varname,pl_ndim,pl_mdv,pl_vardim,
|
|
|
216 |
> pl_varmin,pl_varmax,pl_stag,ierr)
|
|
|
217 |
if (ierr.ne.0) goto 998
|
|
|
218 |
pl_nx =pl_vardim(1)
|
|
|
219 |
pl_ny =pl_vardim(2)
|
|
|
220 |
pl_nz =pl_vardim(3)
|
|
|
221 |
pl_xmin=pl_varmin(1)
|
|
|
222 |
pl_ymin=pl_varmin(2)
|
|
|
223 |
call getlevs(cstid,pl_nz,pl_aklev,pl_bklev,pl_aklay,pl_bklay,ierr)
|
|
|
224 |
call getgrid(cstid,pl_dx,pl_dy,ierr)
|
|
|
225 |
pl_xmax=pl_xmin+real(pl_nx-1)*pl_dx
|
|
|
226 |
pl_ymax=pl_ymin+real(pl_ny-1)*pl_dy
|
|
|
227 |
call gettimes(cdfid,pl_time,pl_ntimes,ierr)
|
|
|
228 |
call getstart(cstid,pl_idate,ierr)
|
|
|
229 |
call getpole(cstid,pl_pollon,pl_pollat,ierr)
|
|
|
230 |
call clscdf(cstid,ierr)
|
|
|
231 |
call clscdf(cdfid,ierr)
|
|
|
232 |
|
|
|
233 |
c Set grid description for output file : height level
|
|
|
234 |
zl_vardim(1) = pl_vardim(1)
|
|
|
235 |
zl_vardim(2) = pl_vardim(2)
|
|
|
236 |
zl_vardim(3) = nlev
|
|
|
237 |
zl_varmin(1) = pl_varmin(1)
|
|
|
238 |
zl_varmin(2) = pl_varmin(2)
|
|
|
239 |
zl_varmin(3) = zlev(1)
|
|
|
240 |
zl_varmax(1) = pl_varmax(1)
|
|
|
241 |
zl_varmax(2) = pl_varmax(2)
|
|
|
242 |
zl_varmax(3) = zlev(nlev)
|
|
|
243 |
do i=1,nlev
|
|
|
244 |
zl_aklay(i) = zlev(i)
|
|
|
245 |
zl_bklay(i) = 0.
|
|
|
246 |
zl_aklev(i) = zlev(i)
|
|
|
247 |
zl_bklev(i) = 0.
|
|
|
248 |
enddo
|
|
|
249 |
zl_dx = pl_dx
|
|
|
250 |
zl_dy = pl_dy
|
|
|
251 |
zl_time = pl_time
|
|
|
252 |
zl_ntimes = pl_ntimes
|
|
|
253 |
zl_ndim = 4
|
|
|
254 |
zl_mdv = pl_mdv
|
|
|
255 |
zl_nx = zl_vardim(1)
|
|
|
256 |
zl_ny = zl_vardim(2)
|
|
|
257 |
zl_nz = zl_vardim(3)
|
|
|
258 |
zl_xmin = zl_varmin(1)
|
|
|
259 |
zl_ymin = zl_varmin(2)
|
|
|
260 |
zl_pollon = ml_pollon
|
|
|
261 |
zl_pollat = ml_pollat
|
|
|
262 |
do i=1,5
|
|
|
263 |
zl_idate(i) = ml_idate(i)
|
|
|
264 |
enddo
|
|
|
265 |
|
|
|
266 |
c Consitency check for the grids
|
|
|
267 |
if ( (ml_nx.ne.pl_nx).or.
|
|
|
268 |
> (ml_ny.ne.pl_ny).or.
|
|
|
269 |
> (abs(ml_xmin-pl_xmin ).gt.eps).or.
|
|
|
270 |
> (abs(ml_ymin-pl_ymin ).gt.eps).or.
|
|
|
271 |
> (abs(ml_xmax-pl_xmax ).gt.eps).or.
|
|
|
272 |
> (abs(ml_ymax-pl_ymax ).gt.eps).or.
|
|
|
273 |
> (abs(ml_dx -pl_dx ).gt.eps).or.
|
|
|
274 |
> (abs(ml_dy -pl_dy ).gt.eps).or.
|
|
|
275 |
> (abs(ml_time-pl_time ).gt.eps).or.
|
|
|
276 |
> (abs(ml_pollon-pl_pollon).gt.eps).or.
|
|
|
277 |
> (abs(ml_pollat-pl_pollat).gt.eps)) then
|
|
|
278 |
print*,'Input P and Z grids are not consistent... Stop'
|
|
|
279 |
print*
|
|
|
280 |
print*,'Xmin : ',ml_xmin,pl_xmin
|
|
|
281 |
print*,'Ymin : ',ml_ymin,pl_ymin
|
|
|
282 |
print*,'Xmax : ',ml_xmax,pl_xmax
|
|
|
283 |
print*,'Ymax : ',ml_ymax,pl_ymax
|
|
|
284 |
print*,'Dx : ',ml_dx, pl_dx
|
|
|
285 |
print*,'Dy : ',ml_dy, pl_dy
|
|
|
286 |
print*,'Time : ',ml_time,pl_time
|
|
|
287 |
print*,'Pollon : ',ml_pollon,pl_pollon
|
|
|
288 |
print*,'Pollat : ',ml_pollat,pl_pollat
|
|
|
289 |
stop
|
|
|
290 |
endif
|
|
|
291 |
|
|
|
292 |
c Allocate memory for all fields
|
|
|
293 |
allocate(ml_ps(ml_nx,ml_ny),stat=stat)
|
|
|
294 |
if (stat.ne.0) print*,'*** error allocating array ml_ps ***'
|
|
|
295 |
allocate(ml_zb(ml_nx,ml_ny),stat=stat)
|
|
|
296 |
if (stat.ne.0) print*,'*** error allocating array ml_zb ***'
|
|
|
297 |
allocate(ml_p3(ml_nx,ml_ny,ml_nz),stat=stat)
|
|
|
298 |
if (stat.ne.0) print*,'*** error allocating array ml_p3 ***'
|
|
|
299 |
allocate(ml_t3(ml_nx,ml_ny,ml_nz),stat=stat)
|
|
|
300 |
if (stat.ne.0) print*,'*** error allocating array ml_t3 ***'
|
|
|
301 |
allocate(ml_q3(ml_nx,ml_ny,ml_nz),stat=stat)
|
|
|
302 |
if (stat.ne.0) print*,'*** error allocating array ml_q3 ***'
|
|
|
303 |
allocate(ml_tv3(ml_nx,ml_ny,ml_nz),stat=stat)
|
|
|
304 |
if (stat.ne.0) print*,'*** error allocating array ml_tv3 ***'
|
|
|
305 |
allocate(ml_z3(ml_nx,ml_ny,ml_nz),stat=stat)
|
|
|
306 |
if (stat.ne.0) print*,'*** error allocating array ml_z3 ***'
|
|
|
307 |
|
|
|
308 |
allocate(in_field(ml_nx,ml_ny,ml_nz),stat=stat)
|
|
|
309 |
if (stat.ne.0) print*,'*** error allocating array in_field ***'
|
|
|
310 |
|
|
|
311 |
allocate(pl_z3(pl_nx,pl_ny,pl_nz),stat=stat)
|
|
|
312 |
if (stat.ne.0) print*,'*** error allocating array pl_z3 ***'
|
|
|
313 |
allocate(pl_p3(pl_nx,pl_ny,pl_nz),stat=stat)
|
|
|
314 |
if (stat.ne.0) print*,'*** error allocating array pl_p3 ***'
|
|
|
315 |
|
|
|
316 |
allocate(zl_field(zl_nx,zl_ny,zl_nz),stat=stat)
|
|
|
317 |
if (stat.ne.0) print*,'*** error allocating array zl_field ***'
|
|
|
318 |
allocate(zl_p(zl_nx,zl_ny,zl_nz),stat=stat)
|
|
|
319 |
if (stat.ne.0) print*,'*** error allocating array zl_p ***'
|
|
|
320 |
|
|
|
321 |
c Read T, Q, PS from P file
|
|
|
322 |
call cdfopn(ml_pfn,cdfid,ierr)
|
|
|
323 |
if (ierr.ne.0) goto 998
|
|
|
324 |
isok=0
|
|
|
325 |
varname='T'
|
|
|
326 |
call check_varok (isok,varname, ml_vnam,ml_nvars)
|
|
|
327 |
varname='Q'
|
|
|
328 |
call check_varok (isok,varname, ml_vnam,ml_nvars)
|
|
|
329 |
varname='PS'
|
|
|
330 |
call check_varok (isok,varname,ml_vnam,ml_nvars)
|
|
|
331 |
if (isok.ne.3) goto 998
|
|
|
332 |
print*,'R T ',trim(ml_pfn)
|
|
|
333 |
call getdat(cdfid,'T',ml_time,0,ml_t3,ierr)
|
|
|
334 |
print*,'R Q ',trim(ml_pfn)
|
|
|
335 |
call getdat(cdfid,'Q',ml_time,0,ml_q3,ierr)
|
|
|
336 |
print*,'R PS ',trim(ml_pfn)
|
|
|
337 |
call getdat(cdfid,'PS',ml_time,1,ml_ps,ierr)
|
|
|
338 |
call clscdf(cdfid,ierr)
|
|
|
339 |
|
|
|
340 |
c Read Z from Z file
|
|
|
341 |
call cdfopn(pl_zfn,cdfid,ierr)
|
|
|
342 |
if (ierr.ne.0) goto 998
|
|
|
343 |
isok=0
|
|
|
344 |
varname='Z'
|
|
|
345 |
call check_varok (isok,varname,pl_vnam,pl_nvars)
|
|
|
346 |
if (isok.ne.1) goto 998
|
|
|
347 |
print*,'R Z ',trim(pl_zfn)
|
|
|
348 |
call getdat(cdfid,varname,pl_time,0,pl_z3,ierr)
|
|
|
349 |
call clscdf(cdfid,ierr)
|
|
|
350 |
|
|
|
351 |
c Set the values for the pressure on the pressure levels
|
|
|
352 |
do i=1,pl_nx
|
|
|
353 |
do j=1,pl_ny
|
|
|
354 |
do k=1,pl_nz
|
|
|
355 |
pl_p3(i,j,k)=pl_aklay(k)
|
|
|
356 |
enddo
|
|
|
357 |
enddo
|
|
|
358 |
enddo
|
|
|
359 |
|
|
|
360 |
c -----------------------------------------------------------------
|
|
|
361 |
c Calculate geopotential on layers
|
|
|
362 |
c -----------------------------------------------------------------
|
|
|
363 |
|
|
|
364 |
c Calculate 3d pressure field
|
|
|
365 |
print*,'C P'
|
|
|
366 |
do k=1,ml_nz
|
|
|
367 |
do i=1,ml_nx
|
|
|
368 |
do j=1,ml_ny
|
|
|
369 |
ml_p3(i,j,k)=ml_aklay(k)+ml_bklay(k)*ml_ps(i,j)
|
|
|
370 |
enddo
|
|
|
371 |
enddo
|
|
|
372 |
enddo
|
|
|
373 |
|
|
|
374 |
c Calculate 3d virtual temperature
|
|
|
375 |
print*,'C TV'
|
|
|
376 |
do k=1,ml_nz
|
|
|
377 |
do i=1,ml_nx
|
|
|
378 |
do j=1,ml_ny
|
|
|
379 |
ml_tv3(i,j,k) = (ml_t3(i,j,k)+tzero)*
|
|
|
380 |
> (1.+kappa*ml_q3(i,j,k))
|
|
|
381 |
enddo
|
|
|
382 |
enddo
|
|
|
383 |
enddo
|
|
|
384 |
|
|
|
385 |
c Loop over all grid points
|
|
|
386 |
print*,'C HYDROSTATIC EQUATION'
|
|
|
387 |
do i=1,ml_nx
|
|
|
388 |
do j=1,ml_ny
|
|
|
389 |
|
|
|
390 |
c Make the virtual temperature profile available
|
|
|
391 |
do k=1,ml_nz
|
|
|
392 |
p1 (ml_nz-k+1)=ml_p3 (i,j,k)
|
|
|
393 |
tv1(ml_nz-k+1)=ml_tv3(i,j,k)
|
|
|
394 |
enddo
|
|
|
395 |
call spline (p1,tv1,ml_nz,1.e30,1.e30,spline_tv1)
|
|
|
396 |
|
|
|
397 |
c Loop over all model levels
|
|
|
398 |
do k=1,ml_nz
|
|
|
399 |
|
|
|
400 |
c Get pressure at the grid point
|
|
|
401 |
p = ml_p3(i,j,k)
|
|
|
402 |
|
|
|
403 |
c Find nearest pressure level which is above topography
|
|
|
404 |
lmin=pl_nz
|
|
|
405 |
do l=1,pl_nz
|
|
|
406 |
if ((abs(p-pl_p3(i,j,l))).lt.(abs(p-pl_p3(i,j,lmin)))
|
|
|
407 |
> .and.
|
|
|
408 |
> (pl_p3(i,j,l).lt.ml_ps(i,j)) ) then
|
|
|
409 |
lmin=l
|
|
|
410 |
endif
|
|
|
411 |
enddo
|
|
|
412 |
|
|
|
413 |
c Integrate hydrostatic equation from this level to the grid point
|
|
|
414 |
p0 = pl_p3(i,j,lmin)
|
|
|
415 |
n = nint(abs(p-p0)/dpmin)
|
|
|
416 |
if (n.lt.1) n=1
|
|
|
417 |
dp = (p-p0)/real(n)
|
|
|
418 |
|
|
|
419 |
pu = p0
|
|
|
420 |
z = pl_z3(i,j,lmin)
|
|
|
421 |
call splint(p1,tv1,spline_tv1,ml_nz,pu,tvu)
|
|
|
422 |
do l=1,n
|
|
|
423 |
po = pu+dp
|
|
|
424 |
call splint(p1,tv1,spline_tv1,ml_nz,po,tvo)
|
|
|
425 |
z = z + rdg*0.5*(tvu+tvo)*alog(pu/po)
|
|
|
426 |
tvu = tvo
|
|
|
427 |
pu = po
|
|
|
428 |
enddo
|
|
|
429 |
|
|
|
430 |
c Set the geopotential at the grid point
|
|
|
431 |
ml_z3(i,j,k) = z
|
|
|
432 |
|
|
|
433 |
enddo
|
|
|
434 |
|
|
|
435 |
enddo
|
|
|
436 |
enddo
|
|
|
437 |
|
|
|
438 |
c -----------------------------------------------------------------
|
|
|
439 |
c Calculate height of topography
|
|
|
440 |
c -----------------------------------------------------------------
|
|
|
441 |
|
|
|
442 |
print*,'C TOPOGRAPHY'
|
|
|
443 |
do i=1,ml_nx
|
|
|
444 |
do j=1,ml_ny
|
|
|
445 |
|
|
|
446 |
c Make the z/p profile available
|
|
|
447 |
do k=1,ml_nz
|
|
|
448 |
p1(ml_nz-k+1)=ml_p3(i,j,k)
|
|
|
449 |
z1(ml_nz-k+1)=ml_z3(i,j,k)
|
|
|
450 |
enddo
|
|
|
451 |
|
|
|
452 |
c Cubic spline interpolation
|
|
|
453 |
call spline (p1,z1,ml_nz,1.e30,1.e30,spline_z1)
|
|
|
454 |
call splint (p1,z1,spline_z1,ml_nz,ml_ps(i,j),ml_zb(i,j))
|
|
|
455 |
|
|
|
456 |
enddo
|
|
|
457 |
enddo
|
|
|
458 |
|
|
|
459 |
c -----------------------------------------------------------------
|
|
|
460 |
c Interpolate pressure onto z levels
|
|
|
461 |
c -----------------------------------------------------------------
|
|
|
462 |
|
|
|
463 |
do i=1,zl_nx
|
|
|
464 |
do j=1,zl_ny
|
|
|
465 |
|
|
|
466 |
c Make 1d profile available
|
|
|
467 |
do k=1,ml_nz
|
|
|
468 |
z1(k)=ml_z3(i,j,k)
|
|
|
469 |
f1(k)=ml_p3(i,j,k)
|
|
|
470 |
enddo
|
|
|
471 |
call spline (z1,f1,ml_nz,1.e30,1.e30,spline_f1)
|
|
|
472 |
do k=1,zl_nz
|
|
|
473 |
if (zl_aklay(k).gt.ml_zb(i,j)) then
|
|
|
474 |
call splint(z1,f1,spline_f1,ml_nz,zl_aklay(k),ff)
|
|
|
475 |
zl_p(i,j,k)=ff
|
|
|
476 |
else
|
|
|
477 |
zl_p(i,j,k)=zl_mdv
|
|
|
478 |
endif
|
|
|
479 |
enddo
|
|
|
480 |
|
|
|
481 |
enddo
|
|
|
482 |
enddo
|
|
|
483 |
|
|
|
484 |
c -----------------------------------------------------------------
|
|
|
485 |
c Write output to netcdf file
|
|
|
486 |
c -----------------------------------------------------------------
|
|
|
487 |
|
|
|
488 |
if (test.eq.1) then
|
|
|
489 |
|
|
|
490 |
c Create the output file (grid taken from temperature)
|
|
|
491 |
cdfname=testfile
|
|
|
492 |
|
|
|
493 |
call cdfopn(ml_pfn,cdfid,ierr)
|
|
|
494 |
if (ierr.ne.0) goto 998
|
|
|
495 |
call getcfn(cdfid,cfn,ierr)
|
|
|
496 |
if (ierr.ne.0) goto 998
|
|
|
497 |
call clscdf(cdfid,ierr)
|
|
|
498 |
|
|
|
499 |
call crecdf(cdfname,cdfid,
|
|
|
500 |
> ml_varmin,ml_varmax,ml_ndim,cfn,ierr)
|
|
|
501 |
if (ierr.ne.0) goto 997
|
|
|
502 |
|
|
|
503 |
c Write the definitions of the output variables
|
|
|
504 |
varname='PS'
|
|
|
505 |
ml_vardim(3)=1
|
|
|
506 |
call putdef(cdfid,varname,ml_ndim,ml_mdv,
|
|
|
507 |
> ml_vardim,ml_varmin,ml_varmax,ml_stag,ierr)
|
|
|
508 |
if (ierr.ne.0) goto 997
|
|
|
509 |
varname='ORO'
|
|
|
510 |
call putdef(cdfid,varname,ml_ndim,ml_mdv,
|
|
|
511 |
> ml_vardim,ml_varmin,ml_varmax,ml_stag,ierr)
|
|
|
512 |
if (ierr.ne.0) goto 997
|
|
|
513 |
ml_vardim(3)=ml_nz
|
|
|
514 |
varname='Z'
|
|
|
515 |
call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
|
|
|
516 |
> ml_varmin,ml_varmax,ml_stag,ierr)
|
|
|
517 |
if (ierr.ne.0) goto 997
|
|
|
518 |
|
|
|
519 |
c Write output fields
|
|
|
520 |
print*,'W PS TEST'
|
|
|
521 |
varname='PS'
|
|
|
522 |
call putdat(cdfid,varname,ml_time,1,ml_ps,ierr)
|
|
|
523 |
if (ierr.ne.0) goto 997
|
|
|
524 |
print*,'W ZB TEST'
|
|
|
525 |
varname='ORO'
|
|
|
526 |
call putdat(cdfid,varname,ml_time,1,ml_zb,ierr)
|
|
|
527 |
if (ierr.ne.0) goto 997
|
|
|
528 |
print*,'W Z TEST'
|
|
|
529 |
varname='Z'
|
|
|
530 |
call putdat(cdfid,varname,ml_time,0,ml_z3,ierr)
|
|
|
531 |
if (ierr.ne.0) goto 997
|
|
|
532 |
|
|
|
533 |
c Close cdf file
|
|
|
534 |
call clscdf(cdfid,ierr)
|
|
|
535 |
|
|
|
536 |
endif
|
|
|
537 |
|
|
|
538 |
c -----------------------------------------------------------------
|
|
|
539 |
c Interpolate fields onto a stack of z levels and write new file
|
|
|
540 |
c -----------------------------------------------------------------
|
|
|
541 |
|
|
|
542 |
c Create output file
|
|
|
543 |
cfn=trim(zl_ofn)//'_cst'
|
|
|
544 |
call crecdf(trim(zl_ofn),cdfid,zl_varmin,zl_varmax,
|
|
|
545 |
> zl_ndim,trim(cfn),ierr)
|
|
|
546 |
if (ierr.ne.0) goto 995
|
|
|
547 |
|
|
|
548 |
c Write topography
|
|
|
549 |
print*,'W ZB ',trim(zl_ofn)
|
|
|
550 |
varname='ORO'
|
|
|
551 |
zl_vardim(3)=1
|
|
|
552 |
call putdef(cdfid,varname,zl_ndim,zl_mdv,zl_vardim,
|
|
|
553 |
> zl_varmin,zl_varmax,zl_stag,ierr)
|
|
|
554 |
zl_vardim(3)=zl_nz
|
|
|
555 |
if (ierr.ne.0) goto 995
|
|
|
556 |
call putdat(cdfid,varname,zl_time,1,ml_zb,ierr)
|
|
|
557 |
if (ierr.ne.0) goto 995
|
|
|
558 |
|
|
|
559 |
c Write pressure
|
|
|
560 |
print*,'W P ',trim(zl_ofn)
|
|
|
561 |
varname='P'
|
|
|
562 |
call putdef(cdfid,varname,4,zl_mdv,zl_vardim,
|
|
|
563 |
> zl_varmin,zl_varmax,zl_stag,ierr)
|
|
|
564 |
if (ierr.ne.0) goto 995
|
|
|
565 |
call putdat(cdfid,varname,zl_time,0,zl_p,ierr)
|
|
|
566 |
if (ierr.ne.0) goto 995
|
|
|
567 |
|
|
|
568 |
c Close output file
|
|
|
569 |
call clscdf(cdfid,ierr)
|
|
|
570 |
|
|
|
571 |
c Loop over all variables
|
|
|
572 |
do l=1,nvar
|
|
|
573 |
|
|
|
574 |
print*,'I ',trim(varinp(l))
|
|
|
575 |
|
|
|
576 |
c Get grid description for variable
|
|
|
577 |
call cdfopn(varsrc(l),cdfid,ierr)
|
|
|
578 |
if (ierr.ne.0) goto 996
|
|
|
579 |
call getcfn(cdfid,cfn,ierr)
|
|
|
580 |
if (ierr.ne.0) goto 996
|
|
|
581 |
call cdfopn(cfn,cstid,ierr)
|
|
|
582 |
if (ierr.ne.0) goto 996
|
|
|
583 |
call getvars(cdfid,in_nvars,in_vnam,ierr)
|
|
|
584 |
isok=0
|
|
|
585 |
varname=varinp(l)
|
|
|
586 |
call check_varok (isok,varname, in_vnam,in_nvars)
|
|
|
587 |
if (isok.ne.1) goto 996
|
|
|
588 |
call getdef(cdfid,varinp(l),in_ndim,in_mdv,in_vardim,
|
|
|
589 |
> in_varmin,in_varmax,in_stag,ierr)
|
|
|
590 |
if (ierr.ne.0) goto 996
|
|
|
591 |
in_nx =in_vardim(1)
|
|
|
592 |
in_ny =in_vardim(2)
|
|
|
593 |
in_nz =in_vardim(3)
|
|
|
594 |
in_xmin=in_varmin(1)
|
|
|
595 |
in_ymin=in_varmin(2)
|
|
|
596 |
call getlevs(cstid,in_nz,in_aklev,in_bklev,
|
|
|
597 |
> in_aklay,in_bklay,ierr)
|
|
|
598 |
call getgrid(cstid,in_dx,in_dy,ierr)
|
|
|
599 |
in_xmax=in_xmin+real(in_nx-1)*in_dx
|
|
|
600 |
in_ymax=in_ymin+real(in_ny-1)*in_dy
|
|
|
601 |
call gettimes(cdfid,in_time,in_ntimes,ierr)
|
|
|
602 |
call getstart(cstid,in_idate,ierr)
|
|
|
603 |
call getpole(cstid,in_pollon,in_pollat,ierr)
|
|
|
604 |
call clscdf(cstid,ierr)
|
|
|
605 |
call clscdf(cdfid,ierr)
|
|
|
606 |
|
|
|
607 |
c Check whether this grid is consistent with P grid
|
|
|
608 |
if ( (ml_nx.ne.in_nx).or.
|
|
|
609 |
> (ml_ny.ne.in_ny).or.
|
|
|
610 |
> (ml_nz.ne.in_nz).or.
|
|
|
611 |
> (abs(ml_xmin-in_xmin ).gt.eps).or.
|
|
|
612 |
> (abs(ml_ymin-in_ymin ).gt.eps).or.
|
|
|
613 |
> (abs(ml_xmax-in_xmax ).gt.eps).or.
|
|
|
614 |
> (abs(ml_ymax-in_ymax ).gt.eps).or.
|
|
|
615 |
> (abs(ml_dx -in_dx ).gt.eps).or.
|
|
|
616 |
> (abs(ml_dy -in_dy ).gt.eps).or.
|
|
|
617 |
> (abs(ml_time-in_time ).gt.eps).or.
|
|
|
618 |
> (abs(ml_stag(1)-in_stag(1)).gt.eps).or.
|
|
|
619 |
> (abs(ml_stag(2)-in_stag(2)).gt.eps).or.
|
|
|
620 |
> (abs(ml_stag(3)-in_stag(3)).gt.eps).or.
|
|
|
621 |
> (abs(ml_pollon-in_pollon ).gt.eps).or.
|
|
|
622 |
> (abs(ml_pollat-in_pollat ).gt.eps)) then
|
|
|
623 |
print*,trim(varinp(l)),
|
|
|
624 |
> ' :Input P and FIELD grids are not consistent... Stop'
|
|
|
625 |
print*,'Stag: ',ml_stag, in_stag
|
|
|
626 |
print*,'Pol: ',ml_pollon,in_pollon,ml_pollat,in_pollat
|
|
|
627 |
stop
|
|
|
628 |
endif
|
|
|
629 |
|
|
|
630 |
c Read variable from file
|
|
|
631 |
call cdfopn(varsrc(l),cdfid,ierr)
|
|
|
632 |
if (ierr.ne.0) goto 996
|
|
|
633 |
call getdat(cdfid,varinp(l),in_time,0,in_field,ierr)
|
|
|
634 |
if (ierr.ne.0) goto 996
|
|
|
635 |
call clscdf(cdfid,ierr)
|
|
|
636 |
|
|
|
637 |
c Write the constants file
|
|
|
638 |
if (l.eq.1) then
|
|
|
639 |
|
|
|
640 |
datar(1)=zl_nx
|
|
|
641 |
datar(2)=zl_ny
|
|
|
642 |
datar(3)=int(1000.*zl_varmax(2))
|
|
|
643 |
datar(4)=int(1000.*zl_varmin(1))
|
|
|
644 |
datar(5)=int(1000.*zl_varmin(2))
|
|
|
645 |
datar(6)=int(1000.*zl_varmax(1))
|
|
|
646 |
datar(7)=int(1000.*zl_dx)
|
|
|
647 |
datar(8)=int(1000.*zl_dy)
|
|
|
648 |
datar(9)=zl_nz
|
|
|
649 |
datar(10)=1
|
|
|
650 |
datar(11)=1
|
|
|
651 |
datar(12)=0
|
|
|
652 |
datar(13)=int(1000.*zl_pollon)
|
|
|
653 |
datar(14)=int(1000.*zl_pollat)
|
|
|
654 |
|
|
|
655 |
cfn=trim(zl_ofn)//'_cst'
|
|
|
656 |
call wricst(cfn,datar,zl_aklev,zl_bklev,
|
|
|
657 |
> zl_aklay,zl_bklay,zl_idate)
|
|
|
658 |
|
|
|
659 |
endif
|
|
|
660 |
|
|
|
661 |
c Do the interpolation
|
|
|
662 |
do i=1,zl_nx
|
|
|
663 |
do j=1,zl_ny
|
|
|
664 |
|
|
|
665 |
c Make 1d profile available
|
|
|
666 |
do k=1,ml_nz
|
|
|
667 |
z1(k)=ml_z3 (i,j,k)
|
|
|
668 |
f1(k)=in_field(i,j,k)
|
|
|
669 |
enddo
|
|
|
670 |
call spline (z1,f1,ml_nz,1.e30,1.e30,spline_f1)
|
|
|
671 |
do k=1,zl_nz
|
|
|
672 |
if (zl_aklay(k).gt.ml_zb(i,j)) then
|
|
|
673 |
call splint(z1,f1,spline_f1,ml_nz,zl_aklay(k),ff)
|
|
|
674 |
zl_field(i,j,k)=ff
|
|
|
675 |
else
|
|
|
676 |
zl_field(i,j,k)=zl_mdv
|
|
|
677 |
endif
|
|
|
678 |
enddo
|
|
|
679 |
|
|
|
680 |
enddo
|
|
|
681 |
enddo
|
|
|
682 |
|
|
|
683 |
c Write the output field
|
|
|
684 |
print*,'W ',trim(varout(l)),' ',trim(zl_ofn)
|
|
|
685 |
call cdfwopn(trim(zl_ofn),cdfid,ierr)
|
|
|
686 |
varname=trim(varout(l))
|
|
|
687 |
call putdef(cdfid,varname,4,zl_mdv,zl_vardim,
|
|
|
688 |
> zl_varmin,zl_varmax,zl_stag,ierr)
|
|
|
689 |
if (ierr.ne.0) goto 995
|
|
|
690 |
call putdat(cdfid,varname,zl_time,0,zl_field,ierr)
|
|
|
691 |
if (ierr.ne.0) goto 995
|
|
|
692 |
call clscdf(cdfid,ierr)
|
|
|
693 |
|
|
|
694 |
enddo
|
|
|
695 |
|
|
|
696 |
c -----------------------------------------------------------------
|
|
|
697 |
c Write topography and geopotential also to the input P file
|
|
|
698 |
c -----------------------------------------------------------------
|
|
|
699 |
|
|
|
700 |
c Open the P file
|
|
|
701 |
call cdfwopn(trim(ml_pfn),cdfid,ierr)
|
|
|
702 |
|
|
|
703 |
c Write topography
|
|
|
704 |
varname='ORO'
|
|
|
705 |
print*,'W ',trim(varname),' ',trim(ml_pfn)
|
|
|
706 |
isok=0
|
|
|
707 |
call check_varok (isok,varname,ml_vnam,ml_nvars)
|
|
|
708 |
if (isok.eq.0) then
|
|
|
709 |
ml_vardim(3)=1
|
|
|
710 |
call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
|
|
|
711 |
> ml_varmin,ml_varmax,ml_stag,ierr)
|
|
|
712 |
ml_vardim(3)=ml_nz
|
|
|
713 |
if (ierr.ne.0) goto 997
|
|
|
714 |
endif
|
|
|
715 |
call putdat(cdfid,varname,ml_time,1,ml_zb,ierr)
|
|
|
716 |
if (ierr.ne.0) goto 997
|
|
|
717 |
|
|
|
718 |
c Write geopotential height
|
|
|
719 |
varname='Z'
|
|
|
720 |
print*,'W ',trim(varname),' ',trim(ml_pfn)
|
|
|
721 |
isok=0
|
|
|
722 |
call check_varok (isok,varname,ml_vnam,ml_nvars)
|
|
|
723 |
if (isok.eq.0) then
|
|
|
724 |
call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
|
|
|
725 |
> ml_varmin,ml_varmax,ml_stag,ierr)
|
|
|
726 |
if (ierr.ne.0) goto 997
|
|
|
727 |
endif
|
|
|
728 |
call putdat(cdfid,varname,ml_time,0,ml_z3,ierr)
|
|
|
729 |
if (ierr.ne.0) goto 997
|
|
|
730 |
|
|
|
731 |
c Close P file
|
|
|
732 |
call clscdf(cdfid,ierr)
|
|
|
733 |
|
|
|
734 |
c -----------------------------------------------------------------
|
|
|
735 |
c Exception handling and format specs
|
|
|
736 |
c -----------------------------------------------------------------
|
|
|
737 |
|
|
|
738 |
stop
|
|
|
739 |
|
|
|
740 |
998 print*,'Z: Problems with input from m level'
|
|
|
741 |
stop
|
|
|
742 |
|
|
|
743 |
997 print*,'Z: Problems with output on m level'
|
|
|
744 |
stop
|
|
|
745 |
|
|
|
746 |
996 print*,'F: Problems with input from m level'
|
|
|
747 |
stop
|
|
|
748 |
|
|
|
749 |
995 print*,'F: Problems with output on z level'
|
|
|
750 |
stop
|
|
|
751 |
|
|
|
752 |
|
|
|
753 |
end
|
|
|
754 |
|
|
|
755 |
c -------------------------------------------------------------
|
|
|
756 |
c Natural cubic spline
|
|
|
757 |
c -------------------------------------------------------------
|
|
|
758 |
|
|
|
759 |
SUBROUTINE spline(x,y,n,yp1,ypn,y2)
|
|
|
760 |
INTEGER n,NMAX
|
|
|
761 |
REAL yp1,ypn,x(n),y(n),y2(n)
|
|
|
762 |
PARAMETER (NMAX=500)
|
|
|
763 |
INTEGER i,k
|
|
|
764 |
REAL p,qn,sig,un,u(NMAX)
|
|
|
765 |
if (yp1.gt..99e30) then
|
|
|
766 |
y2(1)=0.
|
|
|
767 |
u(1)=0.
|
|
|
768 |
else
|
|
|
769 |
y2(1)=-0.5
|
|
|
770 |
u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
|
|
|
771 |
endif
|
|
|
772 |
do 11 i=2,n-1
|
|
|
773 |
sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
|
|
|
774 |
p=sig*y2(i-1)+2.
|
|
|
775 |
y2(i)=(sig-1.)/p
|
|
|
776 |
u(i)=(6.*((y(i+1)-y(i))/(x(i+
|
|
|
777 |
*1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
|
|
|
778 |
*u(i-1))/p
|
|
|
779 |
11 continue
|
|
|
780 |
if (ypn.gt..99e30) then
|
|
|
781 |
qn=0.
|
|
|
782 |
un=0.
|
|
|
783 |
else
|
|
|
784 |
qn=0.5
|
|
|
785 |
un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
|
|
|
786 |
endif
|
|
|
787 |
y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
|
|
|
788 |
do 12 k=n-1,1,-1
|
|
|
789 |
y2(k)=y2(k)*y2(k+1)+u(k)
|
|
|
790 |
12 continue
|
|
|
791 |
return
|
|
|
792 |
END
|
|
|
793 |
|
|
|
794 |
|
|
|
795 |
SUBROUTINE splint(xa,ya,y2a,n,x,y)
|
|
|
796 |
INTEGER n
|
|
|
797 |
REAL x,y,xa(n),y2a(n),ya(n)
|
|
|
798 |
INTEGER k,khi,klo
|
|
|
799 |
REAL a,b,h
|
|
|
800 |
klo=1
|
|
|
801 |
khi=n
|
|
|
802 |
1 if (khi-klo.gt.1) then
|
|
|
803 |
k=(khi+klo)/2
|
|
|
804 |
if(xa(k).gt.x)then
|
|
|
805 |
khi=k
|
|
|
806 |
else
|
|
|
807 |
klo=k
|
|
|
808 |
endif
|
|
|
809 |
goto 1
|
|
|
810 |
endif
|
|
|
811 |
h=xa(khi)-xa(klo)
|
|
|
812 |
if (h.eq.0.) pause 'bad xa input in splint'
|
|
|
813 |
a=(xa(khi)-x)/h
|
|
|
814 |
b=(x-xa(klo))/h
|
|
|
815 |
y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
|
|
|
816 |
*2)/6.
|
|
|
817 |
return
|
|
|
818 |
END
|
|
|
819 |
|
|
|
820 |
c ----------------------------------------------------------------
|
|
|
821 |
c Check whether variable is found on netcdf file
|
|
|
822 |
c ----------------------------------------------------------------
|
|
|
823 |
|
|
|
824 |
subroutine check_varok (isok,varname,varlist,nvars)
|
|
|
825 |
|
|
|
826 |
c Check whether the variable <varname> is in the list <varlist(nvars)>.
|
|
|
827 |
c If this is the case, <isok> is incremented by 1. Otherwise <isok>
|
|
|
828 |
c keeps its value.
|
|
|
829 |
|
|
|
830 |
implicit none
|
|
|
831 |
|
|
|
832 |
c Declaraion of subroutine parameters
|
|
|
833 |
integer isok
|
|
|
834 |
integer nvars
|
|
|
835 |
character*80 varname
|
|
|
836 |
character*80 varlist(nvars)
|
|
|
837 |
|
|
|
838 |
c Auxiliary variables
|
|
|
839 |
integer i
|
|
|
840 |
|
|
|
841 |
c Main
|
|
|
842 |
do i=1,nvars
|
|
|
843 |
if (trim(varname).eq.trim(varlist(i))) isok=isok+1
|
|
|
844 |
enddo
|
|
|
845 |
|
|
|
846 |
end
|