/[MITgcm]/MITgcm_contrib/osse/EnKF/rwObsState3.F
ViewVC logotype

Contents of /MITgcm_contrib/osse/EnKF/rwObsState3.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Tue May 4 18:19:35 2004 UTC (21 years, 2 months ago) by afe
Branch: MAIN
o EnKF stuff

1 c*** Reads iobsloc
2 subroutine ReadObs(obsnum,mobs,iobsloc)
3
4 implicit none
5
6 integer, intent(in) :: obsnum,mobs
7 integer, intent(out) :: iobsloc(mobs)
8 character*(100) fn
9 real*8 r8seg(mobs)
10 write(fn,'(A,I3.3,A)') '../inits/iobsloc-',obsnum,'.txt'
11 write(*,*) fn
12
13 open(unit=3,file=fn,status='old')
14 read(unit=3,'(1I7)') iobsloc
15 close(unit=3)
16 call flush()
17 end subroutine ReadObs
18
19 c*** Call Model Kludge
20 subroutine Model(mem)
21 implicit none
22 integer, intent(in) :: mem
23 character*(1000) fn
24 write(fn,'(A,I2.2,A)')
25 & 'cd ../',mem,'/assimilate;./mitgcmuv>& LOG'
26 call system(fn)
27 end subroutine Model
28
29 c*** Reads a pickup file...
30
31 subroutine ReadPickup(mem,n,nx,ny,nz,y)
32
33 implicit none
34
35 integer, intent(in) :: mem,n,nx,ny,nz
36 real*8, intent(out) :: y(n)
37 real*8 r8seg(nx)
38 integer i, irec,j,k
39 character*(1024) fn
40 write(fn,'(A,I2.2,A)') '../',mem,'/assimilate/pickup.in'
41 open(unit=3,file=fn,status='old',
42 & access='direct',recl=nx*8)
43
44 ccc Set counter for state
45 i =1
46
47 ccc Reading uVel
48 do k = 1,nz
49 do j = 1,ny
50 irec = (k-1)*ny + j
51 read(3,rec=irec) r8seg
52 call DA_BYTESWAPR8(nx,r8seg)
53 y(i:i+nx-1) = r8seg(1:nx)
54 i = i+nx
55 enddo
56 enddo
57
58
59 ccc Skip Next Two gU and gUnm1
60 ccc Read vVel
61 do k = 1,nz
62 do j = 1,ny
63 irec = (k-1)*ny + j + nz*ny*3
64 read(3,rec=irec) r8seg
65 call DA_BYTESWAPR8(nx,r8seg)
66 y(i:i+nx-1) = r8seg(1:nx)
67 i=i+nx
68 enddo
69 enddo
70
71
72 ccc Skip Next two gV and gVnm1
73 ccc Read wVel
74 do k = 1,nz
75 do j = 1,ny
76 irec = (k-1)*ny + j + nz*ny*6
77 read(3,rec=irec) r8seg
78 call DA_BYTESWAPR8(nx,r8seg)
79 y(i:i+nx-1) = r8seg(1:nx)
80 i=i+nx
81 enddo
82 enddo
83
84
85 ccc Read theta
86 do k = 1,nz
87 do j = 1,ny
88 irec = (k-1)*ny + j + nz*ny*7
89 read(3,rec=irec) r8seg
90 call DA_BYTESWAPR8(nx,r8seg)
91 y(i:i+nx-1) = r8seg(1:nx)
92 i=i+nx
93 enddo
94 enddo
95
96 ccc Now close the file
97 close(unit=3)
98 ccc Open Non-Hydrostatic
99 c write(fn,'(A,I2.2,A)') '../',mem,'/assimilate/pickup_nh.in'
100 c open(unit=3,file=fn,
101 c & status='old',access='direct',recl=nx*8)
102
103 ccc Read pressure Non-hydrostatic
104 c do k = 1,nz
105 c do j = 1,ny
106 c irec = (k-1)*ny + j
107 c read(3,rec=irec) r8seg
108 c call DA_BYTESWAPR8(nx,r8seg)
109 c y(i:i+nx-1) = r8seg(1:nx)
110 c i = i+nx
111 c enddo
112 c enddo
113
114 c close(unit=3)
115 ccc Done Reading State
116 call flush()
117 end subroutine ReadPickup
118
119 ccc Writes a pickup
120 subroutine WritePickup(mem,n,nx,ny,nz,y)
121
122 implicit none
123
124 integer, intent(in) :: mem,n,nx,ny,nz
125 real*8, intent(in) :: y(n)
126 real*8 r8seg(nx)
127 integer i, irec,j,k
128 character*(1000) fni, fno
129 write(fni,'(A,I2.2,A)') '../',mem,'/assimilate/pickup.in'
130 open(unit=3,file=fni,status='old',
131 & access='direct',recl=nx*8)
132 write(fni,'(A,I2.2,A)') '../',mem,'/assimilate/pickup.out'
133 open(unit=4,file=fni,status='unknown',
134 & access='direct',recl=nx*8)
135 ccc Set counter for state
136 i =1
137
138 ccc Writing uVel
139 do k = 1,nz
140 do j = 1,ny
141 irec = (k-1)*ny + j
142 r8seg(1:nx) = y(i:i+nx-1)
143 call DA_BYTESWAPR8(nx,r8seg)
144 write(4,rec=irec) r8seg
145 i = i+nx
146 enddo
147 enddo
148
149 ccc Copy gU and gUnm1
150 do k = 1,2*nz
151 do j = 1,ny
152 irec = (k-1)*ny + j+nz*ny
153 read(3,rec=irec) r8seg
154 write(4,rec=irec) r8seg
155 enddo
156 enddo
157 ccc Writing vVel
158 do k = 1,nz
159 do j = 1,ny
160 irec = (k-1)*ny + j + nz*ny*3
161 r8seg(1:nx) = y(i:i+nx-1)
162 call DA_BYTESWAPR8(nx,r8seg)
163 write(4,rec=irec) r8seg
164 i=i+nx
165 enddo
166 enddo
167
168 ccc copy Next two gV and gVnm1
169 do k = 1,2*nz
170 do j = 1,ny
171 irec = (k-1)*ny + j+nz*ny*4
172 read(3,rec=irec) r8seg
173 write(4,rec=irec) r8seg
174 enddo
175 enddo
176 ccc Write wVel and Theta
177 do k = 1,2*nz
178 do j = 1,ny
179 irec = (k-1)*ny + j + nz*ny*6
180 r8seg(1:nx) = y(i:i+nx-1)
181 call DA_BYTESWAPR8(nx,r8seg)
182 write(4,rec=irec) r8seg
183 i=i+nx
184 enddo
185 enddo
186 ccc Copy Remaining fields
187 do k = 1,5*nz
188 do j = 1,ny
189 irec = (k-1)*ny + j+nz*ny*8
190 read(3,rec=irec) r8seg
191 write(4,rec=irec) r8seg
192 enddo
193 enddo
194 do j = 1,ny
195 irec = j+13*nz*ny
196 read(3,rec=irec) r8seg
197 write(4,rec=irec) r8seg
198 enddo
199
200 ccc Now close the file
201 close(unit=3)
202 close(unit=4)
203
204 ccc Open Non-Hydrostatic
205 write(fni,'(A,I2.2,A)') '../',mem,'/assimilate/pickup_nh.in'
206 open(unit=3,file=fni,
207 & status='old',access='direct',recl=nx*8)
208 write(fni,'(A,I2.2,A)') '../',mem,'/assimilate/pickup_nh.out'
209 open(unit=4,file=fni,status='unknown',
210 & access='direct',recl=nx*8)
211
212 ccc Write pressure Non-hydrostatic from State
213 c do k = 1,nz
214 c do j = 1,ny
215 c irec = (k-1)*ny + j
216 c r8seg(1:nx) = y(i:i+nx-1)
217 c call DA_BYTESWAPR8(nx,r8seg)
218 c write(4,rec=irec) r8seg
219 c i = i+nx
220 c enddo
221 c enddo
222 ccc Write out other fields
223 do k = 1,2*nz
224 do j = 1,ny
225 irec = (k-1)*ny + j
226 read(3,rec=irec) r8seg
227 write(4,rec=irec) r8seg
228 enddo
229 enddo
230
231 c close(unit=3)
232 c close(unit=4)
233 ccc Done Reading State
234 call flush()
235 end subroutine WritePickup
236
237
238
239
240 subroutine loadMask(mask,nx,ny)
241 implicit none
242 integer, intent(in) :: nx, ny
243 real*4, intent(out) :: mask(ny,nx)
244 real*4 r4seg(nx)
245 integer i, irec,j,k
246 character*(1000) fni, fno
247
248 write(fni,'(A)') '../inits/bathy61.bin'
249
250 open(unit=3,file=fni,
251 & status='old',access='direct',recl=nx*4)
252
253 do j = 1,ny
254 irec = j
255 read(3,rec=irec) r4seg
256 call DA_BYTESWAPR4(nx,r4seg)
257 mask(j,1:nx) = r4seg(1:nx)
258 enddo
259 close(unit=3)
260 call flush()
261 end subroutine loadMask
262
263 subroutine DA_BYTESWAPR8( na, arr )
264 C IN:
265 C n integer - Number of 8-byte words in arr
266 C IN/OUT:
267 C arr real*8 - Array declared as real*4(n)
268 C
269 C Created: 05/05/99 adcroft@mit.edu (This is an unfortunate hack!!)
270
271 implicit none
272 C Arguments
273 integer na
274 character*(*) arr
275 C Local
276 integer i
277 character*(1) cc
278 C ------------------------------------------------------------------
279 do i=1,8*na,8
280 cc=arr(i:i)
281 arr(i:i)=arr(i+7:i+7)
282 arr(i+7:i+7)=cc
283 cc=arr(i+1:i+1)
284 arr(i+1:i+1)=arr(i+6:i+6)
285 arr(i+6:i+6)=cc
286 cc=arr(i+2:i+2)
287 arr(i+2:i+2)=arr(i+5:i+5)
288 arr(i+5:i+5)=cc
289 cc=arr(i+3:i+3)
290 arr(i+3:i+3)=arr(i+4:i+4)
291 arr(i+4:i+4)=cc
292 enddo
293 C ------------------------------------------------------------------
294 return
295 end
296 C==========================
297
298
299 C=======================================================================
300 subroutine DA_BYTESWAPR4( n, arr )
301 C IN:
302 C n integer - Number of 4-byte words in arr
303 C IN/OUT:
304 C arr real*4 - Array declared as real*4(n)
305 C
306 C Created: 05/05/99 adcroft@mit.edu (This is an unfortunate hack!!)
307
308 implicit none
309 C Arguments
310 integer n
311 character*(*) arr
312 C Local
313 integer i
314 character*(1) cc
315 C ------------------------------------------------------------------
316 do i=1,4*n,4
317 cc=arr(i:i)
318 arr(i:i)=arr(i+3:i+3)
319 arr(i+3:i+3)=cc
320 cc=arr(i+1:i+1)
321 arr(i+1:i+1)=arr(i+2:i+2)
322 arr(i+2:i+2)=cc
323 enddo
324 C ------------------------------------------------------------------
325 return
326 end
327 C=======================================================================
328
329
330
331
332
333

  ViewVC Help
Powered by ViewVC 1.1.22