| 3 |
michaesp |
1 |
PROGRAM z2s
|
|
|
2 |
|
|
|
3 |
c Calculate secondary fields on z levels
|
|
|
4 |
c Michael Sprenger / Summer 2006
|
|
|
5 |
|
|
|
6 |
implicit none
|
|
|
7 |
|
|
|
8 |
c ---------------------------------------------------------------
|
|
|
9 |
c Declaration of parameters
|
|
|
10 |
c ---------------------------------------------------------------
|
|
|
11 |
|
|
|
12 |
real tzero
|
|
|
13 |
parameter (tzero=273.15)
|
|
|
14 |
real kappa
|
|
|
15 |
parameter (kappa=0.6078)
|
|
|
16 |
real rd
|
|
|
17 |
parameter (rd=287.)
|
|
|
18 |
|
|
|
19 |
c ---------------------------------------------------------------
|
|
|
20 |
c Declaration of variables
|
|
|
21 |
c ---------------------------------------------------------------
|
|
|
22 |
|
|
|
23 |
c Variables for output Z file : height level
|
|
|
24 |
character*80 cfn
|
|
|
25 |
real varmin(4),varmax(4),stag(4)
|
|
|
26 |
integer vardim(4)
|
|
|
27 |
real mdv
|
|
|
28 |
integer ndim
|
|
|
29 |
integer nx,ny,nz
|
|
|
30 |
real xmin,xmax,ymin,ymax,dx,dy
|
|
|
31 |
integer ntimes
|
|
|
32 |
real aklev(1000),bklev(1000)
|
|
|
33 |
real aklay(1000),bklay(1000)
|
|
|
34 |
real time
|
|
|
35 |
real pollon,pollat
|
|
|
36 |
integer idate(5)
|
|
|
37 |
integer nfields
|
|
|
38 |
character*80 field_nam(100)
|
|
|
39 |
real,allocatable, dimension (:,:,:,:) :: field_dat
|
|
|
40 |
real,allocatable, dimension (:,:,:) :: z3
|
|
|
41 |
real,allocatable, dimension (:,:) :: x2,y2,f2,oro
|
|
|
42 |
real,allocatable, dimension (:,:,:) :: out1,out2
|
|
|
43 |
real,allocatable, dimension (:,:,:) :: inp
|
|
|
44 |
real,allocatable,dimension (:) :: nsqref
|
|
|
45 |
real,allocatable,dimension (:) :: thetaref
|
|
|
46 |
real,allocatable,dimension (:) :: tref
|
|
|
47 |
real,allocatable,dimension (:) :: rhoref
|
|
|
48 |
real,allocatable,dimension (:) :: pressref
|
|
|
49 |
real,allocatable,dimension (:) :: zref
|
|
|
50 |
real deltax,deltay,deltaz
|
|
|
51 |
integer nvars
|
|
|
52 |
character*80 vnam(100),varname
|
|
|
53 |
integer isok
|
|
|
54 |
integer cdfid,cstid
|
|
|
55 |
character*3 mode
|
|
|
56 |
|
|
|
57 |
c Parameter file
|
|
|
58 |
character*80 fieldname
|
|
|
59 |
integer nfilt
|
|
|
60 |
character*80 ofn,gri
|
|
|
61 |
|
|
|
62 |
c Auxiliary variables
|
|
|
63 |
integer ierr,stat
|
|
|
64 |
integer i,j,k,n
|
|
|
65 |
real,allocatable, dimension (:,:) :: tmp
|
|
|
66 |
character*80 vnam2(100)
|
|
|
67 |
integer nvars2
|
|
|
68 |
|
|
|
69 |
c ---------------------------------------------------------------
|
|
|
70 |
c Preparations
|
|
|
71 |
c ---------------------------------------------------------------
|
|
|
72 |
|
|
|
73 |
print*,'*********************************************************'
|
|
|
74 |
print*,'* qvec_analysis *'
|
|
|
75 |
print*,'*********************************************************'
|
|
|
76 |
|
|
|
77 |
c Read parameter file
|
|
|
78 |
open(10,file='fort.10')
|
|
|
79 |
read(10,*) fieldname
|
|
|
80 |
read(10,*) ofn
|
|
|
81 |
read(10,*) gri
|
|
|
82 |
read(10,*) nfilt
|
|
|
83 |
close(10)
|
|
|
84 |
|
|
|
85 |
c Decide whether to enter ANO or MOD/ORG mode
|
|
|
86 |
mode = ofn(1:3)
|
|
|
87 |
if ( (mode.ne.'ANO').and.
|
|
|
88 |
> (mode.ne.'MOD').and.
|
|
|
89 |
> (mode.ne.'ORG') )
|
|
|
90 |
>then
|
|
|
91 |
print*,'Unsupported mode ',trim(mode)
|
|
|
92 |
stop
|
|
|
93 |
endif
|
|
|
94 |
|
|
|
95 |
c Get grid description for Z file : height level
|
|
|
96 |
call cdfopn(ofn,cdfid,ierr)
|
|
|
97 |
if (ierr.ne.0) goto 998
|
|
|
98 |
call getcfn(cdfid,cfn,ierr)
|
|
|
99 |
if (ierr.ne.0) goto 998
|
|
|
100 |
call cdfopn(cfn,cstid,ierr)
|
|
|
101 |
if (ierr.ne.0) goto 998
|
|
|
102 |
call getdef(cdfid,'T',ndim,mdv,vardim,
|
|
|
103 |
> varmin,varmax,stag,ierr)
|
|
|
104 |
if (ierr.ne.0) goto 998
|
|
|
105 |
nx =vardim(1)
|
|
|
106 |
ny =vardim(2)
|
|
|
107 |
nz =vardim(3)
|
|
|
108 |
xmin=varmin(1)
|
|
|
109 |
ymin=varmin(2)
|
|
|
110 |
call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)
|
|
|
111 |
if (ierr.ne.0) goto 998
|
|
|
112 |
call getgrid(cstid,dx,dy,ierr)
|
|
|
113 |
if (ierr.ne.0) goto 998
|
|
|
114 |
xmax=xmin+real(nx-1)*dx
|
|
|
115 |
ymax=ymin+real(ny-1)*dy
|
|
|
116 |
call gettimes(cdfid,time,ntimes,ierr)
|
|
|
117 |
if (ierr.ne.0) goto 998
|
|
|
118 |
call getstart(cstid,idate,ierr)
|
|
|
119 |
if (ierr.ne.0) goto 998
|
|
|
120 |
call getpole(cstid,pollon,pollat,ierr)
|
|
|
121 |
if (ierr.ne.0) goto 998
|
|
|
122 |
call getvars(cdfid,nvars,vnam,ierr)
|
|
|
123 |
if (ierr.ne.0) goto 998
|
|
|
124 |
call clscdf(cstid,ierr)
|
|
|
125 |
if (ierr.ne.0) goto 998
|
|
|
126 |
call clscdf(cdfid,ierr)
|
|
|
127 |
if (ierr.ne.0) goto 998
|
|
|
128 |
|
|
|
129 |
c Get a list of all variables on the GRID file
|
|
|
130 |
if (trim(gri).ne.trim(ofn)) then
|
|
|
131 |
call cdfopn(gri,cdfid,ierr)
|
|
|
132 |
if (ierr.ne.0) goto 998
|
|
|
133 |
call getvars(cdfid,nvars2,vnam2,ierr)
|
|
|
134 |
if (ierr.ne.0) goto 998
|
|
|
135 |
do i=1,nvars2
|
|
|
136 |
vnam(nvars+i)=vnam2(i)
|
|
|
137 |
enddo
|
|
|
138 |
nvars=nvars+nvars2
|
|
|
139 |
endif
|
|
|
140 |
|
|
|
141 |
c Check whether calculation of <fieldname> is supported
|
|
|
142 |
if ( (fieldname.ne.'GEO' ).and.
|
|
|
143 |
> (fieldname.ne.'AGEO' ).and.
|
|
|
144 |
> (fieldname.ne.'DIV_UV' ).and.
|
|
|
145 |
> (fieldname.ne.'QVEC' ).and.
|
|
|
146 |
> (fieldname.ne.'DIV_Q' ) ) then
|
|
|
147 |
print*,'Variable not supported ',trim(fieldname)
|
|
|
148 |
stop
|
|
|
149 |
endif
|
|
|
150 |
|
|
|
151 |
c Set dependencies
|
|
|
152 |
if (fieldname.eq.'GEO') then
|
|
|
153 |
nfields=2
|
|
|
154 |
field_nam(1)='T'
|
|
|
155 |
field_nam(2)='P'
|
|
|
156 |
elseif (fieldname.eq.'AGEO') then
|
|
|
157 |
nfields=4
|
|
|
158 |
field_nam(1)='U'
|
|
|
159 |
field_nam(2)='UG'
|
|
|
160 |
field_nam(3)='V'
|
|
|
161 |
field_nam(4)='VG'
|
|
|
162 |
else if (fieldname.eq.'DIV_UV') then
|
|
|
163 |
nfields=2
|
|
|
164 |
field_nam(1)='U'
|
|
|
165 |
field_nam(2)='V'
|
|
|
166 |
else if (fieldname.eq.'QVEC') then
|
|
|
167 |
nfields=3
|
|
|
168 |
field_nam(1)='UG'
|
|
|
169 |
field_nam(2)='VG'
|
|
|
170 |
field_nam(3)='TH'
|
|
|
171 |
else if (fieldname.eq.'DIV_Q') then
|
|
|
172 |
nfields=2
|
|
|
173 |
field_nam(1)='QX'
|
|
|
174 |
field_nam(2)='QY'
|
|
|
175 |
endif
|
|
|
176 |
|
|
|
177 |
c Allocate memory
|
|
|
178 |
allocate(field_dat(nfields,nx,ny,nz),stat=stat)
|
|
|
179 |
if (stat.ne.0) print*,'error allocating field_dat'
|
|
|
180 |
allocate(out1(nx,ny,nz),stat=stat)
|
|
|
181 |
if (stat.ne.0) print*,'error allocating out1'
|
|
|
182 |
allocate(out2(nx,ny,nz),stat=stat)
|
|
|
183 |
if (stat.ne.0) print*,'error allocating out2'
|
|
|
184 |
allocate(inp(nx,ny,nz),stat=stat)
|
|
|
185 |
if (stat.ne.0) print*,'error allocating inp'
|
|
|
186 |
allocate(z3(nx,ny,nz),stat=stat)
|
|
|
187 |
if (stat.ne.0) print*,'error allocating z3'
|
|
|
188 |
allocate(x2(nx,ny),stat=stat)
|
|
|
189 |
if (stat.ne.0) print*,'error allocating x2'
|
|
|
190 |
allocate(y2(nx,ny),stat=stat)
|
|
|
191 |
if (stat.ne.0) print*,'error allocating y2'
|
|
|
192 |
allocate(f2(nx,ny),stat=stat)
|
|
|
193 |
if (stat.ne.0) print*,'error allocating f2'
|
|
|
194 |
allocate(oro(nx,ny),stat=stat)
|
|
|
195 |
if (stat.ne.0) print*,'error allocating oro'
|
|
|
196 |
|
|
|
197 |
C Allocate memory for reference profile
|
|
|
198 |
allocate(rhoref (0:2*nz),STAT=stat)
|
|
|
199 |
if (stat.ne.0) print*,'error allocating rhoref'
|
|
|
200 |
allocate(pressref(0:2*nz),STAT=stat)
|
|
|
201 |
if (stat.ne.0) print*,'error allocating pressref'
|
|
|
202 |
allocate(thetaref(0:2*nz),STAT=stat)
|
|
|
203 |
if (stat.ne.0) print*,'error allocating thetaref'
|
|
|
204 |
allocate(nsqref (0:2*nz),STAT=stat)
|
|
|
205 |
if (stat.ne.0) print*,'error allocating nsqref'
|
|
|
206 |
allocate(zref (0:2*nz),STAT=stat)
|
|
|
207 |
if (stat.ne.0) print*,'error allocating zref'
|
|
|
208 |
allocate(tref (0:2*nz),STAT=stat)
|
|
|
209 |
if (stat.ne.0) print*,'error allocating tref'
|
|
|
210 |
|
|
|
211 |
c Allocate auxiliary fields
|
|
|
212 |
allocate(tmp(nx,ny),stat=stat)
|
|
|
213 |
if (stat.ne.0) print*,'error allocating tmp'
|
|
|
214 |
|
|
|
215 |
c Read reference profile and grid parameters
|
|
|
216 |
call read_ref (nsqref,rhoref,thetaref,pressref,zref,
|
|
|
217 |
> nx,ny,nz,deltax,deltay,deltaz,f2,oro,gri)
|
|
|
218 |
|
|
|
219 |
c Read X grid
|
|
|
220 |
varname='X'
|
|
|
221 |
isok=0
|
|
|
222 |
call check_varok (isok,varname,vnam,nvars)
|
|
|
223 |
if (isok.eq.0) then
|
|
|
224 |
print*,'Variable ',trim(varname),' missing... Stop'
|
|
|
225 |
stop
|
|
|
226 |
endif
|
|
|
227 |
call cdfopn(gri,cdfid,ierr)
|
|
|
228 |
if (ierr.ne.0) goto 998
|
|
|
229 |
call getdat(cdfid,varname,time,0,x2,ierr)
|
|
|
230 |
if (ierr.ne.0) goto 998
|
|
|
231 |
print*,'R ',trim(varname),' ',trim(gri)
|
|
|
232 |
call clscdf(cdfid,ierr)
|
|
|
233 |
if (ierr.ne.0) goto 998
|
|
|
234 |
|
|
|
235 |
c Read Y grid
|
|
|
236 |
varname='Y'
|
|
|
237 |
isok=0
|
|
|
238 |
call check_varok (isok,varname,vnam,nvars)
|
|
|
239 |
if (isok.eq.0) then
|
|
|
240 |
print*,'Variable ',trim(varname),' missing... Stop'
|
|
|
241 |
stop
|
|
|
242 |
endif
|
|
|
243 |
call cdfopn(gri,cdfid,ierr)
|
|
|
244 |
if (ierr.ne.0) goto 998
|
|
|
245 |
call getdat(cdfid,varname,time,0,y2,ierr)
|
|
|
246 |
if (ierr.ne.0) goto 998
|
|
|
247 |
print*,'R ',trim(varname),' ',trim(gri)
|
|
|
248 |
call clscdf(cdfid,ierr)
|
|
|
249 |
if (ierr.ne.0) goto 998
|
|
|
250 |
|
|
|
251 |
c Calculate reference profile of temperature
|
|
|
252 |
print*,'C T_ref = TH_ref * ( P_ref/1000)^kappa'
|
|
|
253 |
do i=0,2*nz
|
|
|
254 |
tref(i) = thetaref(i) * ( pressref(i)/1000. ) ** kappa
|
|
|
255 |
enddo
|
|
|
256 |
|
|
|
257 |
c Init height levels
|
|
|
258 |
do i=1,nx
|
|
|
259 |
do j=1,ny
|
|
|
260 |
do k=1,nz
|
|
|
261 |
z3(i,j,k)=zref(2*k-1)
|
|
|
262 |
enddo
|
|
|
263 |
enddo
|
|
|
264 |
enddo
|
|
|
265 |
|
|
|
266 |
c Load needed variables
|
|
|
267 |
do n=1,nfields
|
|
|
268 |
|
|
|
269 |
c Check whether variable is available on file
|
|
|
270 |
varname=field_nam(n)
|
|
|
271 |
isok=0
|
|
|
272 |
call check_varok (isok,varname,vnam,nvars)
|
|
|
273 |
|
|
|
274 |
c Variable is available on file
|
|
|
275 |
if (isok.eq.1) then
|
|
|
276 |
|
|
|
277 |
call cdfopn(ofn,cdfid,ierr)
|
|
|
278 |
if (ierr.ne.0) goto 998
|
|
|
279 |
call getdat(cdfid,varname,time,0,inp,ierr)
|
|
|
280 |
if (ierr.ne.0) goto 998
|
|
|
281 |
print*,'R ',trim(varname),' ',trim(ofn)
|
|
|
282 |
call clscdf(cdfid,ierr)
|
|
|
283 |
if (ierr.ne.0) goto 998
|
|
|
284 |
|
|
|
285 |
do i=1,nx
|
|
|
286 |
do j=1,ny
|
|
|
287 |
do k=1,nz
|
|
|
288 |
field_dat(n,i,j,k)=inp(i,j,k)
|
|
|
289 |
enddo
|
|
|
290 |
enddo
|
|
|
291 |
enddo
|
|
|
292 |
|
|
|
293 |
else
|
|
|
294 |
print*,'Variable missing : ',trim(varname)
|
|
|
295 |
stop
|
|
|
296 |
endif
|
|
|
297 |
|
|
|
298 |
enddo
|
|
|
299 |
|
|
|
300 |
c Add reference profile if ANO file is provided
|
|
|
301 |
if ( mode.eq.'ANO' ) then
|
|
|
302 |
|
|
|
303 |
print*,'C T -> T_ano + T_ref'
|
|
|
304 |
n=0
|
|
|
305 |
do i=1,nfields
|
|
|
306 |
if (field_nam(i).eq.'T') n=i
|
|
|
307 |
enddo
|
|
|
308 |
if (n.ne.0) then
|
|
|
309 |
do i=1,nx
|
|
|
310 |
do j=1,ny
|
|
|
311 |
do k=1,nz
|
|
|
312 |
field_dat(n,i,j,k) = field_dat(n,i,j,k) + tref(2*k-1)
|
|
|
313 |
enddo
|
|
|
314 |
enddo
|
|
|
315 |
enddo
|
|
|
316 |
endif
|
|
|
317 |
|
|
|
318 |
print*,'C TH -> TH_ano + TH_ref'
|
|
|
319 |
n=0
|
|
|
320 |
do i=1,nfields
|
|
|
321 |
if (field_nam(i).eq.'TH') n=i
|
|
|
322 |
enddo
|
|
|
323 |
if (n.ne.0) then
|
|
|
324 |
do i=1,nx
|
|
|
325 |
do j=1,ny
|
|
|
326 |
do k=1,nz
|
|
|
327 |
field_dat(n,i,j,k) = field_dat(n,i,j,k)+thetaref(2*k-1)
|
|
|
328 |
enddo
|
|
|
329 |
enddo
|
|
|
330 |
enddo
|
|
|
331 |
endif
|
|
|
332 |
|
|
|
333 |
print*,'C P -> P_ano + P_ref'
|
|
|
334 |
n=0
|
|
|
335 |
do i=1,nfields
|
|
|
336 |
if (field_nam(i).eq.'P') n=i
|
|
|
337 |
enddo
|
|
|
338 |
if (n.ne.0) then
|
|
|
339 |
do i=1,nx
|
|
|
340 |
do j=1,ny
|
|
|
341 |
do k=1,nz
|
|
|
342 |
field_dat(n,i,j,k) = field_dat(n,i,j,k)+pressref(2*k-1)
|
|
|
343 |
enddo
|
|
|
344 |
enddo
|
|
|
345 |
enddo
|
|
|
346 |
endif
|
|
|
347 |
|
|
|
348 |
endif
|
|
|
349 |
|
|
|
350 |
c Change units: P (hPa -> Pa), T(K -> degC)
|
|
|
351 |
n=0
|
|
|
352 |
do i=1,nfields
|
|
|
353 |
if (field_nam(i).eq.'P') n=i
|
|
|
354 |
enddo
|
|
|
355 |
if (n.ne.0) then
|
|
|
356 |
do i=1,nx
|
|
|
357 |
do j=1,ny
|
|
|
358 |
do k=1,nz
|
|
|
359 |
field_dat(n,i,j,k)=100.*field_dat(n,i,j,k)
|
|
|
360 |
enddo
|
|
|
361 |
enddo
|
|
|
362 |
enddo
|
|
|
363 |
endif
|
|
|
364 |
|
|
|
365 |
n=0
|
|
|
366 |
do i=1,nfields
|
|
|
367 |
if (field_nam(i).eq.'T') n=i
|
|
|
368 |
enddo
|
|
|
369 |
if (n.ne.0) then
|
|
|
370 |
do i=1,nx
|
|
|
371 |
do j=1,ny
|
|
|
372 |
do k=1,nz
|
|
|
373 |
if ( field_dat(n,i,j,1).gt.100. ) then
|
|
|
374 |
field_dat(n,i,j,k)=field_dat(n,i,j,k) - tzero
|
|
|
375 |
endif
|
|
|
376 |
enddo
|
|
|
377 |
enddo
|
|
|
378 |
enddo
|
|
|
379 |
endif
|
|
|
380 |
|
|
|
381 |
c ----------------------------------------------------------------
|
|
|
382 |
c Calculations
|
|
|
383 |
c ----------------------------------------------------------------
|
|
|
384 |
|
|
|
385 |
c Geostrophic wind (GEO)
|
|
|
386 |
if (fieldname.eq.'GEO') then
|
|
|
387 |
|
|
|
388 |
print*,'C ',trim(fieldname)
|
|
|
389 |
field_nam(1)='RHO'
|
|
|
390 |
do i=1,nx
|
|
|
391 |
do j=1,ny
|
|
|
392 |
do k=1,nz
|
|
|
393 |
if ( mode.eq.'ANO' ) then
|
|
|
394 |
field_dat(1,i,j,k)=rhoref(2*k-1)
|
|
|
395 |
else
|
|
|
396 |
field_dat(1,i,j,k)=
|
|
|
397 |
> 1./rd*field_dat(2,i,j,k)/(field_dat(1,i,j,k)+tzero)
|
|
|
398 |
endif
|
|
|
399 |
enddo
|
|
|
400 |
enddo
|
|
|
401 |
enddo
|
|
|
402 |
call calc_geo (out1,out2, ! (Ug,Vg)
|
|
|
403 |
> field_dat(1,:,:,:), ! RHO
|
|
|
404 |
> field_dat(2,:,:,:), ! P
|
|
|
405 |
> z3,x2,y2,f2, ! Z,X,Y,CORIOL
|
|
|
406 |
> nx,ny,nz,mdv)
|
|
|
407 |
|
|
|
408 |
c Ageostrophic wind (AGEO)
|
|
|
409 |
elseif (fieldname.eq.'AGEO') then
|
|
|
410 |
|
|
|
411 |
print*,'C ',trim(fieldname)
|
|
|
412 |
call calc_ageo (out1,out2, ! (Ua,Va)
|
|
|
413 |
> field_dat(1,:,:,:), ! U
|
|
|
414 |
> field_dat(2,:,:,:), ! UG
|
|
|
415 |
> field_dat(3,:,:,:), ! V
|
|
|
416 |
> field_dat(4,:,:,:), ! VG
|
|
|
417 |
> nx,ny,nz,mdv)
|
|
|
418 |
|
|
|
419 |
c Divergence of wind (DIV_UV)
|
|
|
420 |
else if (fieldname.eq.'DIV_UV') then
|
|
|
421 |
|
|
|
422 |
print*,'C ',trim(fieldname)
|
|
|
423 |
call calc_div_uv (out1, ! Div(U,V))
|
|
|
424 |
> field_dat(1,:,:,:), ! U
|
|
|
425 |
> field_dat(2,:,:,:), ! V
|
|
|
426 |
> z3,x2,y2,f2, ! Z,X,Y,CORIOL
|
|
|
427 |
> nx,ny,nz,mdv)
|
|
|
428 |
|
|
|
429 |
c Q vector (QVEC)
|
|
|
430 |
else if (fieldname.eq.'QVEC') then
|
|
|
431 |
|
|
|
432 |
print*,'C ',trim(fieldname)
|
|
|
433 |
call calc_qvec (out1,out2, ! (Qx,Qy)
|
|
|
434 |
> field_dat(1,:,:,:), ! UG
|
|
|
435 |
> field_dat(2,:,:,:), ! VG
|
|
|
436 |
> field_dat(3,:,:,:), ! TH
|
|
|
437 |
> z3,x2,y2,f2, ! Z,X,Y,CORIOL
|
|
|
438 |
> nx,ny,nz,mdv)
|
|
|
439 |
|
|
|
440 |
c Divergence of Q vector (DIV_Q)
|
|
|
441 |
else if (fieldname.eq.'DIV_Q') then
|
|
|
442 |
|
|
|
443 |
print*,'C ',trim(fieldname)
|
|
|
444 |
call calc_div_q (out1, ! Div(U,V))
|
|
|
445 |
> field_dat(1,:,:,:), ! QX
|
|
|
446 |
> field_dat(2,:,:,:), ! QY
|
|
|
447 |
> z3,x2,y2,f2, ! Z,X,Y,CORIOL
|
|
|
448 |
> nx,ny,nz,mdv)
|
|
|
449 |
|
|
|
450 |
endif
|
|
|
451 |
|
|
|
452 |
c Horizontal filtering of the resulting fields
|
|
|
453 |
print*,'F ',trim(fieldname)
|
|
|
454 |
do k=1,nz
|
|
|
455 |
|
|
|
456 |
do i=1,nx
|
|
|
457 |
do j=1,ny
|
|
|
458 |
tmp(i,j)=out1(i,j,k)
|
|
|
459 |
enddo
|
|
|
460 |
enddo
|
|
|
461 |
do n=1,nfilt
|
|
|
462 |
call filt2d (tmp,tmp,nx,ny,1.,mdv,0,0,0,0)
|
|
|
463 |
enddo
|
|
|
464 |
do i=1,nx
|
|
|
465 |
do j=1,ny
|
|
|
466 |
out1(i,j,k)=tmp(i,j)
|
|
|
467 |
enddo
|
|
|
468 |
enddo
|
|
|
469 |
|
|
|
470 |
do i=1,nx
|
|
|
471 |
do j=1,ny
|
|
|
472 |
tmp(i,j)=out2(i,j,k)
|
|
|
473 |
enddo
|
|
|
474 |
enddo
|
|
|
475 |
do n=1,nfilt
|
|
|
476 |
call filt2d (tmp,tmp,nx,ny,1.,mdv,0,0,0,0)
|
|
|
477 |
enddo
|
|
|
478 |
do i=1,nx
|
|
|
479 |
do j=1,ny
|
|
|
480 |
out2(i,j,k)=tmp(i,j)
|
|
|
481 |
enddo
|
|
|
482 |
enddo
|
|
|
483 |
|
|
|
484 |
enddo
|
|
|
485 |
|
|
|
486 |
c ----------------------------------------------------------------
|
|
|
487 |
c Save result onto netcdf file
|
|
|
488 |
c ----------------------------------------------------------------
|
|
|
489 |
|
|
|
490 |
c Open output file
|
|
|
491 |
call cdfwopn(ofn,cdfid,ierr)
|
|
|
492 |
if (ierr.ne.0) goto 998
|
|
|
493 |
|
|
|
494 |
c Save geostrophic wind
|
|
|
495 |
if (fieldname.eq.'GEO') then
|
|
|
496 |
isok=0
|
|
|
497 |
varname='UG'
|
|
|
498 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
499 |
if (isok.eq.0) then
|
|
|
500 |
call putdef(cdfid,varname,ndim,mdv,vardim,
|
|
|
501 |
> varmin,varmax,stag,ierr)
|
|
|
502 |
if (ierr.ne.0) goto 998
|
|
|
503 |
endif
|
|
|
504 |
call putdat(cdfid,varname,time,0,out1,ierr)
|
|
|
505 |
if (ierr.ne.0) goto 998
|
|
|
506 |
print*,'W ',trim(varname),' ',trim(ofn)
|
|
|
507 |
isok=0
|
|
|
508 |
varname='VG'
|
|
|
509 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
510 |
if (isok.eq.0) then
|
|
|
511 |
call putdef(cdfid,varname,ndim,mdv,vardim,
|
|
|
512 |
> varmin,varmax,stag,ierr)
|
|
|
513 |
if (ierr.ne.0) goto 998
|
|
|
514 |
endif
|
|
|
515 |
call putdat(cdfid,varname,time,0,out2,ierr)
|
|
|
516 |
if (ierr.ne.0) goto 998
|
|
|
517 |
print*,'W ',trim(varname),' ',trim(ofn)
|
|
|
518 |
|
|
|
519 |
c Save ageostrophic wind
|
|
|
520 |
elseif (fieldname.eq.'AGEO') then
|
|
|
521 |
isok=0
|
|
|
522 |
varname='UA'
|
|
|
523 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
524 |
if (isok.eq.0) then
|
|
|
525 |
call putdef(cdfid,varname,ndim,mdv,vardim,
|
|
|
526 |
> varmin,varmax,stag,ierr)
|
|
|
527 |
if (ierr.ne.0) goto 998
|
|
|
528 |
endif
|
|
|
529 |
call putdat(cdfid,varname,time,0,out1,ierr)
|
|
|
530 |
if (ierr.ne.0) goto 998
|
|
|
531 |
print*,'W ',trim(varname),' ',trim(ofn)
|
|
|
532 |
isok=0
|
|
|
533 |
varname='VA'
|
|
|
534 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
535 |
if (isok.eq.0) then
|
|
|
536 |
call putdef(cdfid,varname,ndim,mdv,vardim,
|
|
|
537 |
> varmin,varmax,stag,ierr)
|
|
|
538 |
if (ierr.ne.0) goto 998
|
|
|
539 |
endif
|
|
|
540 |
call putdat(cdfid,varname,time,0,out2,ierr)
|
|
|
541 |
if (ierr.ne.0) goto 998
|
|
|
542 |
print*,'W ',trim(varname),' ',trim(ofn)
|
|
|
543 |
|
|
|
544 |
c Save divergence of wind field
|
|
|
545 |
else if (fieldname.eq.'DIV_UV') then
|
|
|
546 |
isok=0
|
|
|
547 |
varname='DIV_UV'
|
|
|
548 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
549 |
if (isok.eq.0) then
|
|
|
550 |
call putdef(cdfid,varname,ndim,mdv,vardim,
|
|
|
551 |
> varmin,varmax,stag,ierr)
|
|
|
552 |
if (ierr.ne.0) goto 998
|
|
|
553 |
endif
|
|
|
554 |
call putdat(cdfid,varname,time,0,out1,ierr)
|
|
|
555 |
if (ierr.ne.0) goto 998
|
|
|
556 |
print*,'W ',trim(varname),' ',trim(ofn)
|
|
|
557 |
|
|
|
558 |
c Save components of Q vector
|
|
|
559 |
else if (fieldname.eq.'QVEC') then
|
|
|
560 |
isok=0
|
|
|
561 |
varname='QX'
|
|
|
562 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
563 |
if (isok.eq.0) then
|
|
|
564 |
call putdef(cdfid,varname,ndim,mdv,vardim,
|
|
|
565 |
> varmin,varmax,stag,ierr)
|
|
|
566 |
if (ierr.ne.0) goto 998
|
|
|
567 |
endif
|
|
|
568 |
call putdat(cdfid,varname,time,0,out1,ierr)
|
|
|
569 |
if (ierr.ne.0) goto 998
|
|
|
570 |
print*,'W ',trim(varname),' ',trim(ofn)
|
|
|
571 |
isok=0
|
|
|
572 |
varname='QY'
|
|
|
573 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
574 |
if (isok.eq.0) then
|
|
|
575 |
call putdef(cdfid,varname,ndim,mdv,vardim,
|
|
|
576 |
> varmin,varmax,stag,ierr)
|
|
|
577 |
if (ierr.ne.0) goto 998
|
|
|
578 |
endif
|
|
|
579 |
call putdat(cdfid,varname,time,0,out2,ierr)
|
|
|
580 |
if (ierr.ne.0) goto 998
|
|
|
581 |
print*,'W ',trim(varname),' ',trim(ofn)
|
|
|
582 |
|
|
|
583 |
c Save divergence of wind field
|
|
|
584 |
else if (fieldname.eq.'DIV_Q') then
|
|
|
585 |
isok=0
|
|
|
586 |
varname='DIV_Q'
|
|
|
587 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
588 |
if (isok.eq.0) then
|
|
|
589 |
call putdef(cdfid,varname,ndim,mdv,vardim,
|
|
|
590 |
> varmin,varmax,stag,ierr)
|
|
|
591 |
if (ierr.ne.0) goto 998
|
|
|
592 |
endif
|
|
|
593 |
call putdat(cdfid,varname,time,0,out1,ierr)
|
|
|
594 |
if (ierr.ne.0) goto 998
|
|
|
595 |
print*,'W ',trim(varname),' ',trim(ofn)
|
|
|
596 |
|
|
|
597 |
endif
|
|
|
598 |
|
|
|
599 |
c Close output file
|
|
|
600 |
call clscdf(cdfid,ierr)
|
|
|
601 |
if (ierr.ne.0) goto 998
|
|
|
602 |
|
|
|
603 |
|
|
|
604 |
c ----------------------------------------------------------------
|
|
|
605 |
c Exception handling
|
|
|
606 |
c ----------------------------------------------------------------
|
|
|
607 |
|
|
|
608 |
stop
|
|
|
609 |
|
|
|
610 |
998 print*,'Problem with netcdf file'
|
|
|
611 |
stop
|
|
|
612 |
|
|
|
613 |
end
|
|
|
614 |
|
|
|
615 |
|
|
|
616 |
|
|
|
617 |
c ****************************************************************
|
|
|
618 |
c * SUBROUTINE SECTION: AUXILIARY ROUTINES *
|
|
|
619 |
c ****************************************************************
|
|
|
620 |
|
|
|
621 |
c ----------------------------------------------------------------
|
|
|
622 |
c Check whether variable is found on netcdf file
|
|
|
623 |
c ----------------------------------------------------------------
|
|
|
624 |
|
|
|
625 |
subroutine check_varok (isok,varname,varlist,nvars)
|
|
|
626 |
|
|
|
627 |
c Check whether the variable <varname> is in the list <varlist(nvars)>.
|
|
|
628 |
c If this is the case, <isok> is incremented by 1. Otherwise <isok>
|
|
|
629 |
c keeps its value.
|
|
|
630 |
|
|
|
631 |
implicit none
|
|
|
632 |
|
|
|
633 |
c Declaraion of subroutine parameters
|
|
|
634 |
integer isok
|
|
|
635 |
integer nvars
|
|
|
636 |
character*80 varname
|
|
|
637 |
character*80 varlist(nvars)
|
|
|
638 |
|
|
|
639 |
c Auxiliary variables
|
|
|
640 |
integer i
|
|
|
641 |
|
|
|
642 |
c Main
|
|
|
643 |
do i=1,nvars
|
|
|
644 |
if (trim(varname).eq.trim(varlist(i))) isok=isok+1
|
|
|
645 |
enddo
|
|
|
646 |
|
|
|
647 |
end
|
|
|
648 |
|
|
|
649 |
|
|
|
650 |
c ****************************************************************
|
|
|
651 |
c * SUBROUTINE SECTION: CALCULATE SECONDARY FIELDS *
|
|
|
652 |
c ****************************************************************
|
|
|
653 |
|
|
|
654 |
c ----------------------------------------------------------------
|
|
|
655 |
c Calculate geostrophic wind components
|
|
|
656 |
c ----------------------------------------------------------------
|
|
|
657 |
|
|
|
658 |
subroutine calc_geo (ug,vg,rho,p,
|
|
|
659 |
> z3,x2,y2,f2,nx,ny,nz,mdv)
|
|
|
660 |
|
|
|
661 |
c Calculate the geostrophic wind components (ug,vg) if the temperature
|
|
|
662 |
c (t) and the pressure (p) are given. The grid and the missing data
|
|
|
663 |
c value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv)
|
|
|
664 |
|
|
|
665 |
implicit none
|
|
|
666 |
|
|
|
667 |
c Declaration of subroutine parameters
|
|
|
668 |
integer nx,ny,nz
|
|
|
669 |
real ug (nx,ny,nz)
|
|
|
670 |
real vg (nx,ny,nz)
|
|
|
671 |
real rho(nx,ny,nz)
|
|
|
672 |
real p (nx,ny,nz)
|
|
|
673 |
real z3 (nx,ny,nz)
|
|
|
674 |
real x2 (nx,ny)
|
|
|
675 |
real y2 (nx,ny)
|
|
|
676 |
real f2 (nx,ny)
|
|
|
677 |
real mdv
|
|
|
678 |
|
|
|
679 |
c Physical parameters and numerical constants
|
|
|
680 |
real eps
|
|
|
681 |
parameter (eps=0.01)
|
|
|
682 |
real g
|
|
|
683 |
parameter (g=9.80616)
|
|
|
684 |
|
|
|
685 |
|
|
|
686 |
c Auxiliray variables
|
|
|
687 |
integer i,j,k
|
|
|
688 |
real dpdx(nx,ny,nz)
|
|
|
689 |
real dpdy(nx,ny,nz)
|
|
|
690 |
|
|
|
691 |
c Calculate horizontal derivatives of pressure
|
|
|
692 |
call deriv(dpdx,p,'x',z3,x2,y2,nx,ny,nz,mdv)
|
|
|
693 |
call deriv(dpdy,p,'y',z3,x2,y2,nx,ny,nz,mdv)
|
|
|
694 |
|
|
|
695 |
c Calculation
|
|
|
696 |
do i=1,nx
|
|
|
697 |
do j=1,ny
|
|
|
698 |
do k=1,nz
|
|
|
699 |
|
|
|
700 |
if ((abs(rho (i,j,k)-mdv).gt.eps).and.
|
|
|
701 |
> (abs(dpdx(i,j,k)-mdv).gt.eps).and.
|
|
|
702 |
> (abs(dpdy(i,j,k)-mdv).gt.eps)) then
|
|
|
703 |
|
|
|
704 |
ug(i,j,k)=-1./(rho(i,j,k)*f2(i,j))*dpdy(i,j,k)
|
|
|
705 |
vg(i,j,k)= 1./(rho(i,j,k)*f2(i,j))*dpdx(i,j,k)
|
|
|
706 |
|
|
|
707 |
endif
|
|
|
708 |
|
|
|
709 |
enddo
|
|
|
710 |
enddo
|
|
|
711 |
enddo
|
|
|
712 |
|
|
|
713 |
end
|
|
|
714 |
|
|
|
715 |
c ----------------------------------------------------------------
|
|
|
716 |
c Calculate ageostrophic wind components
|
|
|
717 |
c ----------------------------------------------------------------
|
|
|
718 |
|
|
|
719 |
subroutine calc_ageo (ua,va,u,ug,v,vg,nx,ny,nz,mdv)
|
|
|
720 |
|
|
|
721 |
c Calculate the geostrophic wind components (ug,vg) if the temperature
|
|
|
722 |
c (t) and the pressure (p) are given. The grid and the missing data
|
|
|
723 |
c value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv)
|
|
|
724 |
|
|
|
725 |
implicit none
|
|
|
726 |
|
|
|
727 |
c Declaration of subroutine parameters
|
|
|
728 |
integer nx,ny,nz
|
|
|
729 |
real ua (nx,ny,nz)
|
|
|
730 |
real va (nx,ny,nz)
|
|
|
731 |
real u (nx,ny,nz)
|
|
|
732 |
real ug (nx,ny,nz)
|
|
|
733 |
real v (nx,ny,nz)
|
|
|
734 |
real vg (nx,ny,nz)
|
|
|
735 |
real mdv
|
|
|
736 |
|
|
|
737 |
c Parameters
|
|
|
738 |
real eps
|
|
|
739 |
parameter (eps=0.01)
|
|
|
740 |
|
|
|
741 |
c Auxiliray variables
|
|
|
742 |
integer i,j,k
|
|
|
743 |
|
|
|
744 |
c Calculation
|
|
|
745 |
do i=1,nx
|
|
|
746 |
do j=1,ny
|
|
|
747 |
do k=1,nz
|
|
|
748 |
|
|
|
749 |
if ((abs(u (i,j,k)-mdv).gt.eps).and.
|
|
|
750 |
> (abs(ug(i,j,k)-mdv).gt.eps).and.
|
|
|
751 |
> (abs(v (i,j,k)-mdv).gt.eps).and.
|
|
|
752 |
> (abs(vg(i,j,k)-mdv).gt.eps)) then
|
|
|
753 |
|
|
|
754 |
ua(i,j,k)=u(i,j,k) - ug(i,j,k)
|
|
|
755 |
va(i,j,k)=v(i,j,k) - vg(i,j,k)
|
|
|
756 |
|
|
|
757 |
endif
|
|
|
758 |
|
|
|
759 |
enddo
|
|
|
760 |
enddo
|
|
|
761 |
enddo
|
|
|
762 |
|
|
|
763 |
end
|
|
|
764 |
|
|
|
765 |
c ----------------------------------------------------------------
|
|
|
766 |
c Calculate divergence of wind field
|
|
|
767 |
c ----------------------------------------------------------------
|
|
|
768 |
|
|
|
769 |
subroutine calc_div_uv (div_uv,u,v,
|
|
|
770 |
> z3,x2,y2,f2,nx,ny,nz,mdv)
|
|
|
771 |
|
|
|
772 |
c Calculate the divergence (div_uv) of the horizontal wind field if+
|
|
|
773 |
c the wind components (u,v) are given. The grid and the missing data
|
|
|
774 |
c value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv)
|
|
|
775 |
|
|
|
776 |
implicit none
|
|
|
777 |
|
|
|
778 |
c Declaration of subroutine parameters
|
|
|
779 |
integer nx,ny,nz
|
|
|
780 |
real div_uv(nx,ny,nz)
|
|
|
781 |
real u (nx,ny,nz)
|
|
|
782 |
real v (nx,ny,nz)
|
|
|
783 |
real z3 (nx,ny,nz)
|
|
|
784 |
real x2 (nx,ny)
|
|
|
785 |
real y2 (nx,ny)
|
|
|
786 |
real f2 (nx,ny)
|
|
|
787 |
real mdv
|
|
|
788 |
|
|
|
789 |
c Physical parameters and numerical constants
|
|
|
790 |
real eps
|
|
|
791 |
parameter (eps=0.01)
|
|
|
792 |
|
|
|
793 |
c Auxiliray variables
|
|
|
794 |
integer i,j,k
|
|
|
795 |
real rho
|
|
|
796 |
real dudx(nx,ny,nz)
|
|
|
797 |
real dvdy(nx,ny,nz)
|
|
|
798 |
|
|
|
799 |
c Calculate horizontal derivatives of pressure
|
|
|
800 |
call deriv(dudx,u,'x',z3,x2,y2,nx,ny,nz,mdv)
|
|
|
801 |
call deriv(dvdy,v,'y',z3,x2,y2,nx,ny,nz,mdv)
|
|
|
802 |
|
|
|
803 |
c Calculation
|
|
|
804 |
do i=1,nx
|
|
|
805 |
do j=1,ny
|
|
|
806 |
do k=1,nz
|
|
|
807 |
|
|
|
808 |
if ((abs(dudx(i,j,k)-mdv).gt.eps).and.
|
|
|
809 |
> (abs(dvdy(i,j,k)-mdv).gt.eps)) then
|
|
|
810 |
|
|
|
811 |
div_uv(i,j,k)=dudx(i,j,k)+dvdy(i,j,k)
|
|
|
812 |
|
|
|
813 |
endif
|
|
|
814 |
|
|
|
815 |
enddo
|
|
|
816 |
enddo
|
|
|
817 |
enddo
|
|
|
818 |
|
|
|
819 |
end
|
|
|
820 |
|
|
|
821 |
c ----------------------------------------------------------------
|
|
|
822 |
c Calculate Q vector
|
|
|
823 |
c ----------------------------------------------------------------
|
|
|
824 |
|
|
|
825 |
subroutine calc_qvec (qx3,qy3,
|
|
|
826 |
> th3,u3,v3,
|
|
|
827 |
> z3,x2,y2,f2,nx,ny,nz,mdv)
|
|
|
828 |
|
|
|
829 |
c Calculate teh Q vector components <qx3> and <qy3>, as well as the divergence
|
|
|
830 |
c of the Q vector <divq3> on the model grid. The grid is specified in the horizontal
|
|
|
831 |
c by <xmin,ymin,dx,dy,nx,ny>. The number of vertical levels is <nz>. The input field
|
|
|
832 |
c are: potential temperature <th3>, horizontal wind <u3> and <v3>, pressure <p3>.
|
|
|
833 |
c The calculation follows the one described in "Weather Analysis, Dusan Djuric"
|
|
|
834 |
|
|
|
835 |
implicit none
|
|
|
836 |
|
|
|
837 |
c Declaration of subroutine parameters
|
|
|
838 |
integer nx,ny,nz
|
|
|
839 |
real qx3(nx,ny,nz)
|
|
|
840 |
real qy3(nx,ny,nz)
|
|
|
841 |
real th3(nx,ny,nz)
|
|
|
842 |
real u3 (nx,ny,nz)
|
|
|
843 |
real v3 (nx,ny,nz)
|
|
|
844 |
real z3 (nx,ny,nz)
|
|
|
845 |
real x2 (nx,ny)
|
|
|
846 |
real y2 (nx,ny)
|
|
|
847 |
real f2 (nx,ny)
|
|
|
848 |
real mdv
|
|
|
849 |
|
|
|
850 |
c Physical and numerical parameters
|
|
|
851 |
real scale1,scale2
|
|
|
852 |
parameter (scale1=1.E10,scale2=1.E14)
|
|
|
853 |
real eps
|
|
|
854 |
parameter (eps=0.01)
|
|
|
855 |
real g
|
|
|
856 |
parameter (g=9.80616)
|
|
|
857 |
real tzero
|
|
|
858 |
parameter (tzero=273.16)
|
|
|
859 |
|
|
|
860 |
c Auxiliary variables
|
|
|
861 |
real dudx(nx,ny,nz)
|
|
|
862 |
real dudy(nx,ny,nz)
|
|
|
863 |
real dvdx(nx,ny,nz)
|
|
|
864 |
real dvdy(nx,ny,nz)
|
|
|
865 |
real dtdx(nx,ny,nz)
|
|
|
866 |
real dtdy(nx,ny,nz)
|
|
|
867 |
integer i,j,k
|
|
|
868 |
|
|
|
869 |
c Needed derivatives
|
|
|
870 |
call deriv (dudx, u3,'x',z3,x2,y2,nx,ny,nz,mdv)
|
|
|
871 |
call deriv (dudy, u3,'y',z3,x2,y2,nx,ny,nz,mdv)
|
|
|
872 |
call deriv (dvdx, v3,'x',z3,x2,y2,nx,ny,nz,mdv)
|
|
|
873 |
call deriv (dvdy, v3,'y',z3,x2,y2,nx,ny,nz,mdv)
|
|
|
874 |
call deriv (dtdx,th3,'x',z3,x2,y2,nx,ny,nz,mdv)
|
|
|
875 |
call deriv (dtdy,th3,'y',z3,x2,y2,nx,ny,nz,mdv)
|
|
|
876 |
|
|
|
877 |
c Calculate vector components of the Q vector
|
|
|
878 |
do i=1,nx
|
|
|
879 |
do j=1,ny
|
|
|
880 |
do k=1,nz
|
|
|
881 |
|
|
|
882 |
c Evaluate Q vector formula with missing data check
|
|
|
883 |
if ((abs(dudx(i,j,k)-mdv).gt.eps).and.
|
|
|
884 |
> (abs(dudy(i,j,k)-mdv).gt.eps).and.
|
|
|
885 |
> (abs(dvdx(i,j,k)-mdv).gt.eps).and.
|
|
|
886 |
> (abs(dvdy(i,j,k)-mdv).gt.eps).and.
|
|
|
887 |
> (abs(dtdx(i,j,k)-mdv).gt.eps).and.
|
|
|
888 |
> (abs(dtdy(i,j,k)-mdv).gt.eps)) then
|
|
|
889 |
|
|
|
890 |
qx3(i,j,k) = -g/tzero * (dudx(i,j,k)*dtdx(i,j,k)
|
|
|
891 |
> +dvdx(i,j,k)*dtdy(i,j,k))
|
|
|
892 |
qy3(i,j,k) = -g/tzero * (dudy(i,j,k)*dtdx(i,j,k)
|
|
|
893 |
> +dvdy(i,j,k)*dtdy(i,j,k))
|
|
|
894 |
|
|
|
895 |
else
|
|
|
896 |
qx3(i,j,k)=mdv
|
|
|
897 |
qy3(i,j,k)=mdv
|
|
|
898 |
endif
|
|
|
899 |
|
|
|
900 |
enddo
|
|
|
901 |
enddo
|
|
|
902 |
enddo
|
|
|
903 |
|
|
|
904 |
c Scale the output
|
|
|
905 |
do i=1,nx
|
|
|
906 |
do j=1,ny
|
|
|
907 |
do k=1,nz
|
|
|
908 |
if (abs(qx3(i,j,k)-mdv).gt.eps) then
|
|
|
909 |
qx3(i,j,k)=scale1*qx3(i,j,k)
|
|
|
910 |
endif
|
|
|
911 |
if (abs(qy3(i,j,k)-mdv).gt.eps) then
|
|
|
912 |
qy3(i,j,k)=scale1*qy3(i,j,k)
|
|
|
913 |
endif
|
|
|
914 |
enddo
|
|
|
915 |
enddo
|
|
|
916 |
enddo
|
|
|
917 |
|
|
|
918 |
end
|
|
|
919 |
|
|
|
920 |
c ----------------------------------------------------------------
|
|
|
921 |
c Calculate divergence of wind field
|
|
|
922 |
c ----------------------------------------------------------------
|
|
|
923 |
|
|
|
924 |
subroutine calc_div_q (div_q,qx,qy,
|
|
|
925 |
> z3,x2,y2,f2,nx,ny,nz,mdv)
|
|
|
926 |
|
|
|
927 |
c Calculate the divergence (div_q) of the Q vector field if
|
|
|
928 |
c the components (qx,qy) are given. The grid and the missing data
|
|
|
929 |
c value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv)
|
|
|
930 |
|
|
|
931 |
implicit none
|
|
|
932 |
|
|
|
933 |
c Declaration of subroutine parameters
|
|
|
934 |
integer nx,ny,nz
|
|
|
935 |
real div_q(nx,ny,nz)
|
|
|
936 |
real qx (nx,ny,nz)
|
|
|
937 |
real qy (nx,ny,nz)
|
|
|
938 |
real z3 (nx,ny,nz)
|
|
|
939 |
real x2 (nx,ny)
|
|
|
940 |
real y2 (nx,ny)
|
|
|
941 |
real f2 (nx,ny)
|
|
|
942 |
real mdv
|
|
|
943 |
|
|
|
944 |
c Physical parameters and numerical constants
|
|
|
945 |
real eps
|
|
|
946 |
parameter (eps=0.01)
|
|
|
947 |
|
|
|
948 |
c Auxiliray variables
|
|
|
949 |
integer i,j,k
|
|
|
950 |
real rho
|
|
|
951 |
real dqxdx(nx,ny,nz)
|
|
|
952 |
real dqydy(nx,ny,nz)
|
|
|
953 |
|
|
|
954 |
c Calculate horizontal derivatives of pressure
|
|
|
955 |
call deriv(dqxdx,qx,'x',z3,x2,y2,nx,ny,nz,mdv)
|
|
|
956 |
call deriv(dqydy,qy,'y',z3,x2,y2,nx,ny,nz,mdv)
|
|
|
957 |
|
|
|
958 |
c Calculation
|
|
|
959 |
do i=1,nx
|
|
|
960 |
do j=1,ny
|
|
|
961 |
do k=1,nz
|
|
|
962 |
|
|
|
963 |
if ((abs(dqxdx(i,j,k)-mdv).gt.eps).and.
|
|
|
964 |
> (abs(dqydy(i,j,k)-mdv).gt.eps)) then
|
|
|
965 |
|
|
|
966 |
div_q(i,j,k)=dqxdx(i,j,k)+dqydy(i,j,k)
|
|
|
967 |
|
|
|
968 |
endif
|
|
|
969 |
|
|
|
970 |
enddo
|
|
|
971 |
enddo
|
|
|
972 |
enddo
|
|
|
973 |
|
|
|
974 |
end
|
|
|
975 |
|
|
|
976 |
|
|
|
977 |
c ****************************************************************
|
|
|
978 |
c * SUBROUTINE SECTION: GRID HANDLING *
|
|
|
979 |
c ****************************************************************
|
|
|
980 |
|
|
|
981 |
c -----------------------------------------------------------------
|
|
|
982 |
c Horizontal and vertical derivatives for 3d fields
|
|
|
983 |
c -----------------------------------------------------------------
|
|
|
984 |
|
|
|
985 |
subroutine deriv (df,f,direction,
|
|
|
986 |
> z3,x2,y2,nx,ny,nz,mdv)
|
|
|
987 |
|
|
|
988 |
c Calculate horizontal and vertical derivatives of the 3d field <f>.
|
|
|
989 |
c The direction of the derivative is specified in <direction>
|
|
|
990 |
c 'x','y' : Horizontal derivative in x and y direction
|
|
|
991 |
c 'p','z','t','m' : Vertical derivative (pressure, height, theta, model)
|
|
|
992 |
c The 3d field <z3> specifies the isosurfaces along which the horizontal
|
|
|
993 |
c derivatives are calculated or the levels for the vertical derivatives.
|
|
|
994 |
|
|
|
995 |
implicit none
|
|
|
996 |
|
|
|
997 |
c Input and output parameters
|
|
|
998 |
integer nx,ny,nz
|
|
|
999 |
real df (nx,ny,nz)
|
|
|
1000 |
real f (nx,ny,nz)
|
|
|
1001 |
real z3 (nx,ny,nz)
|
|
|
1002 |
real x2 (nx,ny)
|
|
|
1003 |
real y2 (nx,ny)
|
|
|
1004 |
character direction
|
|
|
1005 |
real mdv
|
|
|
1006 |
|
|
|
1007 |
c Numerical and physical parameters
|
|
|
1008 |
real pi180
|
|
|
1009 |
parameter (pi180=3.141592654/180.)
|
|
|
1010 |
real deltay
|
|
|
1011 |
parameter (deltay=111.1775E3)
|
|
|
1012 |
real zerodiv
|
|
|
1013 |
parameter (zerodiv=0.00000001)
|
|
|
1014 |
real eps
|
|
|
1015 |
parameter (eps=0.01)
|
|
|
1016 |
|
|
|
1017 |
c Auxiliary variables
|
|
|
1018 |
integer i,j,k
|
|
|
1019 |
real vmin,vmax
|
|
|
1020 |
real scale,lat
|
|
|
1021 |
real vu,vl,vuvl,vlvu
|
|
|
1022 |
integer o,u,w,e,n,s
|
|
|
1023 |
|
|
|
1024 |
c Vertical derivative
|
|
|
1025 |
if ((direction.eq.'z').or.
|
|
|
1026 |
> (direction.eq.'th').or.
|
|
|
1027 |
> (direction.eq.'p').or.
|
|
|
1028 |
> (direction.eq.'m').and.
|
|
|
1029 |
> (nz.gt.1)) then
|
|
|
1030 |
|
|
|
1031 |
do i=1,nx
|
|
|
1032 |
do j=1,ny
|
|
|
1033 |
do k=1,nz
|
|
|
1034 |
|
|
|
1035 |
o=k+1
|
|
|
1036 |
if (o.gt.nz) o=nz
|
|
|
1037 |
u=k-1
|
|
|
1038 |
if (u.lt.1) u=1
|
|
|
1039 |
|
|
|
1040 |
if ((abs(f(i,j,o)-mdv).gt.eps).and.
|
|
|
1041 |
> (abs(f(i,j,u)-mdv).gt.eps).and.
|
|
|
1042 |
> (abs(f(i,j,k)-mdv).gt.eps)) then
|
|
|
1043 |
|
|
|
1044 |
vu = z3(i,j,k)-z3(i,j,o)
|
|
|
1045 |
vl = z3(i,j,u)-z3(i,j,k)
|
|
|
1046 |
vuvl = vu/(vl+zerodiv)
|
|
|
1047 |
vlvu = 1./(vuvl+zerodiv)
|
|
|
1048 |
|
|
|
1049 |
df(i,j,k) = 1./(vu+vl)
|
|
|
1050 |
> * (vuvl*(f(i,j,u)-f(i,j,k))
|
|
|
1051 |
> + vlvu*(f(i,j,k)-f(i,j,o)))
|
|
|
1052 |
|
|
|
1053 |
else
|
|
|
1054 |
df(i,j,k) = mdv
|
|
|
1055 |
endif
|
|
|
1056 |
|
|
|
1057 |
enddo
|
|
|
1058 |
enddo
|
|
|
1059 |
enddo
|
|
|
1060 |
|
|
|
1061 |
c Horizontal derivative in the y direction: 3d
|
|
|
1062 |
elseif (direction.eq.'y') then
|
|
|
1063 |
|
|
|
1064 |
do i=1,nx
|
|
|
1065 |
do j=1,ny
|
|
|
1066 |
do k=1,nz
|
|
|
1067 |
|
|
|
1068 |
s=j-1
|
|
|
1069 |
if (s.lt.1) s=1
|
|
|
1070 |
n=j+1
|
|
|
1071 |
if (n.gt.ny) n=ny
|
|
|
1072 |
|
|
|
1073 |
if ((abs(f(i,n,k)-mdv).gt.eps).and.
|
|
|
1074 |
> (abs(f(i,j,k)-mdv).gt.eps).and.
|
|
|
1075 |
> (abs(f(i,s,k)-mdv).gt.eps)) then
|
|
|
1076 |
|
|
|
1077 |
vu = 1000.*(y2(i,j)-y2(i,n))
|
|
|
1078 |
vl = 1000.*(y2(i,s)-y2(i,j))
|
|
|
1079 |
vuvl = vu/(vl+zerodiv)
|
|
|
1080 |
vlvu = 1./(vuvl+zerodiv)
|
|
|
1081 |
|
|
|
1082 |
df(i,j,k) = 1./(vu+vl)
|
|
|
1083 |
> * (vuvl*(f(i,s,k)-f(i,j,k))
|
|
|
1084 |
> + vlvu*(f(i,j,k)-f(i,n,k)))
|
|
|
1085 |
|
|
|
1086 |
else
|
|
|
1087 |
df(i,j,k) = mdv
|
|
|
1088 |
endif
|
|
|
1089 |
|
|
|
1090 |
enddo
|
|
|
1091 |
enddo
|
|
|
1092 |
enddo
|
|
|
1093 |
|
|
|
1094 |
c Horizontal derivative in the x direction: 3d
|
|
|
1095 |
elseif (direction.eq.'x') then
|
|
|
1096 |
|
|
|
1097 |
do i=1,nx
|
|
|
1098 |
do j=1,ny
|
|
|
1099 |
do k=1,nz
|
|
|
1100 |
|
|
|
1101 |
w=i-1
|
|
|
1102 |
if (w.lt.1) w=1
|
|
|
1103 |
e=i+1
|
|
|
1104 |
if (e.gt.nx) e=nx
|
|
|
1105 |
|
|
|
1106 |
if ((abs(f(w,j,k)-mdv).gt.eps).and.
|
|
|
1107 |
> (abs(f(i,j,k)-mdv).gt.eps).and.
|
|
|
1108 |
> (abs(f(e,j,k)-mdv).gt.eps)) then
|
|
|
1109 |
|
|
|
1110 |
vu = 1000.*(x2(i,j)-x2(e,j))
|
|
|
1111 |
vl = 1000.*(x2(w,j)-x2(i,j))
|
|
|
1112 |
vuvl = vu/(vl+zerodiv)
|
|
|
1113 |
vlvu = 1./(vuvl+zerodiv)
|
|
|
1114 |
|
|
|
1115 |
df(i,j,k) = 1./(vu+vl)
|
|
|
1116 |
> * (vuvl*(f(w,j,k)-f(i,j,k))
|
|
|
1117 |
> + vlvu*(f(i,j,k)-f(e,j,k)))
|
|
|
1118 |
|
|
|
1119 |
else
|
|
|
1120 |
df(i,j,k) = mdv
|
|
|
1121 |
endif
|
|
|
1122 |
|
|
|
1123 |
enddo
|
|
|
1124 |
enddo
|
|
|
1125 |
enddo
|
|
|
1126 |
|
|
|
1127 |
c Undefined direction for derivative
|
|
|
1128 |
else
|
|
|
1129 |
|
|
|
1130 |
print*,'Invalid direction of derivative... Stop'
|
|
|
1131 |
stop
|
|
|
1132 |
|
|
|
1133 |
endif
|
|
|
1134 |
|
|
|
1135 |
end
|
|
|
1136 |
|
|
|
1137 |
c -----------------------------------------------------------------
|
|
|
1138 |
c Horizontal filter
|
|
|
1139 |
c -----------------------------------------------------------------
|
|
|
1140 |
|
|
|
1141 |
subroutine filt2d (a,af,nx,ny,fil,misdat,
|
|
|
1142 |
& iperx,ipery,ispol,inpol)
|
|
|
1143 |
|
|
|
1144 |
c Apply a conservative diffusion operator onto the 2d field a,
|
|
|
1145 |
c with full missing data checking.
|
|
|
1146 |
c
|
|
|
1147 |
c a real inp array to be filtered, dimensioned (nx,ny)
|
|
|
1148 |
c af real out filtered array, dimensioned (nx,ny), can be
|
|
|
1149 |
c equivalenced with array a in the calling routine
|
|
|
1150 |
c f1 real workarray, dimensioned (nx+1,ny)
|
|
|
1151 |
c f2 real workarray, dimensioned (nx,ny+1)
|
|
|
1152 |
c fil real inp filter-coeff., 0<afil<=1. Maximum filtering with afil=1
|
|
|
1153 |
c corresponds to one application of the linear filter.
|
|
|
1154 |
c misdat real inp missing-data value, a(i,j)=misdat indicates that
|
|
|
1155 |
c the corresponding value is not available. The
|
|
|
1156 |
c misdat-checking can be switched off with with misdat=0.
|
|
|
1157 |
c iperx int inp periodic boundaries in the x-direction (1=yes,0=no)
|
|
|
1158 |
c ipery int inp periodic boundaries in the y-direction (1=yes,0=no)
|
|
|
1159 |
c inpol int inp northpole at j=ny (1=yes,0=no)
|
|
|
1160 |
c ispol int inp southpole at j=1 (1=yes,0=no)
|
|
|
1161 |
c
|
|
|
1162 |
c Christoph Schaer, 1993
|
|
|
1163 |
|
|
|
1164 |
c argument declaration
|
|
|
1165 |
integer nx,ny
|
|
|
1166 |
real a(nx,ny),af(nx,ny),fil,misdat
|
|
|
1167 |
integer iperx,ipery,inpol,ispol
|
|
|
1168 |
|
|
|
1169 |
c local variable declaration
|
|
|
1170 |
integer i,j,is
|
|
|
1171 |
real fh
|
|
|
1172 |
real f1(nx+1,ny),f2(nx,ny+1)
|
|
|
1173 |
|
|
|
1174 |
c compute constant fh
|
|
|
1175 |
fh=0.125*fil
|
|
|
1176 |
|
|
|
1177 |
c compute fluxes in x-direction
|
|
|
1178 |
if (misdat.eq.0.) then
|
|
|
1179 |
do j=1,ny
|
|
|
1180 |
do i=2,nx
|
|
|
1181 |
f1(i,j)=a(i-1,j)-a(i,j)
|
|
|
1182 |
enddo
|
|
|
1183 |
enddo
|
|
|
1184 |
else
|
|
|
1185 |
do j=1,ny
|
|
|
1186 |
do i=2,nx
|
|
|
1187 |
if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then
|
|
|
1188 |
f1(i,j)=0.
|
|
|
1189 |
else
|
|
|
1190 |
f1(i,j)=a(i-1,j)-a(i,j)
|
|
|
1191 |
endif
|
|
|
1192 |
enddo
|
|
|
1193 |
enddo
|
|
|
1194 |
endif
|
|
|
1195 |
if (iperx.eq.1) then
|
|
|
1196 |
c do periodic boundaries in the x-direction
|
|
|
1197 |
do j=1,ny
|
|
|
1198 |
f1(1,j)=f1(nx,j)
|
|
|
1199 |
f1(nx+1,j)=f1(2,j)
|
|
|
1200 |
enddo
|
|
|
1201 |
else
|
|
|
1202 |
c set boundary-fluxes to zero
|
|
|
1203 |
do j=1,ny
|
|
|
1204 |
f1(1,j)=0.
|
|
|
1205 |
f1(nx+1,j)=0.
|
|
|
1206 |
enddo
|
|
|
1207 |
endif
|
|
|
1208 |
|
|
|
1209 |
c compute fluxes in y-direction
|
|
|
1210 |
if (misdat.eq.0.) then
|
|
|
1211 |
do j=2,ny
|
|
|
1212 |
do i=1,nx
|
|
|
1213 |
f2(i,j)=a(i,j-1)-a(i,j)
|
|
|
1214 |
enddo
|
|
|
1215 |
enddo
|
|
|
1216 |
else
|
|
|
1217 |
do j=2,ny
|
|
|
1218 |
do i=1,nx
|
|
|
1219 |
if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then
|
|
|
1220 |
f2(i,j)=0.
|
|
|
1221 |
else
|
|
|
1222 |
f2(i,j)=a(i,j-1)-a(i,j)
|
|
|
1223 |
endif
|
|
|
1224 |
enddo
|
|
|
1225 |
enddo
|
|
|
1226 |
endif
|
|
|
1227 |
c set boundary-fluxes to zero
|
|
|
1228 |
do i=1,nx
|
|
|
1229 |
f2(i,1)=0.
|
|
|
1230 |
f2(i,ny+1)=0.
|
|
|
1231 |
enddo
|
|
|
1232 |
if (ipery.eq.1) then
|
|
|
1233 |
c do periodic boundaries in the x-direction
|
|
|
1234 |
do i=1,nx
|
|
|
1235 |
f2(i,1)=f2(i,ny)
|
|
|
1236 |
f2(i,ny+1)=f2(i,2)
|
|
|
1237 |
enddo
|
|
|
1238 |
endif
|
|
|
1239 |
if (iperx.eq.1) then
|
|
|
1240 |
if (ispol.eq.1) then
|
|
|
1241 |
c do south-pole
|
|
|
1242 |
is=(nx-1)/2
|
|
|
1243 |
do i=1,nx
|
|
|
1244 |
f2(i,1)=-f2(mod(i-1+is,nx)+1,2)
|
|
|
1245 |
enddo
|
|
|
1246 |
endif
|
|
|
1247 |
if (inpol.eq.1) then
|
|
|
1248 |
c do north-pole
|
|
|
1249 |
is=(nx-1)/2
|
|
|
1250 |
do i=1,nx
|
|
|
1251 |
f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny)
|
|
|
1252 |
enddo
|
|
|
1253 |
endif
|
|
|
1254 |
endif
|
|
|
1255 |
|
|
|
1256 |
c compute flux-convergence -> filter
|
|
|
1257 |
if (misdat.eq.0.) then
|
|
|
1258 |
do j=1,ny
|
|
|
1259 |
do i=1,nx
|
|
|
1260 |
af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
|
|
|
1261 |
enddo
|
|
|
1262 |
enddo
|
|
|
1263 |
else
|
|
|
1264 |
do j=1,ny
|
|
|
1265 |
do i=1,nx
|
|
|
1266 |
if (a(i,j).eq.misdat) then
|
|
|
1267 |
af(i,j)=misdat
|
|
|
1268 |
else
|
|
|
1269 |
af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
|
|
|
1270 |
endif
|
|
|
1271 |
enddo
|
|
|
1272 |
enddo
|
|
|
1273 |
endif
|
|
|
1274 |
end
|
|
|
1275 |
|
|
|
1276 |
c --------------------------------------------------------------------------------
|
|
|
1277 |
c Read refernece profile from netcdf
|
|
|
1278 |
c --------------------------------------------------------------------------------
|
|
|
1279 |
|
|
|
1280 |
SUBROUTINE read_ref (nsqref,rhoref,thetaref,pressref,zref,
|
|
|
1281 |
> nx,ny,nz,deltax,deltay,deltaz,coriol,oro,
|
|
|
1282 |
> pvsrcfile)
|
|
|
1283 |
|
|
|
1284 |
c Read the reference profile from file
|
|
|
1285 |
c
|
|
|
1286 |
c thetaref : Reference potential temperature (K)
|
|
|
1287 |
c pressref : Reference pressure (Pa)
|
|
|
1288 |
c rhoref : Reference density (kg/m^3)
|
|
|
1289 |
c nsqref : Stratification (s^-1)
|
|
|
1290 |
c zref : Reference height (m)
|
|
|
1291 |
c nx,nny,nz : Grid dimension in x,y,z direction
|
|
|
1292 |
c deltax,deltay,deltaz : Grid spacings used for calculations (m)
|
|
|
1293 |
c coriol : Coriolis parameter (s^-1)
|
|
|
1294 |
c oro : Height of orography (m)
|
|
|
1295 |
c pvsrcfile : Input file
|
|
|
1296 |
|
|
|
1297 |
implicit none
|
|
|
1298 |
|
|
|
1299 |
c Declaration of subroutine parameters
|
|
|
1300 |
integer nx,ny,nz
|
|
|
1301 |
real nsqref (0:2*nz)
|
|
|
1302 |
real thetaref(0:2*nz)
|
|
|
1303 |
real rhoref (0:2*nz)
|
|
|
1304 |
real pressref(0:2*nz)
|
|
|
1305 |
real zref (0:2*nz)
|
|
|
1306 |
real deltax,deltay,deltaz
|
|
|
1307 |
real coriol (0:nx,0:ny)
|
|
|
1308 |
real oro (0:nx,0:ny)
|
|
|
1309 |
character*80 pvsrcfile
|
|
|
1310 |
|
|
|
1311 |
c Numerical and physical parameters
|
|
|
1312 |
real eps
|
|
|
1313 |
parameter (eps=0.01)
|
|
|
1314 |
|
|
|
1315 |
c Auxiliary variables
|
|
|
1316 |
integer cdfid,stat
|
|
|
1317 |
integer vardim(4)
|
|
|
1318 |
real misdat
|
|
|
1319 |
integer ndimin
|
|
|
1320 |
real varmin(4),varmax(4),stag(4)
|
|
|
1321 |
integer i,j,k,nf1
|
|
|
1322 |
integer ntimes
|
|
|
1323 |
real time(200)
|
|
|
1324 |
character*80 vnam(100),varname
|
|
|
1325 |
integer nvars
|
|
|
1326 |
integer isok,ierr
|
|
|
1327 |
real x(0:nx,0:ny),y(0:nx,0:ny)
|
|
|
1328 |
real mean,count
|
|
|
1329 |
|
|
|
1330 |
c Get grid description from topography
|
|
|
1331 |
call cdfopn(pvsrcfile,cdfid,stat)
|
|
|
1332 |
if (stat.ne.0) goto 997
|
|
|
1333 |
call getvars(cdfid,nvars,vnam,stat)
|
|
|
1334 |
if (stat.ne.0) goto 997
|
|
|
1335 |
isok=0
|
|
|
1336 |
varname='ORO'
|
|
|
1337 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
1338 |
if (isok.eq.0) goto 997
|
|
|
1339 |
call getdef(cdfid,varname,ndimin,misdat,vardim,
|
|
|
1340 |
> varmin,varmax,stag,stat)
|
|
|
1341 |
if (stat.ne.0) goto 997
|
|
|
1342 |
time(1)=0.
|
|
|
1343 |
call gettimes(cdfid,time,ntimes,stat)
|
|
|
1344 |
if (stat.ne.0) goto 997
|
|
|
1345 |
call clscdf(cdfid,stat)
|
|
|
1346 |
if (stat.ne.0) goto 997
|
|
|
1347 |
|
|
|
1348 |
c Open output netcdf file
|
|
|
1349 |
call cdfopn(pvsrcfile,cdfid,stat)
|
|
|
1350 |
if (stat.ne.0) goto 997
|
|
|
1351 |
|
|
|
1352 |
c Create the variable if necessary
|
|
|
1353 |
call getvars(cdfid,nvars,vnam,stat)
|
|
|
1354 |
if (stat.ne.0) goto 997
|
|
|
1355 |
|
|
|
1356 |
c Read data from netcdf file
|
|
|
1357 |
isok=0
|
|
|
1358 |
varname='NSQREF'
|
|
|
1359 |
print*,'R ',trim(varname),' ',trim(pvsrcfile)
|
|
|
1360 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
1361 |
if (isok.eq.0) goto 997
|
|
|
1362 |
call getdat(cdfid,varname,time(1),0,nsqref,stat)
|
|
|
1363 |
if (stat.ne.0) goto 997
|
|
|
1364 |
|
|
|
1365 |
isok=0
|
|
|
1366 |
varname='RHOREF'
|
|
|
1367 |
print*,'R ',trim(varname),' ',trim(pvsrcfile)
|
|
|
1368 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
1369 |
if (isok.eq.0) goto 997
|
|
|
1370 |
call getdat(cdfid,varname,time(1),0,rhoref,stat)
|
|
|
1371 |
if (stat.ne.0) goto 997
|
|
|
1372 |
|
|
|
1373 |
isok=0
|
|
|
1374 |
varname='THETAREF'
|
|
|
1375 |
print*,'R ',trim(varname),' ',trim(pvsrcfile)
|
|
|
1376 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
1377 |
if (isok.eq.0) goto 997
|
|
|
1378 |
call getdat(cdfid,varname,time(1),0,thetaref,stat)
|
|
|
1379 |
if (stat.ne.0) goto 997
|
|
|
1380 |
|
|
|
1381 |
isok=0
|
|
|
1382 |
varname='PREREF'
|
|
|
1383 |
print*,'R ',trim(varname),' ',trim(pvsrcfile)
|
|
|
1384 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
1385 |
if (isok.eq.0) goto 997
|
|
|
1386 |
call getdat(cdfid,varname,time(1),0,pressref,stat)
|
|
|
1387 |
if (stat.ne.0) goto 997
|
|
|
1388 |
|
|
|
1389 |
isok=0
|
|
|
1390 |
varname='ZREF'
|
|
|
1391 |
print*,'R ',trim(varname),' ',trim(pvsrcfile)
|
|
|
1392 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
1393 |
if (isok.eq.0) goto 997
|
|
|
1394 |
call getdat(cdfid,varname,time(1),0,zref,stat)
|
|
|
1395 |
if (stat.ne.0) goto 997
|
|
|
1396 |
|
|
|
1397 |
isok=0
|
|
|
1398 |
varname='CORIOL'
|
|
|
1399 |
print*,'R ',trim(varname),' ',trim(pvsrcfile)
|
|
|
1400 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
1401 |
if (isok.eq.0) goto 997
|
|
|
1402 |
call getdat(cdfid,varname,time(1),0,coriol,stat)
|
|
|
1403 |
if (stat.ne.0) goto 997
|
|
|
1404 |
|
|
|
1405 |
isok=0
|
|
|
1406 |
varname='ORO'
|
|
|
1407 |
print*,'R ',trim(varname),' ',trim(pvsrcfile)
|
|
|
1408 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
1409 |
if (isok.eq.0) goto 997
|
|
|
1410 |
call getdat(cdfid,varname,time(1),0,oro,stat)
|
|
|
1411 |
if (stat.ne.0) goto 997
|
|
|
1412 |
|
|
|
1413 |
isok=0
|
|
|
1414 |
varname='X'
|
|
|
1415 |
print*,'R ',trim(varname),' ',trim(pvsrcfile)
|
|
|
1416 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
1417 |
if (isok.eq.0) goto 997
|
|
|
1418 |
call getdat(cdfid,varname,time(1),0,x,stat)
|
|
|
1419 |
if (stat.ne.0) goto 997
|
|
|
1420 |
|
|
|
1421 |
isok=0
|
|
|
1422 |
varname='Y'
|
|
|
1423 |
print*,'R ',trim(varname),' ',trim(pvsrcfile)
|
|
|
1424 |
call check_varok(isok,varname,vnam,nvars)
|
|
|
1425 |
if (isok.eq.0) goto 997
|
|
|
1426 |
call getdat(cdfid,varname,time(1),0,y,stat)
|
|
|
1427 |
if (stat.ne.0) goto 997
|
|
|
1428 |
|
|
|
1429 |
c Close netcdf file
|
|
|
1430 |
call clscdf(cdfid,stat)
|
|
|
1431 |
if (stat.ne.0) goto 997
|
|
|
1432 |
|
|
|
1433 |
c Determine the grid spacings <deltax, deltay, deltaz>
|
|
|
1434 |
mean=0.
|
|
|
1435 |
count=0.
|
|
|
1436 |
do i=1,nx
|
|
|
1437 |
do j=0,ny
|
|
|
1438 |
mean=mean+abs(x(i)-x(i-1))
|
|
|
1439 |
count=count+1.
|
|
|
1440 |
enddo
|
|
|
1441 |
enddo
|
|
|
1442 |
deltax=mean/count
|
|
|
1443 |
|
|
|
1444 |
mean=0.
|
|
|
1445 |
count=0.
|
|
|
1446 |
do j=1,ny
|
|
|
1447 |
do i=0,nx
|
|
|
1448 |
mean=mean+abs(y(j)-y(j-1))
|
|
|
1449 |
count=count+1.
|
|
|
1450 |
enddo
|
|
|
1451 |
enddo
|
|
|
1452 |
deltay=mean/count
|
|
|
1453 |
|
|
|
1454 |
mean=0.
|
|
|
1455 |
count=0.
|
|
|
1456 |
do k=1,nz-1
|
|
|
1457 |
mean=mean+abs(zref(k+1)-zref(k-1))
|
|
|
1458 |
count=count+1.
|
|
|
1459 |
enddo
|
|
|
1460 |
deltaz=mean/count
|
|
|
1461 |
|
|
|
1462 |
return
|
|
|
1463 |
|
|
|
1464 |
c Exception handling
|
|
|
1465 |
997 print*,'Read_Ref: Problem with input netcdf file... Stop'
|
|
|
1466 |
stop
|
|
|
1467 |
|
|
|
1468 |
end
|