| 3 |
michaesp |
1 |
PROGRAM add2p
|
|
|
2 |
|
|
|
3 |
c *********************************************************************
|
|
|
4 |
c * Add output from the PV inversion to the original P file *
|
|
|
5 |
c * Zonal wind (U), meridional wind (V) and temperature (T) are *
|
|
|
6 |
c * modified. All other fields are not affected by the programme *
|
|
|
7 |
c * *
|
|
|
8 |
c * Michael Sprenger / Summer, Autumn 2006
|
|
|
9 |
c *********************************************************************
|
|
|
10 |
|
|
|
11 |
implicit none
|
|
|
12 |
|
|
|
13 |
c ---------------------------------------------------------------------
|
|
|
14 |
c Declaration of parameters and variables
|
|
|
15 |
c ---------------------------------------------------------------------
|
|
|
16 |
|
|
|
17 |
c Input parameters
|
|
|
18 |
character*80 ml_fn
|
|
|
19 |
character*80 zm_fn
|
|
|
20 |
character*80 or_fn
|
|
|
21 |
real transxy
|
|
|
22 |
|
|
|
23 |
c Variables for input P file : model level
|
|
|
24 |
real ml_varmin(4),ml_varmax(4),ml_stag(4)
|
|
|
25 |
integer ml_vardim(4)
|
|
|
26 |
real ml_mdv
|
|
|
27 |
integer ml_ndim
|
|
|
28 |
integer ml_nx,ml_ny,ml_nz
|
|
|
29 |
real ml_xmin,ml_xmax,ml_ymin,ml_ymax,ml_dx,ml_dy
|
|
|
30 |
integer ml_ntimes
|
|
|
31 |
real ml_aklev(500),ml_bklev(500)
|
|
|
32 |
real ml_aklay(500),ml_bklay(500)
|
|
|
33 |
real ml_time
|
|
|
34 |
real ml_pollon,ml_pollat
|
|
|
35 |
integer ml_idate(5)
|
|
|
36 |
integer ml_nvars
|
|
|
37 |
character*80 ml_vnam(80)
|
|
|
38 |
real,allocatable, dimension (:,:) :: ml_ps,ml_oro
|
|
|
39 |
real,allocatable, dimension (:,:,:) :: ml_t3,ml_u3,ml_v3,ml_z3
|
|
|
40 |
real,allocatable, dimension (:,:,:) :: ml_p3,ml_tv3,ml_q3,ml_p0
|
|
|
41 |
integer,allocatable, dimension (:,:,:) :: flag_ml
|
|
|
42 |
|
|
|
43 |
c Variables for GEO.MOD
|
|
|
44 |
real zm_varmin(4),zm_varmax(4),zm_stag(4)
|
|
|
45 |
integer zm_vardim(4)
|
|
|
46 |
real zm_mdv
|
|
|
47 |
integer zm_ndim
|
|
|
48 |
integer zm_nx,zm_ny,zm_nz
|
|
|
49 |
real zm_xmin,zm_xmax,zm_ymin,zm_ymax,zm_dx,zm_dy
|
|
|
50 |
integer zm_ntimes
|
|
|
51 |
real zm_aklev(500),zm_bklev(500)
|
|
|
52 |
real zm_aklay(500),zm_bklay(500)
|
|
|
53 |
real zm_time
|
|
|
54 |
real zm_pollon,zm_pollat
|
|
|
55 |
integer zm_idate(5)
|
|
|
56 |
integer zm_nvars
|
|
|
57 |
character*80 zm_vnam(80)
|
|
|
58 |
real,allocatable, dimension (:,:,:) :: zm_u3,zm_v3,zm_t3,zm_z3
|
|
|
59 |
real,allocatable, dimension (:,:,:) :: zm_p3
|
|
|
60 |
real,allocatable, dimension (:,:) :: zm_rlat,zm_rlon
|
|
|
61 |
integer,allocatable, dimension (:,:,:) :: flag_zm
|
|
|
62 |
|
|
|
63 |
c Variables for GEO.ORG
|
|
|
64 |
real or_varmin(4),or_varmax(4),or_stag(4)
|
|
|
65 |
integer or_vardim(4)
|
|
|
66 |
real or_mdv
|
|
|
67 |
integer or_ndim
|
|
|
68 |
integer or_nx,or_ny,or_nz
|
|
|
69 |
real or_xmin,or_xmax,or_ymin,or_ymax,or_dx,or_dy
|
|
|
70 |
integer or_ntimes
|
|
|
71 |
real or_time
|
|
|
72 |
real or_pollon,or_pollat
|
|
|
73 |
integer or_idate(5)
|
|
|
74 |
integer or_nvars
|
|
|
75 |
character*80 or_vnam(80)
|
|
|
76 |
real,allocatable, dimension (:,:,:) :: or_p3
|
|
|
77 |
real,allocatable, dimension (:,:,:) :: or_u3,or_v3,or_t3
|
|
|
78 |
|
|
|
79 |
c Anomalies
|
|
|
80 |
real,allocatable, dimension (:,:,:) :: an_p3,an_u3,an_v3,an_t3
|
|
|
81 |
|
|
|
82 |
c Array with the weighting function
|
|
|
83 |
real,allocatable, dimension (:,:) :: weight,dist
|
|
|
84 |
|
|
|
85 |
c Flag for test mode
|
|
|
86 |
integer test
|
|
|
87 |
parameter (test=1)
|
|
|
88 |
|
|
|
89 |
c Flag for Poisson filling
|
|
|
90 |
integer poisson
|
|
|
91 |
parameter (poisson=1)
|
|
|
92 |
|
|
|
93 |
c Physical and numerical parameters
|
|
|
94 |
real eps ! Numerical epsilon
|
|
|
95 |
parameter (eps=0.001)
|
|
|
96 |
real dpextra ! Allowed range for extrapolation (% of pressure)
|
|
|
97 |
parameter (dpextra=0.1)
|
|
|
98 |
real g ! Earth's gravity
|
|
|
99 |
parameter (g=9.80665)
|
|
|
100 |
real tzero ! Conversion Celius/Kelvin
|
|
|
101 |
parameter (tzero=273.15)
|
|
|
102 |
real kappa ! Kappa (Definition of potential temperature)
|
|
|
103 |
parameter (kappa=0.6078)
|
|
|
104 |
real zerodiv ! Zero division handling
|
|
|
105 |
parameter (zerodiv=0.0000000001)
|
|
|
106 |
real dpmin ! Pressure step for hydrostatic equation
|
|
|
107 |
parameter (dpmin=10.)
|
|
|
108 |
real rdg ! Ratio gas constant / Earth's gravity
|
|
|
109 |
parameter (rdg=29.271)
|
|
|
110 |
|
|
|
111 |
c Auxiliray variables
|
|
|
112 |
integer i,j,k,l
|
|
|
113 |
character*80 varname
|
|
|
114 |
integer isok
|
|
|
115 |
integer stat
|
|
|
116 |
integer cdfid,ierr,cstid
|
|
|
117 |
character*80 cfn
|
|
|
118 |
real p1(1000),z1(1000),f1(1000),tv1(1000)
|
|
|
119 |
real spline_f1(1000),spline_tv1(1000)
|
|
|
120 |
integer n1
|
|
|
121 |
real hh
|
|
|
122 |
real pmin,pmax
|
|
|
123 |
real boundlat(100000),boundlon(1000000)
|
|
|
124 |
integer nbound
|
|
|
125 |
integer il,ir,ju,jd,im,jm
|
|
|
126 |
real lat,lon
|
|
|
127 |
real distmin,distpos
|
|
|
128 |
character*80 name
|
|
|
129 |
real pu,po,zu,zo,p,z,dp,p0,tvu,tvo,ff
|
|
|
130 |
integer lmin,n
|
|
|
131 |
real mean,nmean
|
|
|
132 |
real psfc
|
|
|
133 |
integer count0,count1
|
|
|
134 |
|
|
|
135 |
c Externals
|
|
|
136 |
real sdis,kink
|
|
|
137 |
external sdis,kink
|
|
|
138 |
|
|
|
139 |
|
|
|
140 |
c -----------------------------------------------------------------
|
|
|
141 |
c Read input fields
|
|
|
142 |
c -----------------------------------------------------------------
|
|
|
143 |
|
|
|
144 |
print*,'*********************************************************'
|
|
|
145 |
print*,'* add2p *'
|
|
|
146 |
print*,'*********************************************************'
|
|
|
147 |
print*
|
|
|
148 |
|
|
|
149 |
c Read in the parameter file
|
|
|
150 |
open(10,file='fort.10')
|
|
|
151 |
read(10,*) ml_fn
|
|
|
152 |
read(10,*) zm_fn
|
|
|
153 |
read(10,*) or_fn
|
|
|
154 |
read(10,*) name,transxy
|
|
|
155 |
if ( trim(name).ne.'TRANS_XY' ) stop
|
|
|
156 |
read(10,*) name,psfc
|
|
|
157 |
if ( trim(name).ne.'PS_CHANGE' ) stop
|
|
|
158 |
close(10)
|
|
|
159 |
print*,trim(ml_fn)
|
|
|
160 |
print*,trim(zm_fn)
|
|
|
161 |
print*,trim(or_fn)
|
|
|
162 |
print*
|
|
|
163 |
|
|
|
164 |
c Get grid description for P file : model level
|
|
|
165 |
call cdfopn(ml_fn,cdfid,ierr)
|
|
|
166 |
if (ierr.ne.0) goto 998
|
|
|
167 |
call getvars(cdfid,ml_nvars,ml_vnam,ierr)
|
|
|
168 |
if (ierr.ne.0) goto 998
|
|
|
169 |
call getcfn(cdfid,cfn,ierr)
|
|
|
170 |
if (ierr.ne.0) goto 998
|
|
|
171 |
call cdfopn(cfn,cstid,ierr)
|
|
|
172 |
if (ierr.ne.0) goto 998
|
|
|
173 |
isok=0
|
|
|
174 |
varname='T'
|
|
|
175 |
call check_varok(isok,varname,ml_vnam,ml_nvars)
|
|
|
176 |
if (isok.ne.1) goto 998
|
|
|
177 |
call getdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
|
|
|
178 |
> ml_varmin,ml_varmax,ml_stag,ierr)
|
|
|
179 |
if (ierr.ne.0) goto 998
|
|
|
180 |
ml_nx =ml_vardim(1)
|
|
|
181 |
ml_ny =ml_vardim(2)
|
|
|
182 |
ml_nz =ml_vardim(3)
|
|
|
183 |
ml_xmin=ml_varmin(1)
|
|
|
184 |
ml_ymin=ml_varmin(2)
|
|
|
185 |
call getlevs(cstid,ml_nz,ml_aklev,ml_bklev,ml_aklay,ml_bklay,ierr)
|
|
|
186 |
if (ierr.ne.0) goto 998
|
|
|
187 |
call getgrid(cstid,ml_dx,ml_dy,ierr)
|
|
|
188 |
ml_xmax=ml_xmin+real(ml_nx-1)*ml_dx
|
|
|
189 |
ml_ymax=ml_ymin+real(ml_ny-1)*ml_dy
|
|
|
190 |
if (ierr.ne.0) goto 998
|
|
|
191 |
call gettimes(cdfid,ml_time,ml_ntimes,ierr)
|
|
|
192 |
if (ierr.ne.0) goto 998
|
|
|
193 |
call getstart(cstid,ml_idate,ierr)
|
|
|
194 |
if (ierr.ne.0) goto 998
|
|
|
195 |
call getpole(cstid,ml_pollon,ml_pollat,ierr)
|
|
|
196 |
if (ierr.ne.0) goto 998
|
|
|
197 |
call clscdf(cstid,ierr)
|
|
|
198 |
if (ierr.ne.0) goto 998
|
|
|
199 |
call clscdf(cdfid,ierr)
|
|
|
200 |
if (ierr.ne.0) goto 998
|
|
|
201 |
|
|
|
202 |
c Get grid description for MOD file
|
|
|
203 |
call cdfopn(zm_fn,cdfid,ierr)
|
|
|
204 |
if (ierr.ne.0) goto 997
|
|
|
205 |
call getvars(cdfid,zm_nvars,zm_vnam,ierr)
|
|
|
206 |
if (ierr.ne.0) goto 997
|
|
|
207 |
call getcfn(cdfid,cfn,ierr)
|
|
|
208 |
if (ierr.ne.0) goto 997
|
|
|
209 |
call cdfopn(cfn,cstid,ierr)
|
|
|
210 |
if (ierr.ne.0) goto 997
|
|
|
211 |
isok=0
|
|
|
212 |
varname='T'
|
|
|
213 |
call check_varok(isok,varname,zm_vnam,zm_nvars)
|
|
|
214 |
if (isok.ne.1) goto 997
|
|
|
215 |
call getdef(cdfid,varname,zm_ndim,zm_mdv,zm_vardim,
|
|
|
216 |
> zm_varmin,zm_varmax,zm_stag,ierr)
|
|
|
217 |
if (ierr.ne.0) goto 997
|
|
|
218 |
zm_nx =zm_vardim(1)
|
|
|
219 |
zm_ny =zm_vardim(2)
|
|
|
220 |
zm_nz =zm_vardim(3)
|
|
|
221 |
zm_xmin=zm_varmin(1)
|
|
|
222 |
zm_ymin=zm_varmin(2)
|
|
|
223 |
call getlevs(cstid,zm_nz,zm_aklev,zm_bklev,zm_aklay,zm_bklay,ierr)
|
|
|
224 |
if (ierr.ne.0) goto 997
|
|
|
225 |
call getgrid(cstid,zm_dx,zm_dy,ierr)
|
|
|
226 |
if (ierr.ne.0) goto 997
|
|
|
227 |
zm_xmax=zm_xmin+real(zm_nx-1)*zm_dx
|
|
|
228 |
zm_ymax=zm_ymin+real(zm_ny-1)*zm_dy
|
|
|
229 |
call gettimes(cdfid,zm_time,zm_ntimes,ierr)
|
|
|
230 |
if (ierr.ne.0) goto 997
|
|
|
231 |
call getstart(cstid,zm_idate,ierr)
|
|
|
232 |
if (ierr.ne.0) goto 997
|
|
|
233 |
call getpole(cstid,zm_pollon,zm_pollat,ierr)
|
|
|
234 |
if (ierr.ne.0) goto 997
|
|
|
235 |
call clscdf(cstid,ierr)
|
|
|
236 |
if (ierr.ne.0) goto 997
|
|
|
237 |
call clscdf(cdfid,ierr)
|
|
|
238 |
if (ierr.ne.0) goto 997
|
|
|
239 |
|
|
|
240 |
c Get grid description for ORG file
|
|
|
241 |
call cdfopn(or_fn,cdfid,ierr)
|
|
|
242 |
if (ierr.ne.0) goto 997
|
|
|
243 |
call getvars(cdfid,or_nvars,or_vnam,ierr)
|
|
|
244 |
if (ierr.ne.0) goto 997
|
|
|
245 |
call getcfn(cdfid,cfn,ierr)
|
|
|
246 |
if (ierr.ne.0) goto 997
|
|
|
247 |
call cdfopn(cfn,cstid,ierr)
|
|
|
248 |
if (ierr.ne.0) goto 997
|
|
|
249 |
isok=0
|
|
|
250 |
varname='P'
|
|
|
251 |
call check_varok(isok,varname,or_vnam,or_nvars)
|
|
|
252 |
if (isok.ne.1) goto 997
|
|
|
253 |
call getdef(cdfid,varname,or_ndim,or_mdv,or_vardim,
|
|
|
254 |
> or_varmin,or_varmax,or_stag,ierr)
|
|
|
255 |
if (ierr.ne.0) goto 997
|
|
|
256 |
or_nx =or_vardim(1)
|
|
|
257 |
or_ny =or_vardim(2)
|
|
|
258 |
or_nz =or_vardim(3)
|
|
|
259 |
or_xmin=or_varmin(1)
|
|
|
260 |
or_ymin=or_varmin(2)
|
|
|
261 |
call getgrid(cstid,or_dx,or_dy,ierr)
|
|
|
262 |
if (ierr.ne.0) goto 997
|
|
|
263 |
or_xmax=or_xmin+real(or_nx-1)*or_dx
|
|
|
264 |
or_ymax=or_ymin+real(or_ny-1)*or_dy
|
|
|
265 |
call gettimes(cdfid,or_time,or_ntimes,ierr)
|
|
|
266 |
if (ierr.ne.0) goto 997
|
|
|
267 |
call getstart(cstid,or_idate,ierr)
|
|
|
268 |
if (ierr.ne.0) goto 997
|
|
|
269 |
call getpole(cstid,or_pollon,or_pollat,ierr)
|
|
|
270 |
if (ierr.ne.0) goto 997
|
|
|
271 |
call clscdf(cstid,ierr)
|
|
|
272 |
if (ierr.ne.0) goto 997
|
|
|
273 |
call clscdf(cdfid,ierr)
|
|
|
274 |
if (ierr.ne.0) goto 997
|
|
|
275 |
|
|
|
276 |
c Consistency check for the grids
|
|
|
277 |
if ( (abs(ml_xmin-zm_xmin).gt.eps).or.
|
|
|
278 |
> (abs(ml_xmax-zm_xmax).gt.eps).or.
|
|
|
279 |
> (abs(ml_ymin-zm_ymin).gt.eps).or.
|
|
|
280 |
> (abs(ml_ymax-zm_ymax).gt.eps).or.
|
|
|
281 |
> (abs(ml_dx -zm_dx ).gt.eps).or.
|
|
|
282 |
> (abs(ml_dy -zm_dy ).gt.eps).or.
|
|
|
283 |
> (abs(ml_time-zm_time).gt.eps).or.
|
|
|
284 |
> (abs(ml_xmin-or_xmin).gt.eps).or.
|
|
|
285 |
> (abs(ml_xmax-or_xmax).gt.eps).or.
|
|
|
286 |
> (abs(ml_ymin-or_ymin).gt.eps).or.
|
|
|
287 |
> (abs(ml_ymax-or_ymax).gt.eps).or.
|
|
|
288 |
> (abs(ml_dx -or_dx ).gt.eps).or.
|
|
|
289 |
> (abs(ml_dy -or_dy ).gt.eps) )
|
|
|
290 |
>then
|
|
|
291 |
print*,'P, GEO.MOD, GEO.ORG grids not consistent... Stop'
|
|
|
292 |
print*,ml_xmin,zm_xmin,or_xmin
|
|
|
293 |
print*,ml_ymin,zm_ymin,or_ymin
|
|
|
294 |
print*,ml_dx, zm_dx ,or_dx
|
|
|
295 |
print*,ml_dy, zm_dy ,or_dy
|
|
|
296 |
print*,ml_time,zm_time,or_time
|
|
|
297 |
stop
|
|
|
298 |
endif
|
|
|
299 |
|
|
|
300 |
c Allocate memory for input and output P file
|
|
|
301 |
allocate(ml_p3 (ml_nx,ml_ny,ml_nz),stat=stat)
|
|
|
302 |
if (stat.ne.0) print*,'error allocating ml_p3'
|
|
|
303 |
allocate(ml_u3 (ml_nx,ml_ny,ml_nz),stat=stat)
|
|
|
304 |
if (stat.ne.0) print*,'error allocating ml_u3'
|
|
|
305 |
allocate(ml_v3 (ml_nx,ml_ny,ml_nz),stat=stat)
|
|
|
306 |
if (stat.ne.0) print*,'error allocating ml_v3'
|
|
|
307 |
allocate(ml_t3 (ml_nx,ml_ny,ml_nz),stat=stat)
|
|
|
308 |
if (stat.ne.0) print*,'error allocating ml_t3'
|
|
|
309 |
allocate(ml_tv3 (ml_nx,ml_ny,ml_nz),stat=stat)
|
|
|
310 |
if (stat.ne.0) print*,'error allocating ml_tv3'
|
|
|
311 |
allocate(ml_z3 (ml_nx,ml_ny,ml_nz),stat=stat)
|
|
|
312 |
if (stat.ne.0) print*,'error allocating ml_z3'
|
|
|
313 |
allocate(ml_q3 (ml_nx,ml_ny,ml_nz),stat=stat)
|
|
|
314 |
if (stat.ne.0) print*,'error allocating ml_q3'
|
|
|
315 |
allocate(ml_ps (ml_nx,ml_ny ),stat=stat)
|
|
|
316 |
if (stat.ne.0) print*,'error allocating ml_ps'
|
|
|
317 |
allocate(ml_oro (ml_nx,ml_ny ),stat=stat)
|
|
|
318 |
if (stat.ne.0) print*,'error allocating ml_oro'
|
|
|
319 |
allocate(ml_p0 (ml_nx,ml_ny,ml_nz),stat=stat)
|
|
|
320 |
if (stat.ne.0) print*,'error allocating ml_p0'
|
|
|
321 |
|
|
|
322 |
c Allocate memory for input GEO.MOD file
|
|
|
323 |
allocate(zm_p3(zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
324 |
if (stat.ne.0) print*,'error allocating zm_p3'
|
|
|
325 |
allocate(zm_u3(zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
326 |
if (stat.ne.0) print*,'error allocating zm_u3'
|
|
|
327 |
allocate(zm_v3(zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
328 |
if (stat.ne.0) print*,'error allocating zm_v3'
|
|
|
329 |
allocate(zm_t3(zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
330 |
if (stat.ne.0) print*,'error allocating zm_t3'
|
|
|
331 |
allocate(zm_z3(zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
332 |
if (stat.ne.0) print*,'error allocating zm_z3'
|
|
|
333 |
allocate(zm_rlat(zm_nx,zm_ny),stat=stat)
|
|
|
334 |
if (stat.ne.0) print*,'error allocating zm_rlat'
|
|
|
335 |
allocate(zm_rlon(zm_nx,zm_ny),stat=stat)
|
|
|
336 |
if (stat.ne.0) print*,'error allocating zm_rlon'
|
|
|
337 |
|
|
|
338 |
c Allocate memory for input GEO.ORG file
|
|
|
339 |
allocate(or_p3(zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
340 |
if (stat.ne.0) print*,'error allocating zm_p3'
|
|
|
341 |
allocate(or_t3(zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
342 |
if (stat.ne.0) print*,'error allocating zm_t3'
|
|
|
343 |
allocate(or_u3(zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
344 |
if (stat.ne.0) print*,'error allocating zm_u3'
|
|
|
345 |
allocate(or_v3(zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
346 |
if (stat.ne.0) print*,'error allocating zm_v3'
|
|
|
347 |
|
|
|
348 |
c Allocate memory for weighting function, flag and anomalies
|
|
|
349 |
allocate(weight(zm_nx,zm_ny),stat=stat)
|
|
|
350 |
if (stat.ne.0) print*,'error allocating weight'
|
|
|
351 |
allocate(dist (zm_nx,zm_ny),stat=stat)
|
|
|
352 |
if (stat.ne.0) print*,'error allocating dist'
|
|
|
353 |
allocate(flag_zm(zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
354 |
if (stat.ne.0) print*,'error allocating flag'
|
|
|
355 |
allocate(flag_ml(zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
356 |
if (stat.ne.0) print*,'error allocating flag'
|
|
|
357 |
allocate(an_p3 (zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
358 |
if (stat.ne.0) print*,'error allocating an_p3'
|
|
|
359 |
allocate(an_t3 (zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
360 |
if (stat.ne.0) print*,'error allocating an_t3'
|
|
|
361 |
allocate(an_u3 (zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
362 |
if (stat.ne.0) print*,'error allocating an_u3'
|
|
|
363 |
allocate(an_v3 (zm_nx,zm_ny,zm_nz),stat=stat)
|
|
|
364 |
if (stat.ne.0) print*,'error allocating an_v3'
|
|
|
365 |
|
|
|
366 |
c Read variables from P file : model levels
|
|
|
367 |
call cdfopn(ml_fn,cdfid,ierr)
|
|
|
368 |
if (ierr.ne.0) goto 998
|
|
|
369 |
varname='T'
|
|
|
370 |
call getdat(cdfid,varname,ml_time,0,ml_t3,ierr)
|
|
|
371 |
if (ierr.ne.0) goto 998
|
|
|
372 |
print*,'R T ',trim(ml_fn)
|
|
|
373 |
varname='U'
|
|
|
374 |
call getdat(cdfid,varname,ml_time,0,ml_u3,ierr)
|
|
|
375 |
if (ierr.ne.0) goto 998
|
|
|
376 |
print*,'R U ',trim(ml_fn)
|
|
|
377 |
varname='V'
|
|
|
378 |
call getdat(cdfid,varname,ml_time,0,ml_v3,ierr)
|
|
|
379 |
if (ierr.ne.0) goto 998
|
|
|
380 |
print*,'R V ',trim(ml_fn)
|
|
|
381 |
varname='Z'
|
|
|
382 |
call getdat(cdfid,varname,ml_time,0,ml_z3,ierr)
|
|
|
383 |
if (ierr.ne.0) goto 998
|
|
|
384 |
print*,'R Z ',trim(ml_fn)
|
|
|
385 |
varname='Q'
|
|
|
386 |
call getdat(cdfid,varname,ml_time,0,ml_q3,ierr)
|
|
|
387 |
if (ierr.ne.0) goto 998
|
|
|
388 |
print*,'R Q ',trim(ml_fn)
|
|
|
389 |
varname='PS'
|
|
|
390 |
call getdat(cdfid,varname,ml_time,0,ml_ps,ierr)
|
|
|
391 |
if (ierr.ne.0) goto 998
|
|
|
392 |
print*,'R PS ',trim(ml_fn)
|
|
|
393 |
varname='ORO'
|
|
|
394 |
call getdat(cdfid,varname,ml_time,0,ml_oro,ierr)
|
|
|
395 |
if (ierr.ne.0) goto 998
|
|
|
396 |
print*,'R ORO ',trim(ml_fn)
|
|
|
397 |
call clscdf(cdfid,ierr)
|
|
|
398 |
if (ierr.ne.0) goto 998
|
|
|
399 |
|
|
|
400 |
c Read variables from GEO.MOD
|
|
|
401 |
call cdfopn(zm_fn,cdfid,ierr)
|
|
|
402 |
if (ierr.ne.0) goto 997
|
|
|
403 |
|
|
|
404 |
varname='T'
|
|
|
405 |
call getdat(cdfid,varname,zm_time,0,zm_t3,ierr)
|
|
|
406 |
if (ierr.ne.0) goto 997
|
|
|
407 |
print*,'R T ',trim(zm_fn)
|
|
|
408 |
varname='U'
|
|
|
409 |
call getdat(cdfid,varname,zm_time,0,zm_u3,ierr)
|
|
|
410 |
if (ierr.ne.0) goto 997
|
|
|
411 |
print*,'R U ',trim(zm_fn)
|
|
|
412 |
varname='V'
|
|
|
413 |
call getdat(cdfid,varname,zm_time,0,zm_v3,ierr)
|
|
|
414 |
if (ierr.ne.0) goto 997
|
|
|
415 |
print*,'R V ',trim(zm_fn)
|
|
|
416 |
varname='P'
|
|
|
417 |
call getdat(cdfid,varname,zm_time,0,zm_p3,ierr)
|
|
|
418 |
if (ierr.ne.0) goto 997
|
|
|
419 |
print*,'R P ',trim(zm_fn)
|
|
|
420 |
varname='RLAT'
|
|
|
421 |
call getdat(cdfid,varname,zm_time,0,zm_rlat,ierr)
|
|
|
422 |
if (ierr.ne.0) goto 997
|
|
|
423 |
print*,'R RLAT ',trim(zm_fn)
|
|
|
424 |
varname='RLON'
|
|
|
425 |
call getdat(cdfid,varname,zm_time,0,zm_rlon,ierr)
|
|
|
426 |
if (ierr.ne.0) goto 997
|
|
|
427 |
print*,'R RLON ',trim(zm_fn)
|
|
|
428 |
|
|
|
429 |
call clscdf(cdfid,ierr)
|
|
|
430 |
if (ierr.ne.0) goto 997
|
|
|
431 |
|
|
|
432 |
c Read variables from GEO.ORG
|
|
|
433 |
call cdfopn(or_fn,cdfid,ierr)
|
|
|
434 |
if (ierr.ne.0) goto 998
|
|
|
435 |
varname='P'
|
|
|
436 |
call getdat(cdfid,varname,or_time,0,or_p3,ierr)
|
|
|
437 |
if (ierr.ne.0) goto 998
|
|
|
438 |
print*,'R P ',trim(or_fn)
|
|
|
439 |
varname='U'
|
|
|
440 |
call getdat(cdfid,varname,or_time,0,or_u3,ierr)
|
|
|
441 |
if (ierr.ne.0) goto 998
|
|
|
442 |
print*,'R U ',trim(or_fn)
|
|
|
443 |
varname='V'
|
|
|
444 |
call getdat(cdfid,varname,or_time,0,or_v3,ierr)
|
|
|
445 |
if (ierr.ne.0) goto 998
|
|
|
446 |
print*,'R V ',trim(or_fn)
|
|
|
447 |
varname='T'
|
|
|
448 |
call getdat(cdfid,varname,or_time,0,or_t3,ierr)
|
|
|
449 |
if (ierr.ne.0) goto 998
|
|
|
450 |
print*,'R T ',trim(or_fn)
|
|
|
451 |
call clscdf(cdfid,ierr)
|
|
|
452 |
if (ierr.ne.0) goto 998
|
|
|
453 |
|
|
|
454 |
c Initialize the height levels of the Z file
|
|
|
455 |
do i=1,zm_nx
|
|
|
456 |
do j=1,zm_ny
|
|
|
457 |
do k=1,zm_nz
|
|
|
458 |
zm_z3(i,j,k)=zm_aklay(k)
|
|
|
459 |
enddo
|
|
|
460 |
enddo
|
|
|
461 |
enddo
|
|
|
462 |
|
|
|
463 |
c Calculate 3d pressure field on model levels
|
|
|
464 |
print*,'C P (ORIGINAL)'
|
|
|
465 |
do k=1,ml_nz
|
|
|
466 |
do i=1,ml_nx
|
|
|
467 |
do j=1,ml_ny
|
|
|
468 |
if (abs(ml_stag(3)+0.5).lt.eps) then
|
|
|
469 |
ml_p0(i,j,k)=ml_aklay(k)+ml_bklay(k)*ml_ps(i,j)
|
|
|
470 |
else
|
|
|
471 |
ml_p0(i,j,k)=ml_aklev(k)+ml_bklev(k)*ml_ps(i,j)
|
|
|
472 |
endif
|
|
|
473 |
enddo
|
|
|
474 |
enddo
|
|
|
475 |
enddo
|
|
|
476 |
|
|
|
477 |
c Calculate 3d anomalies due to inversion
|
|
|
478 |
print*,'C DP (MODIFIED-ORIGINAL)'
|
|
|
479 |
print*,'C DU (MODIFIED-ORIGINAL)'
|
|
|
480 |
print*,'C DV (MODIFIED-ORIGINAL)'
|
|
|
481 |
print*,'C DT (MODIFIED-ORIGINAL)'
|
|
|
482 |
do k=1,zm_nz
|
|
|
483 |
do i=1,zm_nx
|
|
|
484 |
do j=1,zm_ny
|
|
|
485 |
|
|
|
486 |
c P
|
|
|
487 |
if ( ( abs(zm_p3(i,j,k)-zm_mdv).gt.eps).and.
|
|
|
488 |
> ( abs(or_p3(i,j,k)-or_mdv).gt.eps) )
|
|
|
489 |
> then
|
|
|
490 |
an_p3(i,j,k) = zm_p3(i,j,k) - or_p3(i,j,k)
|
|
|
491 |
elseif ( poisson.eq.1 ) then
|
|
|
492 |
an_p3(i,j,k) = zm_mdv
|
|
|
493 |
else
|
|
|
494 |
an_p3(i,j,k) = 0.
|
|
|
495 |
endif
|
|
|
496 |
|
|
|
497 |
c T
|
|
|
498 |
if ( ( abs(zm_t3(i,j,k)-zm_mdv).gt.eps).and.
|
|
|
499 |
> ( abs(or_t3(i,j,k)-or_mdv).gt.eps) )
|
|
|
500 |
> then
|
|
|
501 |
an_t3(i,j,k) = zm_t3(i,j,k) - or_t3(i,j,k)
|
|
|
502 |
elseif ( poisson.eq.1 ) then
|
|
|
503 |
an_t3(i,j,k) = zm_mdv
|
|
|
504 |
else
|
|
|
505 |
an_t3(i,j,k) = 0.
|
|
|
506 |
endif
|
|
|
507 |
|
|
|
508 |
c U
|
|
|
509 |
if ( ( abs(zm_u3(i,j,k)-zm_mdv).gt.eps).and.
|
|
|
510 |
> ( abs(or_u3(i,j,k)-or_mdv).gt.eps) )
|
|
|
511 |
> then
|
|
|
512 |
an_u3(i,j,k) = zm_u3(i,j,k) - or_u3(i,j,k)
|
|
|
513 |
elseif ( poisson.eq.1 ) then
|
|
|
514 |
an_u3(i,j,k) = zm_mdv
|
|
|
515 |
else
|
|
|
516 |
an_u3(i,j,k) = 0.
|
|
|
517 |
endif
|
|
|
518 |
|
|
|
519 |
c V
|
|
|
520 |
if ( ( abs(zm_v3(i,j,k)-zm_mdv).gt.eps).and.
|
|
|
521 |
> ( abs(or_v3(i,j,k)-or_mdv).gt.eps) )
|
|
|
522 |
> then
|
|
|
523 |
an_v3(i,j,k) = zm_v3(i,j,k) - or_v3(i,j,k)
|
|
|
524 |
elseif ( poisson.eq.1 ) then
|
|
|
525 |
an_v3(i,j,k) = zm_mdv
|
|
|
526 |
else
|
|
|
527 |
an_v3(i,j,k) = 0.
|
|
|
528 |
endif
|
|
|
529 |
|
|
|
530 |
enddo
|
|
|
531 |
enddo
|
|
|
532 |
enddo
|
|
|
533 |
|
|
|
534 |
c ----------------------------------------------------------------
|
|
|
535 |
c Get the weight function for the inset (from 1 inside to 0 outside)
|
|
|
536 |
c ----------------------------------------------------------------
|
|
|
537 |
|
|
|
538 |
print*,'C BOUNDARY FILTER'
|
|
|
539 |
|
|
|
540 |
c Init the weight function (1 inside, 0 outside)
|
|
|
541 |
do i=1,ml_nx
|
|
|
542 |
do j=1,ml_ny
|
|
|
543 |
|
|
|
544 |
if ( (zm_rlat(i,j).lt. -90.).or.
|
|
|
545 |
> (zm_rlat(i,j).gt. 90.).or.
|
|
|
546 |
> (zm_rlon(i,j).lt.-180.).or.
|
|
|
547 |
> (zm_rlon(i,j).gt. 180.) ) then
|
|
|
548 |
weight(i,j)=0.
|
|
|
549 |
else
|
|
|
550 |
weight(i,j)=1.
|
|
|
551 |
endif
|
|
|
552 |
|
|
|
553 |
enddo
|
|
|
554 |
enddo
|
|
|
555 |
|
|
|
556 |
c Get a list of all boundary points
|
|
|
557 |
nbound=0
|
|
|
558 |
do i=1,ml_nx
|
|
|
559 |
do j=1,ml_ny
|
|
|
560 |
|
|
|
561 |
c Get neighbouring points
|
|
|
562 |
ir=i+1
|
|
|
563 |
if (ir.gt.ml_nx) ir=ml_nx
|
|
|
564 |
il=i-1
|
|
|
565 |
if (il.lt. 1) il=1
|
|
|
566 |
ju=j+1
|
|
|
567 |
if (ju.gt.ml_ny) ju=ml_ny
|
|
|
568 |
jd=j-1
|
|
|
569 |
if (jd.lt. 1) jd=1
|
|
|
570 |
|
|
|
571 |
c A boundary point has a 0/1 switch
|
|
|
572 |
if (abs(weight(i,j)-0.).lt.eps) then
|
|
|
573 |
|
|
|
574 |
if ( (abs(weight(ir, j)-1.).lt.eps).or.
|
|
|
575 |
> (abs(weight(il, j)-1.).lt.eps).or.
|
|
|
576 |
> (abs(weight(i ,ju)-1.).lt.eps).or.
|
|
|
577 |
> (abs(weight(i ,jd)-1.).lt.eps).or.
|
|
|
578 |
> (abs(weight(ir,ju)-1.).lt.eps).or.
|
|
|
579 |
> (abs(weight(ir,jd)-1.).lt.eps).or.
|
|
|
580 |
> (abs(weight(il,ju)-1.).lt.eps).or.
|
|
|
581 |
> (abs(weight(il,jd)-1.).lt.eps).or.
|
|
|
582 |
> (abs(weight(ir,ju)-1.).lt.eps).or.
|
|
|
583 |
> (abs(weight(il,ju)-1.).lt.eps).or.
|
|
|
584 |
> (abs(weight(ir,jd)-1.).lt.eps).or.
|
|
|
585 |
> (abs(weight(il,jd)-1.).lt.eps) ) then
|
|
|
586 |
|
|
|
587 |
nbound=nbound+1
|
|
|
588 |
boundlon(nbound)=ml_xmin+real(i-1)*ml_dx
|
|
|
589 |
boundlat(nbound)=ml_ymin+real(j-1)*ml_dy
|
|
|
590 |
|
|
|
591 |
endif
|
|
|
592 |
|
|
|
593 |
endif
|
|
|
594 |
|
|
|
595 |
enddo
|
|
|
596 |
enddo
|
|
|
597 |
|
|
|
598 |
c Get the distance from the subdomain
|
|
|
599 |
do i=1,ml_nx
|
|
|
600 |
do j=1,ml_ny
|
|
|
601 |
|
|
|
602 |
lon=ml_xmin+real(i-1)*ml_dx
|
|
|
603 |
lat=ml_ymin+real(j-1)*ml_dy
|
|
|
604 |
|
|
|
605 |
distmin=sdis(lon,lat,boundlon(1),boundlat(1))
|
|
|
606 |
do k=2,nbound
|
|
|
607 |
distpos=sdis(lon,lat,boundlon(k),boundlat(k))
|
|
|
608 |
if (distpos.lt.distmin) distmin=distpos
|
|
|
609 |
enddo
|
|
|
610 |
|
|
|
611 |
if ( abs(weight(i,j)-1.).lt.eps) then
|
|
|
612 |
dist(i,j)=distmin
|
|
|
613 |
else
|
|
|
614 |
dist(i,j)=-distmin
|
|
|
615 |
endif
|
|
|
616 |
|
|
|
617 |
enddo
|
|
|
618 |
enddo
|
|
|
619 |
|
|
|
620 |
c Set the new weights
|
|
|
621 |
do i=1,ml_nx
|
|
|
622 |
do j=1,ml_ny
|
|
|
623 |
|
|
|
624 |
if (weight(i,j).gt.0.) then
|
|
|
625 |
weight(i,j)=kink(dist(i,j),transxy)
|
|
|
626 |
endif
|
|
|
627 |
|
|
|
628 |
enddo
|
|
|
629 |
enddo
|
|
|
630 |
|
|
|
631 |
c ----------------------------------------------------------------
|
|
|
632 |
c Remove MDV regions in input field
|
|
|
633 |
c ----------------------------------------------------------------
|
|
|
634 |
|
|
|
635 |
if ( poisson.ne.1 ) goto 120
|
|
|
636 |
|
|
|
637 |
c Define region for mdv filling
|
|
|
638 |
do i=1,zm_nx
|
|
|
639 |
do j=1,zm_ny
|
|
|
640 |
do k=1,zm_nz
|
|
|
641 |
|
|
|
642 |
c Get neighbour of grid point
|
|
|
643 |
il = i-1
|
|
|
644 |
if (il.lt.1) il=1
|
|
|
645 |
ir = i+1
|
|
|
646 |
if ( ir.gt.zm_nx ) ir=zm_nx
|
|
|
647 |
jd = j-1
|
|
|
648 |
if (jd.lt.1) jd=1
|
|
|
649 |
ju = j+1
|
|
|
650 |
if ( ju.gt.zm_ny ) ju=zm_ny
|
|
|
651 |
|
|
|
652 |
c Set flag 2 for boundary and exterior points
|
|
|
653 |
if ( (abs(zm_rlat(il, j)-zm_mdv).lt.eps).or.
|
|
|
654 |
> (abs(zm_rlat(ir, j)-zm_mdv).lt.eps).or.
|
|
|
655 |
> (abs(zm_rlat(i , j)-zm_mdv).lt.eps).or.
|
|
|
656 |
> (abs(zm_rlat(il,ju)-zm_mdv).lt.eps).or.
|
|
|
657 |
> (abs(zm_rlat(ir,ju)-zm_mdv).lt.eps).or.
|
|
|
658 |
> (abs(zm_rlat(i ,ju)-zm_mdv).lt.eps).or.
|
|
|
659 |
> (abs(zm_rlat(il,jd)-zm_mdv).lt.eps).or.
|
|
|
660 |
> (abs(zm_rlat(ir,jd)-zm_mdv).lt.eps).or.
|
|
|
661 |
> (abs(zm_rlat(i ,jd)-zm_mdv).lt.eps) )
|
|
|
662 |
> then
|
|
|
663 |
flag_zm(i,j,k) = 2
|
|
|
664 |
an_p3(i,j,k) = 0.
|
|
|
665 |
an_u3(i,j,k) = 0.
|
|
|
666 |
an_v3(i,j,k) = 0.
|
|
|
667 |
an_t3(i,j,k) = 0.
|
|
|
668 |
|
|
|
669 |
c Set flag 1 for interior MDV points
|
|
|
670 |
elseif ( abs(an_p3(i,j,k)-zm_mdv).le.eps ) then
|
|
|
671 |
flag_zm(i,j,k) = 1
|
|
|
672 |
|
|
|
673 |
c Set flag 0 for interior valid points
|
|
|
674 |
else
|
|
|
675 |
flag_zm(i,j,k) = 0
|
|
|
676 |
endif
|
|
|
677 |
|
|
|
678 |
enddo
|
|
|
679 |
enddo
|
|
|
680 |
enddo
|
|
|
681 |
|
|
|
682 |
c Apply mdv filling
|
|
|
683 |
print*,'C DP (POISSON FILLING)'
|
|
|
684 |
call mdvfill(an_p3,an_p3,flag_zm,zm_nx,zm_ny,zm_nz,100)
|
|
|
685 |
print*,'C DT (POISSON FILLING)'
|
|
|
686 |
call mdvfill(an_t3,an_t3,flag_zm,zm_nx,zm_ny,zm_nz,100)
|
|
|
687 |
print*,'C DU (POISSON FILLING)'
|
|
|
688 |
call mdvfill(an_u3,an_u3,flag_zm,zm_nx,zm_ny,zm_nz,100)
|
|
|
689 |
print*,'C DV (POISSON FILLING)'
|
|
|
690 |
call mdvfill(an_v3,an_v3,flag_zm,zm_nx,zm_ny,zm_nz,100)
|
|
|
691 |
|
|
|
692 |
c Special treatment: if the number of missing values surpasses 50%
|
|
|
693 |
c on a level, then the anomaly is imported from the level above
|
|
|
694 |
do k=zm_nz-1,1,-1
|
|
|
695 |
|
|
|
696 |
count1 = 0
|
|
|
697 |
count0 = 0
|
|
|
698 |
do i=1,zm_nx
|
|
|
699 |
do j= 1,zm_ny
|
|
|
700 |
if ( flag_zm(i,j,k).eq.1 ) count1 = count1 + 1
|
|
|
701 |
if ( flag_zm(i,j,k).eq.0 ) count0 = count0 + 1
|
|
|
702 |
enddo
|
|
|
703 |
enddo
|
|
|
704 |
|
|
|
705 |
if ( count1.gt.count0 ) then
|
|
|
706 |
print*,'C P (IMPORTING FROM LEVEL ABOVE)',k
|
|
|
707 |
do i=1,zm_nx
|
|
|
708 |
do j= 1,zm_ny
|
|
|
709 |
an_p3(i,j,k) = an_p3(i,j,k+1)
|
|
|
710 |
flag_zm(i,j,k) = flag_zm(i,j,k+1)
|
|
|
711 |
enddo
|
|
|
712 |
enddo
|
|
|
713 |
print*,'C U (IMPORTING FROM LEVEL ABOVE)',k
|
|
|
714 |
do i=1,zm_nx
|
|
|
715 |
do j= 1,zm_ny
|
|
|
716 |
an_u3(i,j,k) = an_u3(i,j,k+1)
|
|
|
717 |
flag_zm(i,j,k) = flag_zm(i,j,k+1)
|
|
|
718 |
enddo
|
|
|
719 |
enddo
|
|
|
720 |
|
|
|
721 |
print*,'C V (IMPORTING FROM LEVEL ABOVE)',k
|
|
|
722 |
do i=1,zm_nx
|
|
|
723 |
do j= 1,zm_ny
|
|
|
724 |
an_v3(i,j,k) = an_v3(i,j,k+1)
|
|
|
725 |
flag_zm(i,j,k) = flag_zm(i,j,k+1)
|
|
|
726 |
enddo
|
|
|
727 |
enddo
|
|
|
728 |
|
|
|
729 |
print*,'C T (IMPORTING FROM LEVEL ABOVE)',k
|
|
|
730 |
do i=1,zm_nx
|
|
|
731 |
do j= 1,zm_ny
|
|
|
732 |
an_t3(i,j,k) = an_t3(i,j,k+1)
|
|
|
733 |
flag_zm(i,j,k) = flag_zm(i,j,k+1)
|
|
|
734 |
enddo
|
|
|
735 |
enddo
|
|
|
736 |
|
|
|
737 |
endif
|
|
|
738 |
|
|
|
739 |
enddo
|
|
|
740 |
|
|
|
741 |
c Write new fields for tests
|
|
|
742 |
if ( test.eq.1 ) then
|
|
|
743 |
|
|
|
744 |
call cdfwopn(zm_fn,cdfid,ierr)
|
|
|
745 |
if (ierr.ne.0) goto 994
|
|
|
746 |
|
|
|
747 |
isok=0
|
|
|
748 |
varname='P_ANO'
|
|
|
749 |
call check_varok(isok,varname,zm_vnam,zm_nvars)
|
|
|
750 |
if (isok.eq.0) then
|
|
|
751 |
call putdef(cdfid,varname,zm_ndim,zm_mdv,zm_vardim,
|
|
|
752 |
> zm_varmin,zm_varmax,zm_stag,ierr)
|
|
|
753 |
if (ierr.ne.0) goto 994
|
|
|
754 |
endif
|
|
|
755 |
call putdat(cdfid,varname,zm_time,0,an_p3,ierr)
|
|
|
756 |
if (ierr.ne.0) goto 994
|
|
|
757 |
print*,'W P_ANO ',trim(zm_fn)
|
|
|
758 |
|
|
|
759 |
isok=0
|
|
|
760 |
varname='T_ANO'
|
|
|
761 |
call check_varok(isok,varname,zm_vnam,zm_nvars)
|
|
|
762 |
if (isok.eq.0) then
|
|
|
763 |
call putdef(cdfid,varname,zm_ndim,zm_mdv,zm_vardim,
|
|
|
764 |
> zm_varmin,zm_varmax,zm_stag,ierr)
|
|
|
765 |
if (ierr.ne.0) goto 994
|
|
|
766 |
endif
|
|
|
767 |
call putdat(cdfid,varname,zm_time,0,an_t3,ierr)
|
|
|
768 |
if (ierr.ne.0) goto 994
|
|
|
769 |
print*,'W T_ANO ',trim(zm_fn)
|
|
|
770 |
|
|
|
771 |
isok=0
|
|
|
772 |
varname='U_ANO'
|
|
|
773 |
call check_varok(isok,varname,zm_vnam,zm_nvars)
|
|
|
774 |
if (isok.eq.0) then
|
|
|
775 |
call putdef(cdfid,varname,zm_ndim,zm_mdv,zm_vardim,
|
|
|
776 |
> zm_varmin,zm_varmax,zm_stag,ierr)
|
|
|
777 |
if (ierr.ne.0) goto 994
|
|
|
778 |
endif
|
|
|
779 |
call putdat(cdfid,varname,zm_time,0,an_u3,ierr)
|
|
|
780 |
if (ierr.ne.0) goto 994
|
|
|
781 |
print*,'W U_ANO ',trim(zm_fn)
|
|
|
782 |
|
|
|
783 |
isok=0
|
|
|
784 |
varname='V_ANO'
|
|
|
785 |
call check_varok(isok,varname,zm_vnam,zm_nvars)
|
|
|
786 |
if (isok.eq.0) then
|
|
|
787 |
call putdef(cdfid,varname,zm_ndim,zm_mdv,zm_vardim,
|
|
|
788 |
> zm_varmin,zm_varmax,zm_stag,ierr)
|
|
|
789 |
if (ierr.ne.0) goto 994
|
|
|
790 |
endif
|
|
|
791 |
call putdat(cdfid,varname,zm_time,0,an_v3,ierr)
|
|
|
792 |
if (ierr.ne.0) goto 994
|
|
|
793 |
print*,'W V_ANO ',trim(zm_fn)
|
|
|
794 |
|
|
|
795 |
call clscdf(cdfid,ierr)
|
|
|
796 |
if (ierr.ne.0) goto 994
|
|
|
797 |
|
|
|
798 |
endif
|
|
|
799 |
|
|
|
800 |
c Exit point for SOR solver
|
|
|
801 |
120 continue
|
|
|
802 |
|
|
|
803 |
c ----------------------------------------------------------------
|
|
|
804 |
c Change surface pressure and get 3d presure field
|
|
|
805 |
c ----------------------------------------------------------------
|
|
|
806 |
|
|
|
807 |
print*,'C PS'
|
|
|
808 |
|
|
|
809 |
c Change surface pressure: based on PV inversion on GEO
|
|
|
810 |
do i=1,ml_nx
|
|
|
811 |
do j=1,ml_ny
|
|
|
812 |
|
|
|
813 |
c Make vertical profile of pressure available
|
|
|
814 |
n1=0
|
|
|
815 |
do k=1,zm_nz
|
|
|
816 |
n1=n1+1
|
|
|
817 |
p1(n1)=an_p3(i,j,k)
|
|
|
818 |
z1(n1)=zm_z3(i,j,k)
|
|
|
819 |
enddo
|
|
|
820 |
|
|
|
821 |
if ( n1.ne.0 ) then
|
|
|
822 |
|
|
|
823 |
c Keep the original surface pressure (psfc=0)
|
|
|
824 |
if ( abs(psfc).lt.eps ) then
|
|
|
825 |
ml_ps(i,j) = ml_ps(i,j)
|
|
|
826 |
|
|
|
827 |
c Take the full change of surface pressure (psfc=1);
|
|
|
828 |
c Interpolation/extrapolation with a natural cubic spline
|
|
|
829 |
elseif ( abs(psfc-1.).lt.eps ) then
|
|
|
830 |
call spline (z1,p1,n1,1.e30,1.e30,spline_f1)
|
|
|
831 |
call splint (z1,p1,spline_f1,n1,ml_oro(i,j),hh)
|
|
|
832 |
ml_ps(i,j)=ml_ps(i,j) + hh*weight(i,j)
|
|
|
833 |
|
|
|
834 |
c Only take a fractional change of surface pressure
|
|
|
835 |
c Interpolation/extrapolation with a natural cubic spline
|
|
|
836 |
else
|
|
|
837 |
call spline (z1,p1,n1,1.e30,1.e30,spline_f1)
|
|
|
838 |
call splint (z1,p1,spline_f1,n1,ml_oro(i,j),hh)
|
|
|
839 |
ml_ps(i,j)=ml_ps(i,j) + hh*psfc*weight(i,j)
|
|
|
840 |
|
|
|
841 |
endif
|
|
|
842 |
|
|
|
843 |
endif
|
|
|
844 |
|
|
|
845 |
enddo
|
|
|
846 |
enddo
|
|
|
847 |
|
|
|
848 |
c Calculate 3d pressure field on model levels
|
|
|
849 |
print*,'C P'
|
|
|
850 |
do k=1,ml_nz
|
|
|
851 |
do i=1,ml_nx
|
|
|
852 |
do j=1,ml_ny
|
|
|
853 |
if (abs(ml_stag(3)+0.5).lt.eps) then
|
|
|
854 |
ml_p3(i,j,k)=ml_aklay(k)+ml_bklay(k)*ml_ps(i,j)
|
|
|
855 |
else
|
|
|
856 |
ml_p3(i,j,k)=ml_aklev(k)+ml_bklev(k)*ml_ps(i,j)
|
|
|
857 |
endif
|
|
|
858 |
enddo
|
|
|
859 |
enddo
|
|
|
860 |
enddo
|
|
|
861 |
|
|
|
862 |
c ----------------------------------------------------------------
|
|
|
863 |
c Get T,U,V at the new pressure levels, based on the original P file
|
|
|
864 |
c ----------------------------------------------------------------
|
|
|
865 |
|
|
|
866 |
c Interpolate T from original P file to new pressure levels
|
|
|
867 |
print*,'C T (ORIGINAL)'
|
|
|
868 |
do i=1,ml_nx
|
|
|
869 |
do j=1,ml_ny
|
|
|
870 |
|
|
|
871 |
n1=0
|
|
|
872 |
do k=ml_nz,1,-1
|
|
|
873 |
n1=n1+1
|
|
|
874 |
f1(n1)=ml_t3(i,j,k)
|
|
|
875 |
p1(n1)=ml_p0(i,j,k)
|
|
|
876 |
enddo
|
|
|
877 |
call spline (p1,f1,n1,1.e30,1.e30,spline_f1)
|
|
|
878 |
do k=1,ml_nz
|
|
|
879 |
call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh)
|
|
|
880 |
ml_t3(i,j,k)=hh
|
|
|
881 |
enddo
|
|
|
882 |
|
|
|
883 |
enddo
|
|
|
884 |
enddo
|
|
|
885 |
|
|
|
886 |
c Interpolate U from original P file to new pressure levels
|
|
|
887 |
print*,'C U (ORIGINAL)'
|
|
|
888 |
do i=1,ml_nx
|
|
|
889 |
do j=1,ml_ny
|
|
|
890 |
|
|
|
891 |
n1=0
|
|
|
892 |
do k=ml_nz,1,-1
|
|
|
893 |
n1=n1+1
|
|
|
894 |
f1(n1)=ml_u3(i,j,k)
|
|
|
895 |
p1(n1)=ml_p0(i,j,k)
|
|
|
896 |
enddo
|
|
|
897 |
call spline (p1,f1,n1,1.e30,1.e30,spline_f1)
|
|
|
898 |
do k=1,ml_nz
|
|
|
899 |
call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh)
|
|
|
900 |
ml_u3(i,j,k)=hh
|
|
|
901 |
enddo
|
|
|
902 |
|
|
|
903 |
enddo
|
|
|
904 |
enddo
|
|
|
905 |
|
|
|
906 |
c Interpolate V from original P file to new pressure levels
|
|
|
907 |
print*,'C V (ORIGINAL)'
|
|
|
908 |
do i=1,ml_nx
|
|
|
909 |
do j=1,ml_ny
|
|
|
910 |
|
|
|
911 |
n1=0
|
|
|
912 |
do k=ml_nz,1,-1
|
|
|
913 |
n1=n1+1
|
|
|
914 |
f1(n1)=ml_v3(i,j,k)
|
|
|
915 |
p1(n1)=ml_p0(i,j,k)
|
|
|
916 |
enddo
|
|
|
917 |
call spline (p1,f1,n1,1.e30,1.e30,spline_f1)
|
|
|
918 |
do k=1,ml_nz
|
|
|
919 |
call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh)
|
|
|
920 |
ml_v3(i,j,k)=hh
|
|
|
921 |
enddo
|
|
|
922 |
|
|
|
923 |
enddo
|
|
|
924 |
enddo
|
|
|
925 |
|
|
|
926 |
c ----------------------------------------------------------------
|
|
|
927 |
c Add T,U,V anomalies at the new pressure levels
|
|
|
928 |
c ----------------------------------------------------------------
|
|
|
929 |
|
|
|
930 |
c Change temperature field
|
|
|
931 |
print*,'C T (ANOMALY)'
|
|
|
932 |
do i=1,ml_nx
|
|
|
933 |
do j=1,ml_ny
|
|
|
934 |
|
|
|
935 |
c Make vertical profile of temperature available
|
|
|
936 |
n1=0
|
|
|
937 |
pmax=-100.
|
|
|
938 |
pmin=2000.
|
|
|
939 |
do k=zm_nz,1,-1
|
|
|
940 |
if ((abs(an_t3(i,j,k)-zm_mdv).gt.eps).and.
|
|
|
941 |
> (zm_z3(i,j,k).gt.ml_oro(i,j)).and.
|
|
|
942 |
> (abs(zm_p3(i,j,k)-zm_mdv).gt.eps)) then
|
|
|
943 |
n1=n1+1
|
|
|
944 |
f1(n1)=an_t3(i,j,k)
|
|
|
945 |
p1(n1)=zm_p3(i,j,k)
|
|
|
946 |
if (p1(n1).lt.pmin) pmin=p1(n1)
|
|
|
947 |
if (p1(n1).gt.pmax) pmax=p1(n1)
|
|
|
948 |
endif
|
|
|
949 |
enddo
|
|
|
950 |
pmin=pmin-dpextra*pmin
|
|
|
951 |
pmax=pmax+dpextra+pmax
|
|
|
952 |
|
|
|
953 |
c Cubic spline interpolation
|
|
|
954 |
if (n1.gt.0) then
|
|
|
955 |
call spline (p1,f1,n1,1.e30,1.e30,spline_f1)
|
|
|
956 |
do k=1,ml_nz
|
|
|
957 |
if ( (ml_p3(i,j,k).gt.pmin).and.
|
|
|
958 |
> (ml_p3(i,j,k).lt.pmax) ) then
|
|
|
959 |
call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh)
|
|
|
960 |
ml_t3(i,j,k)=ml_t3(i,j,k) + hh*weight(i,j)
|
|
|
961 |
endif
|
|
|
962 |
enddo
|
|
|
963 |
endif
|
|
|
964 |
|
|
|
965 |
enddo
|
|
|
966 |
enddo
|
|
|
967 |
|
|
|
968 |
c Change zonal velocity field
|
|
|
969 |
print*,'C U'
|
|
|
970 |
do i=1,ml_nx
|
|
|
971 |
do j=1,ml_ny
|
|
|
972 |
|
|
|
973 |
c Make vertical profile of zonal velocity available
|
|
|
974 |
n1=0
|
|
|
975 |
pmax=-100.
|
|
|
976 |
pmin=2000.
|
|
|
977 |
do k=zm_nz,1,-1
|
|
|
978 |
if ((abs(an_u3(i,j,k)-zm_mdv).gt.eps).and.
|
|
|
979 |
> (zm_z3(i,j,k).gt.ml_oro(i,j)).and.
|
|
|
980 |
> (abs(zm_p3(i,j,k)-zm_mdv).gt.eps)) then
|
|
|
981 |
n1=n1+1
|
|
|
982 |
f1(n1)=an_u3(i,j,k)
|
|
|
983 |
p1(n1)=zm_p3(i,j,k)
|
|
|
984 |
if (p1(n1).lt.pmin) pmin=p1(n1)
|
|
|
985 |
if (p1(n1).gt.pmax) pmax=p1(n1)
|
|
|
986 |
endif
|
|
|
987 |
enddo
|
|
|
988 |
pmin=pmin-dpextra*pmin
|
|
|
989 |
pmax=pmax+dpextra*pmax
|
|
|
990 |
|
|
|
991 |
c Cubic spline interpolation
|
|
|
992 |
if (n1.gt.0) then
|
|
|
993 |
call spline (p1,f1,n1,1.e30,1.e30,spline_f1)
|
|
|
994 |
do k=1,ml_nz
|
|
|
995 |
if ( (ml_p3(i,j,k).gt.pmin).and.
|
|
|
996 |
> (ml_p3(i,j,k).lt.pmax) ) then
|
|
|
997 |
call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh)
|
|
|
998 |
ml_u3(i,j,k)=ml_u3(i,j,k) + hh*weight(i,j)
|
|
|
999 |
endif
|
|
|
1000 |
enddo
|
|
|
1001 |
endif
|
|
|
1002 |
|
|
|
1003 |
enddo
|
|
|
1004 |
enddo
|
|
|
1005 |
|
|
|
1006 |
c Change meridional velocity field
|
|
|
1007 |
print*,'C V'
|
|
|
1008 |
do i=1,ml_nx
|
|
|
1009 |
do j=1,ml_ny
|
|
|
1010 |
|
|
|
1011 |
c Make vertical profile of zonal velocity available
|
|
|
1012 |
n1=0
|
|
|
1013 |
pmax=-100.
|
|
|
1014 |
pmin=2000.
|
|
|
1015 |
do k=zm_nz,1,-1
|
|
|
1016 |
if ((abs(an_v3(i,j,k)-zm_mdv).gt.eps).and.
|
|
|
1017 |
> (zm_z3(i,j,k).gt.ml_oro(i,j)).and.
|
|
|
1018 |
> (abs(zm_p3(i,j,k)-zm_mdv).gt.eps)) then
|
|
|
1019 |
n1=n1+1
|
|
|
1020 |
f1(n1)=an_v3(i,j,k)
|
|
|
1021 |
p1(n1)=zm_p3(i,j,k)
|
|
|
1022 |
if (p1(n1).lt.pmin) pmin=p1(n1)
|
|
|
1023 |
if (p1(n1).gt.pmax) pmax=p1(n1)
|
|
|
1024 |
endif
|
|
|
1025 |
enddo
|
|
|
1026 |
pmin=pmin-dpextra*pmin
|
|
|
1027 |
pmax=pmax+dpextra*pmax
|
|
|
1028 |
|
|
|
1029 |
c Cubic spline interpolation
|
|
|
1030 |
if (n1.gt.0) then
|
|
|
1031 |
call spline (p1,f1,n1,1.e30,1.e30,spline_f1)
|
|
|
1032 |
do k=1,ml_nz
|
|
|
1033 |
if ( (ml_p3(i,j,k).gt.pmin).and.
|
|
|
1034 |
> (ml_p3(i,j,k).lt.pmax) ) then
|
|
|
1035 |
call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh)
|
|
|
1036 |
ml_v3(i,j,k)=ml_v3(i,j,k) + hh*weight(i,j)
|
|
|
1037 |
endif
|
|
|
1038 |
enddo
|
|
|
1039 |
endif
|
|
|
1040 |
|
|
|
1041 |
enddo
|
|
|
1042 |
enddo
|
|
|
1043 |
|
|
|
1044 |
c ---------------------------------------------------------------------
|
|
|
1045 |
c Change geopotential height
|
|
|
1046 |
c ---------------------------------------------------------------------
|
|
|
1047 |
|
|
|
1048 |
c Interpolate Z from original P file to new pressure levels
|
|
|
1049 |
print*,'C Z (ORIGINAL)'
|
|
|
1050 |
do i=1,ml_nx
|
|
|
1051 |
do j=1,ml_ny
|
|
|
1052 |
|
|
|
1053 |
n1=0
|
|
|
1054 |
do k=ml_nz,1,-1
|
|
|
1055 |
n1=n1+1
|
|
|
1056 |
f1(n1)=ml_z3(i,j,k)
|
|
|
1057 |
p1(n1)=ml_p0(i,j,k)
|
|
|
1058 |
enddo
|
|
|
1059 |
call spline (p1,f1,n1,1.e30,1.e30,spline_f1)
|
|
|
1060 |
do k=1,ml_nz
|
|
|
1061 |
call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh)
|
|
|
1062 |
ml_z3(i,j,k)=hh
|
|
|
1063 |
enddo
|
|
|
1064 |
|
|
|
1065 |
enddo
|
|
|
1066 |
enddo
|
|
|
1067 |
|
|
|
1068 |
c Calculate 3d virtual temperature
|
|
|
1069 |
print*,'C TV'
|
|
|
1070 |
do k=1,ml_nz
|
|
|
1071 |
do i=1,ml_nx
|
|
|
1072 |
do j=1,ml_ny
|
|
|
1073 |
ml_tv3(i,j,k) = (ml_t3(i,j,k)+tzero)*
|
|
|
1074 |
> (1.+kappa*ml_q3(i,j,k))
|
|
|
1075 |
enddo
|
|
|
1076 |
enddo
|
|
|
1077 |
enddo
|
|
|
1078 |
|
|
|
1079 |
c Add geopotential height anomaly: first, the pressure anomaly is
|
|
|
1080 |
c interpolated to the new position, then it is transformed into
|
|
|
1081 |
c a correction of geopotential height with the hydrostatic equation
|
|
|
1082 |
print*,'C DZ (MODIFIED -ORIGINAL)'
|
|
|
1083 |
do i=1,ml_nx
|
|
|
1084 |
do j=1,ml_ny
|
|
|
1085 |
|
|
|
1086 |
c Make vertical profile of pressure available
|
|
|
1087 |
n1=0
|
|
|
1088 |
pmax=-100.
|
|
|
1089 |
pmin=2000.
|
|
|
1090 |
do k=zm_nz,1,-1
|
|
|
1091 |
if ((abs(an_p3(i,j,k)-zm_mdv).gt.eps).and.
|
|
|
1092 |
> (zm_z3(i,j,k).gt.ml_oro(i,j)).and.
|
|
|
1093 |
> (abs(zm_p3(i,j,k)-zm_mdv).gt.eps)) then
|
|
|
1094 |
n1=n1+1
|
|
|
1095 |
f1(n1)=an_p3(i,j,k)
|
|
|
1096 |
p1(n1)=zm_p3(i,j,k)
|
|
|
1097 |
if (p1(n1).lt.pmin) pmin=p1(n1)
|
|
|
1098 |
if (p1(n1).gt.pmax) pmax=p1(n1)
|
|
|
1099 |
endif
|
|
|
1100 |
enddo
|
|
|
1101 |
pmin=pmin-dpextra*pmin
|
|
|
1102 |
pmax=pmax+dpextra*pmax
|
|
|
1103 |
|
|
|
1104 |
c Cubic spline interpolation and conversion dp -> dz
|
|
|
1105 |
if (n1.gt.0) then
|
|
|
1106 |
call spline (p1,f1,n1,1.e30,1.e30,spline_f1)
|
|
|
1107 |
do k=1,ml_nz
|
|
|
1108 |
if ( (ml_p3(i,j,k).gt.pmin).and.
|
|
|
1109 |
> (ml_p3(i,j,k).lt.pmax) ) then
|
|
|
1110 |
call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh)
|
|
|
1111 |
hh = -rdg * ml_tv3(i,j,k) * hh/ml_p3(i,j,k)
|
|
|
1112 |
ml_z3(i,j,k)=ml_z3(i,j,k) + hh*weight(i,j)
|
|
|
1113 |
endif
|
|
|
1114 |
enddo
|
|
|
1115 |
endif
|
|
|
1116 |
|
|
|
1117 |
enddo
|
|
|
1118 |
enddo
|
|
|
1119 |
|
|
|
1120 |
c ----------------------------------------------------------------
|
|
|
1121 |
c Remove unrealistic values from final fields
|
|
|
1122 |
c ----------------------------------------------------------------
|
|
|
1123 |
|
|
|
1124 |
if ( poisson.ne.1 ) goto 120
|
|
|
1125 |
|
|
|
1126 |
c Define region for mdv filling
|
|
|
1127 |
count1 = 0
|
|
|
1128 |
|
|
|
1129 |
do i=1,ml_nx
|
|
|
1130 |
do j=1,ml_ny
|
|
|
1131 |
do k=1,ml_nz
|
|
|
1132 |
|
|
|
1133 |
c Get neighbour of grid point
|
|
|
1134 |
il = i-1
|
|
|
1135 |
if (il.lt.1) il=1
|
|
|
1136 |
ir = i+1
|
|
|
1137 |
if ( ir.gt.ml_nx ) ir=ml_nx
|
|
|
1138 |
jd = j-1
|
|
|
1139 |
if (jd.lt.1) jd=1
|
|
|
1140 |
ju = j+1
|
|
|
1141 |
if ( ju.gt.ml_ny ) ju=ml_ny
|
|
|
1142 |
|
|
|
1143 |
c Set flag 2 for boundary and exterior points
|
|
|
1144 |
if ( (abs(zm_rlat(il, j)-zm_mdv).lt.eps).or.
|
|
|
1145 |
> (abs(zm_rlat(ir, j)-zm_mdv).lt.eps).or.
|
|
|
1146 |
> (abs(zm_rlat(i , j)-zm_mdv).lt.eps).or.
|
|
|
1147 |
> (abs(zm_rlat(il,ju)-zm_mdv).lt.eps).or.
|
|
|
1148 |
> (abs(zm_rlat(ir,ju)-zm_mdv).lt.eps).or.
|
|
|
1149 |
> (abs(zm_rlat(i ,ju)-zm_mdv).lt.eps).or.
|
|
|
1150 |
> (abs(zm_rlat(il,jd)-zm_mdv).lt.eps).or.
|
|
|
1151 |
> (abs(zm_rlat(ir,jd)-zm_mdv).lt.eps).or.
|
|
|
1152 |
> (abs(zm_rlat(i ,jd)-zm_mdv).lt.eps) )
|
|
|
1153 |
> then
|
|
|
1154 |
flag_ml(i,j,k) = 2
|
|
|
1155 |
|
|
|
1156 |
c Set flag 1 for interior unphysical points
|
|
|
1157 |
elseif ( ( abs(ml_t3(i,j,k)).gt.500. ).or.
|
|
|
1158 |
> ( abs(ml_u3(i,j,k)).gt.500. ).or.
|
|
|
1159 |
> ( abs(ml_v3(i,j,k)).gt.500. ) )
|
|
|
1160 |
> then
|
|
|
1161 |
flag_ml(i,j,k) = 1
|
|
|
1162 |
count1 = count1 + 1
|
|
|
1163 |
|
|
|
1164 |
c Set flag 0 for interior valid points
|
|
|
1165 |
else
|
|
|
1166 |
flag_ml(i,j,k) = 0
|
|
|
1167 |
endif
|
|
|
1168 |
|
|
|
1169 |
enddo
|
|
|
1170 |
enddo
|
|
|
1171 |
enddo
|
|
|
1172 |
|
|
|
1173 |
print*,'C MASK UNREALISTIC VALUES FOR T,U,V',count1
|
|
|
1174 |
|
|
|
1175 |
c Apply mdv filling
|
|
|
1176 |
print*,'C T (POISSON FILLING)'
|
|
|
1177 |
call mdvfill(ml_t3,ml_t3,flag_ml,ml_nx,ml_ny,ml_nz,10)
|
|
|
1178 |
print*,'C U (POISSON FILLING)'
|
|
|
1179 |
call mdvfill(ml_u3,ml_u3,flag_ml,ml_nx,ml_ny,ml_nz,10)
|
|
|
1180 |
print*,'C V (POISSON FILLING)'
|
|
|
1181 |
call mdvfill(ml_v3,ml_v3,flag_ml,ml_nx,ml_ny,ml_nz,10)
|
|
|
1182 |
|
|
|
1183 |
c ----------------------------------------------------------------
|
|
|
1184 |
c Write modified fields to P file
|
|
|
1185 |
c ----------------------------------------------------------------
|
|
|
1186 |
|
|
|
1187 |
c Open output file
|
|
|
1188 |
call cdfwopn(ml_fn,cdfid,ierr)
|
|
|
1189 |
if (ierr.ne.0) goto 994
|
|
|
1190 |
|
|
|
1191 |
c Write surface pressure
|
|
|
1192 |
isok=0
|
|
|
1193 |
varname='PS'
|
|
|
1194 |
call check_varok(isok,varname,ml_vnam,ml_nvars)
|
|
|
1195 |
if (isok.eq.0) then
|
|
|
1196 |
ml_vardim(3)=3
|
|
|
1197 |
call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
|
|
|
1198 |
> ml_varmin,ml_varmax,ml_stag,ierr)
|
|
|
1199 |
if (ierr.ne.0) goto 994
|
|
|
1200 |
ml_vardim(3)=4
|
|
|
1201 |
endif
|
|
|
1202 |
call putdat(cdfid,varname,ml_time,0,ml_ps,ierr)
|
|
|
1203 |
if (ierr.ne.0) goto 994
|
|
|
1204 |
print*,'W PS ',trim(ml_fn)
|
|
|
1205 |
|
|
|
1206 |
c Write temperature
|
|
|
1207 |
isok=0
|
|
|
1208 |
varname='T'
|
|
|
1209 |
call check_varok(isok,varname,ml_vnam,ml_nvars)
|
|
|
1210 |
if (isok.eq.0) then
|
|
|
1211 |
call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
|
|
|
1212 |
> ml_varmin,ml_varmax,ml_stag,ierr)
|
|
|
1213 |
if (ierr.ne.0) goto 994
|
|
|
1214 |
endif
|
|
|
1215 |
call putdat(cdfid,varname,ml_time,0,ml_t3,ierr)
|
|
|
1216 |
if (ierr.ne.0) goto 994
|
|
|
1217 |
print*,'W T ',trim(ml_fn)
|
|
|
1218 |
|
|
|
1219 |
c Write zonal wind
|
|
|
1220 |
isok=0
|
|
|
1221 |
varname='U'
|
|
|
1222 |
call check_varok(isok,varname,ml_vnam,ml_nvars)
|
|
|
1223 |
if (isok.eq.0) then
|
|
|
1224 |
call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
|
|
|
1225 |
> ml_varmin,ml_varmax,ml_stag,ierr)
|
|
|
1226 |
if (ierr.ne.0) goto 994
|
|
|
1227 |
endif
|
|
|
1228 |
call putdat(cdfid,varname,ml_time,0,ml_u3,ierr)
|
|
|
1229 |
if (ierr.ne.0) goto 994
|
|
|
1230 |
print*,'W U ',trim(ml_fn)
|
|
|
1231 |
|
|
|
1232 |
c Write meridional wind
|
|
|
1233 |
isok=0
|
|
|
1234 |
varname='V'
|
|
|
1235 |
call check_varok(isok,varname,ml_vnam,ml_nvars)
|
|
|
1236 |
if (isok.eq.0) then
|
|
|
1237 |
call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
|
|
|
1238 |
> ml_varmin,ml_varmax,ml_stag,ierr)
|
|
|
1239 |
if (ierr.ne.0) goto 994
|
|
|
1240 |
endif
|
|
|
1241 |
call putdat(cdfid,varname,ml_time,0,ml_v3,ierr)
|
|
|
1242 |
if (ierr.ne.0) goto 994
|
|
|
1243 |
print*,'W V ',trim(ml_fn)
|
|
|
1244 |
|
|
|
1245 |
c Write geopotential height
|
|
|
1246 |
isok=0
|
|
|
1247 |
varname='Z'
|
|
|
1248 |
call check_varok(isok,varname,ml_vnam,ml_nvars)
|
|
|
1249 |
if (isok.eq.0) then
|
|
|
1250 |
call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
|
|
|
1251 |
> ml_varmin,ml_varmax,ml_stag,ierr)
|
|
|
1252 |
if (ierr.ne.0) goto 994
|
|
|
1253 |
endif
|
|
|
1254 |
call putdat(cdfid,varname,ml_time,0,ml_z3,ierr)
|
|
|
1255 |
if (ierr.ne.0) goto 994
|
|
|
1256 |
print*,'W Z ',trim(ml_fn)
|
|
|
1257 |
|
|
|
1258 |
c Filter matrix (only in test mode)
|
|
|
1259 |
if (test.eq.1) then
|
|
|
1260 |
isok=0
|
|
|
1261 |
varname='DIST'
|
|
|
1262 |
call check_varok(isok,varname,ml_vnam,ml_nvars)
|
|
|
1263 |
if (isok.eq.0) then
|
|
|
1264 |
ml_vardim(3)=1
|
|
|
1265 |
call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
|
|
|
1266 |
> ml_varmin,ml_varmax,ml_stag,ierr)
|
|
|
1267 |
if (ierr.ne.0) goto 994
|
|
|
1268 |
endif
|
|
|
1269 |
call putdat(cdfid,varname,ml_time,1,dist,ierr)
|
|
|
1270 |
if (ierr.ne.0) goto 994
|
|
|
1271 |
print*,'W DIST ',trim(ml_fn)
|
|
|
1272 |
|
|
|
1273 |
isok=0
|
|
|
1274 |
varname='WEIGHT'
|
|
|
1275 |
call check_varok(isok,varname,ml_vnam,ml_nvars)
|
|
|
1276 |
if (isok.eq.0) then
|
|
|
1277 |
ml_vardim(3)=1
|
|
|
1278 |
call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,
|
|
|
1279 |
> ml_varmin,ml_varmax,ml_stag,ierr)
|
|
|
1280 |
if (ierr.ne.0) goto 994
|
|
|
1281 |
endif
|
|
|
1282 |
call putdat(cdfid,varname,ml_time,1,weight,ierr)
|
|
|
1283 |
if (ierr.ne.0) goto 994
|
|
|
1284 |
print*,'W WEIGHT ',trim(ml_fn)
|
|
|
1285 |
|
|
|
1286 |
endif
|
|
|
1287 |
|
|
|
1288 |
c Close output file
|
|
|
1289 |
call clscdf(cdfid,ierr)
|
|
|
1290 |
if (ierr.ne.0) goto 994
|
|
|
1291 |
|
|
|
1292 |
c ----------------------------------------------------------------
|
|
|
1293 |
c Exception handling
|
|
|
1294 |
c ----------------------------------------------------------------
|
|
|
1295 |
|
|
|
1296 |
stop
|
|
|
1297 |
|
|
|
1298 |
998 print*,'Problem with input netcdf file ',trim(ml_fn)
|
|
|
1299 |
stop
|
|
|
1300 |
|
|
|
1301 |
997 print*,'Problem with input netcdf file ',trim(zm_fn)
|
|
|
1302 |
stop
|
|
|
1303 |
|
|
|
1304 |
994 print*,'Problem with output netcdf file ',trim(ml_fn)
|
|
|
1305 |
stop
|
|
|
1306 |
|
|
|
1307 |
end
|
|
|
1308 |
|
|
|
1309 |
|
|
|
1310 |
c ****************************************************************
|
|
|
1311 |
c * SUBROUTINE SECTION: AUXILIARY ROUTINES *
|
|
|
1312 |
c ****************************************************************
|
|
|
1313 |
|
|
|
1314 |
c ----------------------------------------------------------------
|
|
|
1315 |
c Check whether variable is found on netcdf file
|
|
|
1316 |
c ----------------------------------------------------------------
|
|
|
1317 |
|
|
|
1318 |
subroutine check_varok (isok,varname,varlist,nvars)
|
|
|
1319 |
|
|
|
1320 |
c Check whether the variable <varname> is in the list <varlist(nvars)>.
|
|
|
1321 |
c If this is the case, <isok> is incremented by 1. Otherwise <isok>
|
|
|
1322 |
c keeps its value.
|
|
|
1323 |
|
|
|
1324 |
implicit none
|
|
|
1325 |
|
|
|
1326 |
c Declaraion of subroutine parameters
|
|
|
1327 |
integer isok
|
|
|
1328 |
integer nvars
|
|
|
1329 |
character*80 varname
|
|
|
1330 |
character*80 varlist(nvars)
|
|
|
1331 |
|
|
|
1332 |
c Auxiliary variables
|
|
|
1333 |
integer i
|
|
|
1334 |
|
|
|
1335 |
c Main
|
|
|
1336 |
do i=1,nvars
|
|
|
1337 |
if (trim(varname).eq.trim(varlist(i))) isok=isok+1
|
|
|
1338 |
enddo
|
|
|
1339 |
|
|
|
1340 |
end
|
|
|
1341 |
|
|
|
1342 |
c -------------------------------------------------------------
|
|
|
1343 |
c Natural cubic spline /from Numerical Recipes)
|
|
|
1344 |
c -------------------------------------------------------------
|
|
|
1345 |
|
|
|
1346 |
SUBROUTINE spline(x,y,n,yp1,ypn,y2)
|
|
|
1347 |
INTEGER n,NMAX
|
|
|
1348 |
REAL yp1,ypn,x(n),y(n),y2(n)
|
|
|
1349 |
PARAMETER (NMAX=500)
|
|
|
1350 |
INTEGER i,k
|
|
|
1351 |
REAL p,qn,sig,un,u(NMAX)
|
|
|
1352 |
if (yp1.gt..99e30) then
|
|
|
1353 |
y2(1)=0.
|
|
|
1354 |
u(1)=0.
|
|
|
1355 |
else
|
|
|
1356 |
y2(1)=-0.5
|
|
|
1357 |
u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
|
|
|
1358 |
endif
|
|
|
1359 |
do 11 i=2,n-1
|
|
|
1360 |
sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
|
|
|
1361 |
p=sig*y2(i-1)+2.
|
|
|
1362 |
y2(i)=(sig-1.)/p
|
|
|
1363 |
u(i)=(6.*((y(i+1)-y(i))/(x(i+
|
|
|
1364 |
*1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
|
|
|
1365 |
*u(i-1))/p
|
|
|
1366 |
11 continue
|
|
|
1367 |
if (ypn.gt..99e30) then
|
|
|
1368 |
qn=0.
|
|
|
1369 |
un=0.
|
|
|
1370 |
else
|
|
|
1371 |
qn=0.5
|
|
|
1372 |
un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
|
|
|
1373 |
endif
|
|
|
1374 |
y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
|
|
|
1375 |
do 12 k=n-1,1,-1
|
|
|
1376 |
y2(k)=y2(k)*y2(k+1)+u(k)
|
|
|
1377 |
12 continue
|
|
|
1378 |
return
|
|
|
1379 |
END
|
|
|
1380 |
|
|
|
1381 |
|
|
|
1382 |
SUBROUTINE splint(xa,ya,y2a,n,x,y)
|
|
|
1383 |
INTEGER n
|
|
|
1384 |
REAL x,y,xa(n),y2a(n),ya(n)
|
|
|
1385 |
INTEGER k,khi,klo
|
|
|
1386 |
REAL a,b,h
|
|
|
1387 |
klo=1
|
|
|
1388 |
khi=n
|
|
|
1389 |
1 if (khi-klo.gt.1) then
|
|
|
1390 |
k=(khi+klo)/2
|
|
|
1391 |
if(xa(k).gt.x)then
|
|
|
1392 |
khi=k
|
|
|
1393 |
else
|
|
|
1394 |
klo=k
|
|
|
1395 |
endif
|
|
|
1396 |
goto 1
|
|
|
1397 |
endif
|
|
|
1398 |
h=xa(khi)-xa(klo)
|
|
|
1399 |
if (h.eq.0.) pause 'bad xa input in splint'
|
|
|
1400 |
a=(xa(khi)-x)/h
|
|
|
1401 |
b=(x-xa(klo))/h
|
|
|
1402 |
y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
|
|
|
1403 |
*2)/6.
|
|
|
1404 |
return
|
|
|
1405 |
END
|
|
|
1406 |
|
|
|
1407 |
|
|
|
1408 |
c ---------------------------------------------------------
|
|
|
1409 |
c Spherical distance between two lat/lon points
|
|
|
1410 |
c ---------------------------------------------------------
|
|
|
1411 |
|
|
|
1412 |
real function sdis(xp,yp,xq,yq)
|
|
|
1413 |
c
|
|
|
1414 |
c calculates spherical distance (in km) between two points given
|
|
|
1415 |
c by their spherical coordinates (xp,yp) and (xq,yq), respectively.
|
|
|
1416 |
c
|
|
|
1417 |
real re
|
|
|
1418 |
parameter (re=6370.)
|
|
|
1419 |
real xp,yp,xq,yq,arg
|
|
|
1420 |
|
|
|
1421 |
arg=sind(yp)*sind(yq)+cosd(yp)*cosd(yq)*cosd(xp-xq)
|
|
|
1422 |
if (arg.lt.-1.) arg=-1.
|
|
|
1423 |
if (arg.gt.1.) arg=1.
|
|
|
1424 |
sdis=re*acos(arg)
|
|
|
1425 |
|
|
|
1426 |
end
|
|
|
1427 |
|
|
|
1428 |
C -----------------------------------------------------------------
|
|
|
1429 |
c Rene's KINK function (for smoothing at bounadries)
|
|
|
1430 |
C -----------------------------------------------------------------
|
|
|
1431 |
|
|
|
1432 |
real function kink (x,a)
|
|
|
1433 |
|
|
|
1434 |
implicit none
|
|
|
1435 |
|
|
|
1436 |
c declaration of parameters
|
|
|
1437 |
real x,a
|
|
|
1438 |
|
|
|
1439 |
c parameters
|
|
|
1440 |
real pi
|
|
|
1441 |
parameter (pi=3.1415926535)
|
|
|
1442 |
|
|
|
1443 |
if (x.lt.0.) then
|
|
|
1444 |
kink=0.
|
|
|
1445 |
elseif (x.gt.a) then
|
|
|
1446 |
kink=1.
|
|
|
1447 |
else
|
|
|
1448 |
kink=(1.+cos(pi*(x-a)/a))/2.
|
|
|
1449 |
endif
|
|
|
1450 |
|
|
|
1451 |
return
|
|
|
1452 |
end
|
|
|
1453 |
|
|
|
1454 |
C -----------------------------------------------------------------
|
|
|
1455 |
c Rene's KINK function (for smoothing at bounadries)
|
|
|
1456 |
C -----------------------------------------------------------------
|
|
|
1457 |
|
|
|
1458 |
subroutine mdvfill (out,inp,flag,nx,ny,nz,maxiter)
|
|
|
1459 |
|
|
|
1460 |
implicit none
|
|
|
1461 |
|
|
|
1462 |
c Declaration of subroutine parameters
|
|
|
1463 |
integer nx,ny,nz
|
|
|
1464 |
real inp (nx,ny,nz)
|
|
|
1465 |
real out (nx,ny,nz)
|
|
|
1466 |
integer flag(nx,ny,nz)
|
|
|
1467 |
integer maxiter
|
|
|
1468 |
|
|
|
1469 |
c Parameters
|
|
|
1470 |
real omega ! Omega fopr SOR
|
|
|
1471 |
parameter (omega=1.5)
|
|
|
1472 |
|
|
|
1473 |
c Auxiliary variables
|
|
|
1474 |
integer i,j,k
|
|
|
1475 |
integer iter
|
|
|
1476 |
real tmp0(nx,ny,nz)
|
|
|
1477 |
real tmp1(nx,ny,nz)
|
|
|
1478 |
integer il,ir,ju,jd
|
|
|
1479 |
integer count
|
|
|
1480 |
real mean(nz)
|
|
|
1481 |
|
|
|
1482 |
c Calculate mean of variable for all levels
|
|
|
1483 |
do k=1,nz
|
|
|
1484 |
mean(k) = 0.
|
|
|
1485 |
count = 0
|
|
|
1486 |
do i=1,nx
|
|
|
1487 |
do j=1,ny
|
|
|
1488 |
if ( flag(i,j,k).eq.0 ) then
|
|
|
1489 |
count = count + 1
|
|
|
1490 |
mean(k) = mean(k) + inp(i,j,k)
|
|
|
1491 |
endif
|
|
|
1492 |
enddo
|
|
|
1493 |
enddo
|
|
|
1494 |
if ( count.ne.0 ) then
|
|
|
1495 |
mean(k) = mean(k)/real(count)
|
|
|
1496 |
else
|
|
|
1497 |
mean(k) = 0.
|
|
|
1498 |
endif
|
|
|
1499 |
enddo
|
|
|
1500 |
|
|
|
1501 |
c Create first guess
|
|
|
1502 |
do k=1,nz
|
|
|
1503 |
do i=1,nx
|
|
|
1504 |
do j=1,ny
|
|
|
1505 |
if ( flag(i,j,k).eq.0 ) then
|
|
|
1506 |
tmp0(i,j,k) = inp(i,j,k)
|
|
|
1507 |
else
|
|
|
1508 |
tmp0(i,j,k) = mean(k)
|
|
|
1509 |
endif
|
|
|
1510 |
enddo
|
|
|
1511 |
enddo
|
|
|
1512 |
enddo
|
|
|
1513 |
|
|
|
1514 |
c SOR iterations
|
|
|
1515 |
iter = 0
|
|
|
1516 |
122 continue
|
|
|
1517 |
|
|
|
1518 |
c Loop over all points
|
|
|
1519 |
do k=1,nz
|
|
|
1520 |
do i=1,nx
|
|
|
1521 |
do j=1,ny
|
|
|
1522 |
|
|
|
1523 |
c Apply the updating only for specified points
|
|
|
1524 |
if ( flag(i,j,k).ne.1 ) goto 121
|
|
|
1525 |
|
|
|
1526 |
c Get neighbouring points - no handling of periodic domains!
|
|
|
1527 |
il = i-1
|
|
|
1528 |
if (il.lt.1) il=1
|
|
|
1529 |
ir = i+1
|
|
|
1530 |
if ( ir.gt.nx ) ir=nx
|
|
|
1531 |
jd = j-1
|
|
|
1532 |
if (jd.lt.1) jd=1
|
|
|
1533 |
ju = j+1
|
|
|
1534 |
if ( ju.gt.ny ) ju=ny
|
|
|
1535 |
|
|
|
1536 |
c Update field
|
|
|
1537 |
tmp1(i,j,k) = 0.25 * ( tmp0(il,j,k) + tmp0(ir,j,k) +
|
|
|
1538 |
> tmp0(i,ju,k) + tmp0(i,jd,k) )
|
|
|
1539 |
|
|
|
1540 |
tmp0(i,j,k) = omega * tmp1(i,j,k) +
|
|
|
1541 |
> (1. - omega) * tmp1(i,j,k)
|
|
|
1542 |
|
|
|
1543 |
c Exit point for loop
|
|
|
1544 |
121 continue
|
|
|
1545 |
|
|
|
1546 |
enddo
|
|
|
1547 |
enddo
|
|
|
1548 |
enddo
|
|
|
1549 |
|
|
|
1550 |
c Decide whether further iterations are needed
|
|
|
1551 |
iter = iter + 1
|
|
|
1552 |
if ( iter.lt.maxiter) goto 122
|
|
|
1553 |
|
|
|
1554 |
c Save output
|
|
|
1555 |
do i=1,nx
|
|
|
1556 |
do j=1,ny
|
|
|
1557 |
do k=1,nz
|
|
|
1558 |
if ( flag(i,j,k).eq.1 ) then
|
|
|
1559 |
out(i,j,k) = tmp0(i,j,k)
|
|
|
1560 |
endif
|
|
|
1561 |
enddo
|
|
|
1562 |
enddo
|
|
|
1563 |
enddo
|
|
|
1564 |
|
|
|
1565 |
end
|
|
|
1566 |
|
|
|
1567 |
c
|