1 |
subroutine applyH(mem,xy,iobsloc,currobs) |
2 |
implicit none |
3 |
#include "osse_parms.h" |
4 |
! Arguments |
5 |
integer, intent(in) :: iobsloc(mobs,5) |
6 |
integer, intent(in) :: currobs |
7 |
real*8, intent(in) :: mem(n) |
8 |
real*8, intent(out) :: xy(mobs) |
9 |
|
10 |
integer modelloc(mobs,5) |
11 |
real*8 ampin(currobs/2) |
12 |
real*8 ampout(currobs/2) |
13 |
real*8 amperr(currobs/2) |
14 |
real*8 correction |
15 |
integer i,indx |
16 |
integer method |
17 |
real theta |
18 |
|
19 |
method=2 |
20 |
c do i=1,n |
21 |
c write(*,*) mem(i) |
22 |
c enddo |
23 |
|
24 |
c write(*,*) 'in applyH' |
25 |
call getmodelloc(modelloc,iobsloc,currobs) |
26 |
c write(*,*) 'past getmodelloc' |
27 |
|
28 |
if (method==1) then |
29 |
|
30 |
do i=1,currobs |
31 |
indx=modelloc(i,1) |
32 |
xy(i)=mem(indx) |
33 |
c write(*,*) xy(i), modelloc(i,5),modelloc(i,4),modelloc(i,3), |
34 |
c & i |
35 |
enddo |
36 |
|
37 |
else if (method==2) then |
38 |
do i=1,currobs |
39 |
c write(*,*) mem(i) |
40 |
indx=modelloc(i,1) |
41 |
c theta=-atan2(real(iobsloc(i,5))-,real(iobsloc(i,5))); |
42 |
c theta=-atan2(real(iobsloc(i,5))-,real(iobsloc(i,5))); |
43 |
theta=-atan2(real(iobsloc(i,5))-xorigin, |
44 |
& real(iobsloc(i,4))-yorigin) |
45 |
c theta=theta+pi |
46 |
c if (modelloc(i,5)>nx .or. modelloc(i,4)>ny .or. |
47 |
c & modelloc(i,4)<9) then |
48 |
c write(*,*) i, modelloc(i,5),modelloc(i,4) |
49 |
c xy(i)=0.0 |
50 |
c elseif (iobsloc(i,2)==0) then ! azim |
51 |
c following imaged with imagesc(rot90(vect(:,:,5,2)));axis xy |
52 |
if (iobsloc(i,2)==0) then ! U |
53 |
c xy(i)=0 |
54 |
c xy(i)=cos(theta) ! azim component (data flipped) |
55 |
c xy(i)=-sin(theta) ! radial component (data flipped) |
56 |
c xy(i)=cos(theta)*mem(indx) |
57 |
c & -sin(theta)*mem(indx+nndata) |
58 |
c xy(i)=sin(theta) ! azim component (data unflipped) |
59 |
c xy(i)=cos(theta) ! radial component (data unflipped) |
60 |
c xy(i)=sin(theta)+cos(theta) ! (data unflipped) |
61 |
c xy(i)=sin(theta)*mem(indx) |
62 |
c & +cos(theta)*mem(indx+nn) ! (data unflipped) |
63 |
c xy(i)=sin(theta)*mem(indx) |
64 |
c & +cos(theta)*mem(modelloc(i+currobs/2,1)) ! (data unflipped) |
65 |
xy(i)=-cos(theta)*mem(indx) |
66 |
& -sin(theta)*mem(modelloc(i+currobs/2,1)) ! (hard way) |
67 |
elseif (iobsloc(i,2)==1) then ! radial |
68 |
c xy(i)=0 |
69 |
c xy(i)=-sin(theta) ! azim (data flipped) |
70 |
c xy(i)=-cos(theta) ! radial (data flipped) |
71 |
c xy(i)=mem(indx-nndata)*-sin(theta) + |
72 |
c & mem(indx)*-cos(theta) |
73 |
c xy(i)=cos(theta) ! azim (data unflipped) |
74 |
c xy(i)=-sin(theta) ! radial (data unflipped) |
75 |
c xy(i)=cos(theta)-sin(theta) ! (data unflipped) |
76 |
c xy(i)=cos(theta)*mem(modelloc(i-currobs/2,1)) - |
77 |
c & sin(theta)*mem(indx) ! (data unflipped) |
78 |
xy(i)=-sin(theta)*mem(modelloc(i-currobs/2,1)) |
79 |
& +cos(theta)*mem(indx) ! (hard way) |
80 |
else |
81 |
xy(i)=mem(indx) |
82 |
endif |
83 |
|
84 |
c write(*,*) xy(i) |
85 |
c if (xy(i)>100 .or. xy(i)<-100) then |
86 |
c write(*,*) modelloc(i,5), modelloc(i,4), modelloc(i,3), |
87 |
c & modelloc(i,2), xy(i), mem(indx+nndata), |
88 |
c & mem(indx),-sin(theta),-cos(theta),theta |
89 |
c stop |
90 |
c write(*,*) i |
91 |
c write(*,*) modelloc(i,5), modelloc(i,4), modelloc(i,3), |
92 |
c & modelloc(i,2), mem(modelloc(i,1)) |
93 |
c & modelloc(i,2), xy(i) |
94 |
c endif |
95 |
c write(*,*) iobsloc(i,5), iobsloc(i,4), iobsloc(i,3), |
96 |
c & iobsloc(i,2), xy(i) |
97 |
enddo |
98 |
do i=1,currobs/2 |
99 |
indx=modelloc(i,1) |
100 |
ampin(i)=sqrt(mem(indx)**2 + mem(indx+nn)**2) |
101 |
ampout(i)=sqrt(xy(i)**2 + xy(i+(currobs/2))**2) |
102 |
if (ampin(i)==0) then |
103 |
xy(i)=0 |
104 |
xy(i+(currobs/2))=0 |
105 |
else |
106 |
correction=ampout(i)/ampin(i) |
107 |
xy(i)=xy(i)/correction |
108 |
xy(i+(currobs/2))=xy(i+(currobs/2))/correction |
109 |
endif |
110 |
ampin(i)=sqrt(mem(indx)**2 + mem(indx+nn)**2) |
111 |
ampout(i)=sqrt(xy(i)**2 + xy(i+(currobs/2))**2) |
112 |
c write(*,*) i, ampin(i),ampout(i) |
113 |
enddo |
114 |
c amperr=ampout-ampin |
115 |
c write(*,*) maxval(ampin),maxval(ampout),maxval(amperr) |
116 |
|
117 |
|
118 |
else |
119 |
stop 'bad H method' |
120 |
endif |
121 |
c stop |
122 |
end subroutine applyH |
123 |
|
124 |
|
125 |
subroutine unH(mem,xy,iobsloc,currobs) |
126 |
implicit none |
127 |
#include "osse_parms.h" |
128 |
! Arguments |
129 |
integer, intent(in) :: iobsloc(mobs,5) |
130 |
integer, intent(in) :: currobs |
131 |
real*8, intent(in) :: mem(n) |
132 |
real*8, intent(out) :: xy(mobs) |
133 |
|
134 |
integer modelloc(mobs,5) |
135 |
real*8 ampin(currobs/2) |
136 |
real*8 ampout(currobs/2) |
137 |
real*8 amperr(currobs/2) |
138 |
real*8 correction |
139 |
integer i,indx |
140 |
integer method |
141 |
real theta |
142 |
|
143 |
method=2 |
144 |
c do i=1,n |
145 |
c write(*,*) mem(i) |
146 |
c enddo |
147 |
|
148 |
c write(*,*) 'in applyH' |
149 |
call getmodelloc(modelloc,iobsloc,currobs) |
150 |
c write(*,*) 'past getmodelloc' |
151 |
|
152 |
if (method==1) then |
153 |
|
154 |
do i=1,currobs |
155 |
indx=modelloc(i,1) |
156 |
xy(i)=mem(indx) |
157 |
c write(*,*) xy(i), modelloc(i,5),modelloc(i,4),modelloc(i,3), |
158 |
c & i |
159 |
enddo |
160 |
|
161 |
else if (method==2) then |
162 |
do i=1,currobs |
163 |
c write(*,*) mem(i) |
164 |
c indx=modelloc(i,1) |
165 |
c theta=-atan2(real(iobsloc(i,5))-,real(iobsloc(i,5))); |
166 |
c theta=-atan2(real(iobsloc(i,5))-,real(iobsloc(i,5))); |
167 |
theta=-atan2(real(iobsloc(i,5))-xorigin, |
168 |
& real(iobsloc(i,4))-yorigin) |
169 |
c theta=theta+pi |
170 |
c if (modelloc(i,5)>nx .or. modelloc(i,4)>ny .or. |
171 |
c & modelloc(i,4)<9) then |
172 |
c write(*,*) i, modelloc(i,5),modelloc(i,4) |
173 |
c xy(i)=0.0 |
174 |
c elseif (iobsloc(i,2)==0) then ! azim |
175 |
c following imaged with imagesc(rot90(vect(:,:,5,2)));axis xy |
176 |
if (iobsloc(i,2)==0) then ! azimuthal (model U) |
177 |
xy(i)=-cos(theta)*mem(i) |
178 |
& -sin(theta)*mem(i+currobs/2) ! (hard way) |
179 |
elseif (iobsloc(i,2)==1) then ! radial (model V) |
180 |
xy(i)=-sin(theta)*mem(i-currobs/2) |
181 |
& +cos(theta)*mem(i) ! (hard way) |
182 |
else |
183 |
xy(i)=mem(i) |
184 |
endif |
185 |
|
186 |
c write(*,*) xy(i) |
187 |
c if (xy(i)>100 .or. xy(i)<-100) then |
188 |
c write(*,*) modelloc(i,5), modelloc(i,4), modelloc(i,3), |
189 |
c & modelloc(i,2), xy(i), mem(indx+nndata), |
190 |
c & mem(indx),-sin(theta),-cos(theta),theta |
191 |
c stop |
192 |
c write(*,*) i |
193 |
c write(*,*) modelloc(i,5), modelloc(i,4), modelloc(i,3), |
194 |
c & modelloc(i,2), mem(modelloc(i,1)) |
195 |
c & modelloc(i,2), xy(i) |
196 |
c endif |
197 |
c write(*,*) iobsloc(i,5), iobsloc(i,4), iobsloc(i,3), |
198 |
c & iobsloc(i,2), xy(i) |
199 |
enddo |
200 |
if(0) then |
201 |
do i=1,currobs/2 |
202 |
indx=modelloc(i,1) |
203 |
ampin(i)=sqrt(mem(indx)**2 + mem(indx+nn)**2) |
204 |
ampout(i)=sqrt(xy(i)**2 + xy(i+(currobs/2))**2) |
205 |
if (ampin(i)==0) then |
206 |
xy(i)=0 |
207 |
xy(i+(currobs/2))=0 |
208 |
else |
209 |
correction=ampout(i)/ampin(i) |
210 |
xy(i)=xy(i)/correction |
211 |
xy(i+(currobs/2))=xy(i+(currobs/2))/correction |
212 |
endif |
213 |
ampin(i)=sqrt(mem(indx)**2 + mem(indx+nn)**2) |
214 |
ampout(i)=sqrt(xy(i)**2 + xy(i+(currobs/2))**2) |
215 |
c write(*,*) i, ampin(i),ampout(i) |
216 |
enddo |
217 |
endif |
218 |
c amperr=ampout-ampin |
219 |
c write(*,*) maxval(ampin),maxval(ampout),maxval(amperr) |
220 |
|
221 |
|
222 |
else |
223 |
stop 'bad H method' |
224 |
endif |
225 |
c stop |
226 |
end subroutine unH |
227 |
|
228 |
|
229 |
subroutine cyl2cart(mem,Hmem) |
230 |
implicit none |
231 |
#include "osse_parms.h" |
232 |
! Arguments |
233 |
real*8, intent(in) :: mem(n) |
234 |
real*8, intent(out) :: Hmem(n) |
235 |
|
236 |
real*8 correction |
237 |
integer i,v,z,y,x |
238 |
real theta |
239 |
|
240 |
do i=1,n |
241 |
|
242 |
v=(i-1)/nn |
243 |
z= 1+mod(i-1,nn)/(nx*ny) |
244 |
y= 1+((mod(i-1,nn)-nx*ny*(z-1))/nx) |
245 |
x=1+mod(i-1,nn)-(z-1)*ny*nx-(y-1)*nx |
246 |
|
247 |
theta=2.0*pi*(x-1)/nx |
248 |
|
249 |
if (v==0) then ! azimuthal (model U) |
250 |
Hmem(i)=-cos(theta)*mem(i) |
251 |
& -sin(theta)*mem(i+nn) |
252 |
elseif (v==1) then ! radial (model V) |
253 |
Hmem(i)=-sin(theta)*mem(i-nn) |
254 |
& +cos(theta)*mem(i) ! (hard way) |
255 |
else |
256 |
Hmem(i)=mem(i) |
257 |
endif |
258 |
|
259 |
enddo |
260 |
end subroutine cyl2cart |
261 |
|
262 |
|
263 |
|
264 |
subroutine cart2cyl(Hmem,mem) |
265 |
implicit none |
266 |
#include "osse_parms.h" |
267 |
! Arguments |
268 |
real*8, intent(in) :: Hmem(n) |
269 |
real*8, intent(out) :: mem(n) |
270 |
|
271 |
real*8 correction |
272 |
integer i,v,z,y,x |
273 |
real theta |
274 |
|
275 |
do i=1,n |
276 |
|
277 |
v=(i-1)/nn |
278 |
z= 1+mod(i-1,nn)/(nx*ny) |
279 |
y= 1+((mod(i-1,nn)-nx*ny*(z-1))/nx) |
280 |
x=1+mod(i-1,nn)-(z-1)*ny*nx-(y-1)*nx |
281 |
|
282 |
theta=2.0*pi*(x-1)/nx |
283 |
|
284 |
if (v==0) then ! azimuthal (model U) |
285 |
mem(i)=-cos(theta)*Hmem(i) |
286 |
& -sin(theta)*Hmem(i+nn) |
287 |
elseif (v==1) then ! radial (model V) |
288 |
mem(i)=-sin(theta)*Hmem(i-nn) |
289 |
& +cos(theta)*Hmem(i) ! (hard way) |
290 |
else |
291 |
mem(i)=Hmem(i) |
292 |
endif |
293 |
|
294 |
enddo |
295 |
end subroutine cart2cyl |
296 |
|
297 |
|
298 |
|
299 |
subroutine getmodelloc(modelloc,iobsloc,currobs) |
300 |
|
301 |
implicit none |
302 |
#include "osse_parms.h" |
303 |
! Arguments |
304 |
integer, intent(in) :: iobsloc(mobs,5) |
305 |
integer, intent(in) :: currobs |
306 |
integer, intent(out) :: modelloc(mobs,5) |
307 |
|
308 |
integer i,j,loc,k,method |
309 |
integer yob(mobs),xob(mobs),zob(mobs),vob(mobs) |
310 |
integer Xu(xdata),Xv(xdata),Xt(xdata) |
311 |
integer Yu(ydata),Yv(ydata),Yt(ydata) |
312 |
real Xufract(xdata),Xvfract(xdata) |
313 |
real Yufract(ydata),Yvfract(ydata) |
314 |
real theta(xdata,ydata),thetamod(xdata,ydata) |
315 |
real rho(xdata,ydata),rhomod(xdata,ydata) |
316 |
integer X,Y |
317 |
integer, external :: round |
318 |
|
319 |
c write(*,*) ' in applyH' |
320 |
method=2 ! nearest neighbor, same coordinates |
321 |
|
322 |
c get nearest neighbor locations in model space for all of observation space |
323 |
c tuned to value location on grid box -- u, v and t values |
324 |
c (this is for a C-grid) |
325 |
if (method==1) then |
326 |
|
327 |
|
328 |
c find fractional gridpoint coords in model space of observation locations |
329 |
do k=1,xdata |
330 |
Xufract(k)=(real(k-1)/real(xdata))*real(nx)+1.0 |
331 |
Xvfract(k)=((real(k)-0.5)/real(xdata))*real(nx) |
332 |
enddo |
333 |
|
334 |
do k=1,ydata |
335 |
Yufract(k)=((real(k)-0.5)/real(ydata))*real(ny) |
336 |
Yvfract(k)=(real((k-1))/real(ydata))*real(ny)+1.0 |
337 |
enddo |
338 |
|
339 |
c now find nearest model grid cell locations of observations |
340 |
do k=1,xdata |
341 |
Xu(k)=round(Xufract(k)) |
342 |
Xv(k)=ceiling(Xvfract(k)) |
343 |
enddo |
344 |
Xt=Xv |
345 |
do k=1,ydata |
346 |
Yu(k)=ceiling(Yufract(k)) |
347 |
Yv(k)=round(Yvfract(k)) |
348 |
c write(*,*) k, Yu(k), Yv(k) |
349 |
enddo |
350 |
Yt=Yu |
351 |
|
352 |
c for each obervation, get nearest neighbor value |
353 |
do k=1,currobs |
354 |
|
355 |
modelloc(k,2)=iobsloc(k,2) |
356 |
modelloc(k,3)=iobsloc(k,3) |
357 |
c vob accessed more than once to make reasoning clear, heh |
358 |
if (iobsloc(k,2)==0) then ! u location |
359 |
modelloc(k,1)=(iobsloc(k,2))*nn+(iobsloc(k,3)-1)*nx*ny+ |
360 |
& (Yu(iobsloc(k,4))-1)*nx+Xu(iobsloc(k,5)) |
361 |
modelloc(k,4)=Yu(iobsloc(k,4)) |
362 |
modelloc(k,5)=Xu(iobsloc(k,5)) |
363 |
elseif (iobsloc(k,2)==1) then ! v location |
364 |
modelloc(k,1)=(iobsloc(k,2))*nn+(iobsloc(k,3)-1)*nx*ny+ |
365 |
& (Yv(iobsloc(k,4))-1)*nx+Xv(iobsloc(k,5)) |
366 |
modelloc(k,4)=Yv(iobsloc(k,4)) |
367 |
modelloc(k,5)=Xv(iobsloc(k,5)) |
368 |
else ! vob > 1 |
369 |
modelloc(k,1)=(iobsloc(k,2))*nn+(iobsloc(k,3)-1)*nx*ny+ |
370 |
& (Yt(iobsloc(k,4))-1)*nx+Xt(iobsloc(k,5)) |
371 |
modelloc(k,4)=Yt(iobsloc(k,4)) |
372 |
modelloc(k,5)=Xt(iobsloc(k,5)) |
373 |
endif |
374 |
enddo |
375 |
|
376 |
|
377 |
elseif (method==2) then |
378 |
|
379 |
c the following get only the cell centered values since this is coming from |
380 |
c PIV |
381 |
do j=1,ydata |
382 |
do i=1,xdata |
383 |
theta(i,j)= |
384 |
& -atan2(real(i)-xorigin-0.5,real(j)-yorigin-0.5) |
385 |
c & -atan2(real(i)-xorigin-0.5,real(j)-yorigin-0.5)+pi |
386 |
thetamod(i,j)=1.0+((real(theta(i,j))+pi)/(2.0*pi)* |
387 |
& real(nx-1)) |
388 |
rho(i,j)=sqrt((real(i)-xorigin-0.5)**2 + |
389 |
& (real(j)-yorigin-0.5)**2) |
390 |
c rhomod(i,j)=1.0+rho(i,j)*(real(ny)-1)/(real(xdata)/2) |
391 |
c rhomod(i,j)=rho(i,j)*(44.0/(37/2)) |
392 |
c rhomod(i,j)=rho(i,j)*(30.0/(38.0/2.0)) |
393 |
rhomod(i,j)=rho(i,j)*radrat |
394 |
enddo |
395 |
enddo |
396 |
|
397 |
|
398 |
c for each obervation, get nearest neighbor value |
399 |
do k=1,currobs |
400 |
|
401 |
modelloc(k,2)=iobsloc(k,2) |
402 |
modelloc(k,3)=iobsloc(k,3) |
403 |
c vob accessed more than once to make reasoning clear, heh |
404 |
if (iobsloc(k,2)==0) then ! u location |
405 |
X=round(thetamod(iobsloc(k,5),iobsloc(k,4))) |
406 |
Y=ceiling(rhomod(iobsloc(k,5),iobsloc(k,4))) |
407 |
modelloc(k,1)=(iobsloc(k,2))*nn+(iobsloc(k,3)-1)*nx*ny+ |
408 |
& (Y-1)*nx+X |
409 |
modelloc(k,4)=Y |
410 |
modelloc(k,5)=X |
411 |
elseif (iobsloc(k,2)==1) then ! v location |
412 |
X=ceiling(thetamod(iobsloc(k,5),iobsloc(k,4))) |
413 |
Y=round(rhomod(iobsloc(k,5),iobsloc(k,4))) |
414 |
modelloc(k,1)=(iobsloc(k,2))*nn+(iobsloc(k,3)-1)*nx*ny+ |
415 |
& (Y-1)*nx+X |
416 |
modelloc(k,4)=Y |
417 |
modelloc(k,5)=X |
418 |
else ! vob > 1 |
419 |
X=ceiling(thetamod(iobsloc(k,5),iobsloc(k,4))) |
420 |
Y=ceiling(rhomod(iobsloc(k,5),iobsloc(k,4))) |
421 |
modelloc(k,1)=(iobsloc(k,2))*nn+(iobsloc(k,3)-1)*nx*ny+ |
422 |
& (Y-1)*nx+X |
423 |
modelloc(k,4)=Y |
424 |
modelloc(k,5)=X |
425 |
endif |
426 |
|
427 |
c write(*,*) modelloc(k,5), modelloc(k,4), modelloc(k,3), |
428 |
c & modelloc(k,1) |
429 |
|
430 |
enddo |
431 |
c stop |
432 |
else |
433 |
stop 'no valid applyH method specified' |
434 |
endif |
435 |
|
436 |
|
437 |
end subroutine getmodelloc |
438 |
|
439 |
function round(somereal) |
440 |
implicit none |
441 |
integer :: round |
442 |
real, intent(in) :: somereal |
443 |
|
444 |
if (modulo(somereal,1.0) .lt. 0.5) then |
445 |
round=floor(somereal) |
446 |
else |
447 |
round=ceiling(somereal) |
448 |
endif |
449 |
end function round |