/[MITgcm]/MITgcm/pkg/ecco/cost_obcs_ageos.F
ViewVC logotype

Contents of /MITgcm/pkg/ecco/cost_obcs_ageos.F

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


Revision 1.2 - (show annotations) (download)
Tue Oct 9 00:02:50 2007 UTC (16 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint60, checkpoint61, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59k, checkpoint59j, checkpoint61b, checkpoint61a
Changes since 1.1: +17 -16 lines
add missing cvs $Header:$ or $Name:$

1 C $Header: $
2 C $Name: $
3
4 #include "COST_CPPOPTIONS.h"
5 #ifdef ALLOW_OBCS
6 # include "OBCS_OPTIONS.h"
7 #endif
8
9 subroutine cost_obcs_ageos( myiter, mytime, mythid )
10
11 c ==================================================================
12 c SUBROUTINE cost_obcs_ageos
13 c ==================================================================
14 c
15 c o cost function contribution obc -- Ageostrophic boundary flow.
16 c
17 c started: G. Gebbie gebbie@mit.edu 4-Feb-2003
18 c
19 c warning: masks redundantly assume no gradient of topography at
20 c boundary.
21 c
22 c ==================================================================
23 c SUBROUTINE cost_obcs_ageos
24 c ==================================================================
25
26 implicit none
27
28 c == global variables ==
29
30 #include "EEPARAMS.h"
31 #include "SIZE.h"
32 #include "GRID.h"
33 #include "DYNVARS.h"
34 #include "PARAMS.h"
35 #ifdef ALLOW_OBCS
36 # include "OBCS.h"
37 #endif
38
39 #include "cal.h"
40 #include "ecco_cost.h"
41 #include "ctrl.h"
42 #include "ctrl_dummy.h"
43 #include "optim.h"
44
45 c == routine arguments ==
46
47 integer myiter
48 _RL mytime
49 integer mythid
50
51 c == local variables ==
52
53 _RS one_rs
54 parameter( one_rs = 1. )
55
56 integer bi,bj
57 integer i,j,k
58 integer itlo,ithi
59 integer jtlo,jthi
60 integer jmin,jmax
61 integer imin,imax
62 integer irec
63 integer levmon
64 integer levoff
65 integer iltheta
66 integer ilsalt
67 integer iluvel
68 integer ilvvel
69 integer ip1, jp1
70
71 _RL fctile
72 _RL fcthread
73
74 _RL rholoc (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
75 _RL xzgrdrho(1-olx:snx+olx,Nr,nsx,nsy)
76 _RL yzgrdrho(1-oly:sny+oly,Nr,nsx,nsy)
77 _RL xzdvel1 (1-olx:snx+olx,nr,nsx,nsy)
78 _RL xzdvel2 (1-olx:snx+olx,nr,nsx,nsy)
79 _RL yzdvel1 (1-oly:sny+oly,nr,nsx,nsy)
80 _RL yzdvel2 (1-oly:sny+oly,nr,nsx,nsy)
81 _RL maskxzageos (1-olx:snx+olx,nr,nsx,nsy)
82 _RL maskyzageos (1-oly:sny+oly,nr,nsx,nsy)
83 _RL dummy
84
85 character*(80) fnametheta
86 character*(80) fnamesalt
87 character*(80) fnameuvel
88 character*(80) fnamevvel
89
90 logical doglobalread
91 logical ladinit
92
93 character*(MAX_LEN_MBUF) msgbuf
94
95 c == external functions ==
96
97 integer ilnblnk
98 external ilnblnk
99
100 c == end of interface ==
101
102 jtlo = mybylo(mythid)
103 jthi = mybyhi(mythid)
104 itlo = mybxlo(mythid)
105 ithi = mybxhi(mythid)
106 jmin = 1
107 jmax = sny
108 imin = 1
109 imax = snx
110
111 c-- Read tiled data.
112 doglobalread = .false.
113 ladinit = .false.
114
115 #ifdef OBCS_AGEOS_COST_CONTRIBUTION
116
117 #ifdef ECCO_VERBOSE
118 _BEGIN_MASTER( mythid )
119 write(msgbuf,'(a)') ' '
120 call print_message( msgbuf, standardmessageunit,
121 & SQUEEZE_RIGHT , mythid)
122 write(msgbuf,'(a,i8.8)')
123 & ' cost_obcs_ageos: number of records to process = ',nmonsrec
124 call print_message( msgbuf, standardmessageunit,
125 & SQUEEZE_RIGHT , mythid)
126 write(msgbuf,'(a)') ' '
127 call print_message( msgbuf, standardmessageunit,
128 & SQUEEZE_RIGHT , mythid)
129 _END_MASTER( mythid )
130 #endif
131
132 if (optimcycle .ge. 0) then
133 iltheta = ilnblnk( tbarfile )
134 write(fnametheta(1:80),'(2a,i10.10)')
135 & tbarfile(1:iltheta),'.',optimcycle
136 ilsalt = ilnblnk( sbarfile )
137 write(fnamesalt(1:80),'(2a,i10.10)')
138 & sbarfile(1:ilsalt),'.',optimcycle
139 iluvel = ilnblnk( ubarfile )
140 write(fnameuvel(1:80),'(2a,i10.10)')
141 & ubarfile(1:iluvel),'.',optimcycle
142 ilvvel = ilnblnk( vbarfile )
143 write(fnamevvel(1:80),'(2a,i10.10)')
144 & vbarfile(1:ilvvel),'.',optimcycle
145 endif
146
147 fcthread = 0. _d 0
148 fctile = 0. _d 0
149
150 cgg Code safe: always initialize to zero.
151 do bj = jtlo,jthi
152 do bi = itlo,ithi
153 do k = 1,nr
154 do i = 1-olx,snx+olx
155 maskxzageos(i,k,bi,bj) = 0. _d 0
156 xzdvel1(i,k,bi,bj) = 0. _d 0
157 xzdvel2(i,k,bi,bj) = 0. _d 0
158 xzgrdrho(i,k,bi,bj) = 0. _d 0
159 enddo
160 do j = 1-oly,sny+oly
161 maskyzageos(j,k,bi,bj) = 0. _d 0
162 yzdvel1(j,k,bi,bj) = 0. _d 0
163 yzdvel2(j,k,bi,bj) = 0. _d 0
164 yzgrdrho(j,k,bi,bj) = 0. _d 0
165 enddo
166 enddo
167 enddo
168 enddo
169
170 do bj = jtlo,jthi
171 do bi = itlo,ithi
172 do j = 1-oly,sny+oly
173 do i = 1-olx,snx+olx
174 rholoc(i,j,bi,bj) = 0. _d 0
175 enddo
176 enddo
177 enddo
178 enddo
179
180 c-- Loop over records.
181 do irec = 1,nmonsrec
182
183 c-- Read time averages and the monthly mean data.
184 call active_read_xyz( fnametheta, tbar, irec,
185 & doglobalread, ladinit,
186 & optimcycle, mythid,
187 & xx_tbar_mean_dummy )
188
189 c-- Read time averages and the monthly mean data.
190 call active_read_xyz( fnamesalt, sbar, irec,
191 & doglobalread, ladinit,
192 & optimcycle, mythid,
193 & xx_sbar_mean_dummy )
194
195 c-- Read time averages and the monthly mean data.
196 call active_read_xyz( fnameuvel, ubar, irec,
197 & doglobalread, ladinit,
198 & optimcycle, mythid,
199 & xx_ubar_mean_dummy )
200
201 c-- Read time averages and the monthly mean data.
202 call active_read_xyz( fnamevvel, vbar, irec,
203 & doglobalread, ladinit,
204 & optimcycle, mythid,
205 & xx_vbar_mean_dummy )
206
207 cgg Minor problem : grad T,S needs overlap values.
208 _BARRIER
209 _EXCH_XYZ_R8(tbar,myThid)
210 _EXCH_XYZ_R8(sbar,myThid)
211
212
213 do bj = jtlo,jthi
214 do bi = itlo,ithi
215
216 #ifdef ALLOW_OBCSN_CONTROL
217 jp1 = 0
218
219 cgg Make a mask for the velocity shear comparison.
220 do k = 1,nr-1
221 do i = imin, imax
222 j = ob_jn(i,bi,bj)
223 cgg All these points need to be wet.
224 if (j .eq. 0) then
225 maskxzageos(i,k,bi,bj) = 0.
226 else
227 maskxzageos(i,k,bi,bj) =
228 & hfacC(i,j+jp1,k,bi,bj)*hfacC(i+1,j+jp1,k,bi,bj) *
229 & hfacC(i-1,j+jp1,k,bi,bj)*hfacC(i,j+jp1,k+1,bi,bj)*
230 & hfacC(i-1,j+jp1,k+1,bi,bj)*hfacC(i+1,j+jp1,k+1,bi,bj)*
231 & hfacS(i,j+jp1,k,bi,bj)*hfacS(i,j+jp1,k+1,bi,bj)
232 endif
233 enddo
234 enddo
235
236 do k = 1,nr
237 call find_rho(bi,bj,imin,imax,jmin,jmax,k,k,
238 & tbar,sbar,rholoc,mythid)
239 _EXCH_XY_R8(rholoc , myThid)
240
241 cgg Compute centered difference horizontal gradient on bdy.
242 do i = imin, imax
243 j = ob_jn(i,bi,bj)
244 xzgrdrho(i,k,bi,bj) =
245 & (rholoc(i-1,j+jp1,bi,bj)-rholoc(i+1,j+jp1,bi,bj))
246 & /(2.*dxc(i,j+jp1,bi,bj))
247 enddo
248 enddo
249
250 cgg Compute vertical shear from geostrophy/thermal wind.
251 cgg Above level 4 needs not be geostrophic.
252 cgg Please get rid of the "4" ASAP. Ridiculous.
253 do k = 4,Nr-1
254 do i = imin,imax
255 j = ob_jn(i,bi,bj)
256 xzdvel1(i,k,bi,bj) = vbar(i,j+jp1,k ,bi,bj)
257 & - vbar(i,j+jp1,k+1,bi,bj)
258 xzdvel2(i,k,bi,bj)=((xzgrdrho(i,k,bi,bj)*delz(k)/2.)+
259 & (xzgrdrho(i,k+1,bi,bj)*delz(k+1)/2.))
260 & * gravity / f0 / rhonil
261
262 fctile = fctile + 100*wcurrent(k,bi,bj) *
263 & maskxzageos(i,k,bi,bj)*
264 & (xzdvel2(i,k,bi,bj) - xzdvel1(i,k,bi,bj))*
265 & (xzdvel2(i,k,bi,bj) - xzdvel1(i,k,bi,bj))
266 if (maskxzageos(i,k,bi,bj) .ne. 0) then
267 cgg print*,'XX shear,geos shear',xzdvel1(i,k,bi,bj),xzdvel2(i,k,bi,bj)
268 cgg print*,'i,j,k,fctile N',i,j,k,fctile
269 endif
270 enddo
271 enddo
272 c-- End of loop over layers.
273 #endif /* ALLOW_OBCSN_CONTROL */
274
275 #ifdef ALLOW_OBCSS_CONTROL
276 jp1 = 1
277
278 cgg Make a mask for the velocity shear comparison.
279 do k = 1,nr-1
280 do i = imin, imax
281 j = ob_js(i,bi,bj)
282 if (j .eq. 0) then
283 maskxzageos(i,k,bi,bj) = 0.
284 else
285 cgg All these points need to be wet.
286 maskxzageos(i,k,bi,bj) =
287 & hfacC(i,j+jp1,k,bi,bj)*hfacC(i+1,j+jp1,k,bi,bj) *
288 & hfacC(i-1,j+jp1,k,bi,bj)*hfacC(i,j+jp1,k+1,bi,bj)*
289 & hfacC(i-1,j+jp1,k+1,bi,bj)*hfacC(i+1,j+jp1,k+1,bi,bj)*
290 & hfacS(i,j+jp1,k,bi,bj)*hfacS(i,j+jp1,k+1,bi,bj)
291 endif
292 enddo
293 enddo
294
295 do k = 1,nr
296
297 call find_rho(bi,bj,imin,imax,jmin,jmax,k,k,
298 & tbar,sbar,rholoc,mythid)
299 _EXCH_XY_R8(rholoc , myThid)
300
301 cgg Compute centered difference horizontal gradient on bdy.
302 do i = imin, imax
303 j = ob_js(i,bi,bj)
304 xzgrdrho(i,k,bi,bj) =
305 & (rholoc(i-1,j+jp1,bi,bj)-rholoc(i+1,j+jp1,bi,bj))
306 & /(2.*dxc(i,j+jp1,bi,bj))
307 enddo
308 enddo
309
310 cgg Compute vertical shear from geostrophy/thermal wind.
311 do k = 4,Nr-1
312 do i = imin,imax
313 j = ob_js(i,bi,bj)
314 cgg Retrieve the model vertical shear.
315 xzdvel1(i,k,bi,bj) = vbar(i,j+jp1,k ,bi,bj)
316 & - vbar(i,j+jp1,k+1,bi,bj)
317
318 cgg Compute vertical shear from geostrophy/thermal wind.
319 xzdvel2(i,k,bi,bj) =((xzgrdrho(i,k,bi,bj)*delz(k)/2.)+
320 & (xzgrdrho(i,k+1,bi,bj)*delz(k+1)/2.))
321 & * gravity / f0 /rhonil
322
323 cgg Make a comparison.
324 fctile = fctile + 100*wcurrent(k,bi,bj) *
325 & maskxzageos(i,k,bi,bj)*
326 & (xzdvel2(i,k,bi,bj) - xzdvel1(i,k,bi,bj))*
327 & (xzdvel2(i,k,bi,bj) - xzdvel1(i,k,bi,bj))
328 cgg print*,'fctile S',fctile
329 enddo
330 enddo
331 c-- End of loop over layers.
332 #endif
333
334 #ifdef ALLOW_OBCSW_CONTROL
335 ip1 = 1
336
337 cgg Make a mask for the velocity shear comparison.
338 do k = 1,nr-1
339 do j = jmin, jmax
340 i = ob_iw(j,bi,bj)
341 cgg All these points need to be wet.
342 if (i.eq. 0) then
343 maskyzageos(j,k,bi,bj) = 0.
344 else
345 maskyzageos(j,k,bi,bj) =
346 & hfacC(i+ip1,j,k,bi,bj)*hfacC(i+ip1,j+1,k,bi,bj) *
347 & hfacC(i+ip1,j-1,k,bi,bj)*hfacC(i+ip1,j,k+1,bi,bj)*
348 & hfacC(i+ip1,j-1,k+1,bi,bj)*hfacC(i+ip1,j+1,k+1,bi,bj)*
349 & hfacW(i+ip1,j,k,bi,bj)*hfacW(i+ip1,j,k+1,bi,bj)
350 endif
351 enddo
352 enddo
353
354 do k = 1,nr
355
356 call find_rho(bi,bj,imin,imax,jmin,jmax,k,k,
357 & tbar,sbar,rholoc,mythid)
358 _EXCH_XY_R8(rholoc , myThid)
359
360 cgg Compute centered difference horizontal gradient on bdy.
361 do j = jmin, jmax
362 i = ob_iw(j,bi,bj)
363 cgg Negative sign due to geostrophy.
364 yzgrdrho(j,k,bi,bj) =
365 & (rholoc(i+ip1,j+1,bi,bj)-rholoc(i+ip1,j-1,bi,bj))
366 & /(2.*dyc(i+ip1,j,bi,bj))
367 enddo
368 enddo
369
370 cgg Compute vertical shear from geostrophy/thermal wind.
371 do k = 4,Nr-1
372 do j = jmin,jmax
373 i = ob_iw(j,bi,bj)
374 cgg Retrieve the model vertical shear.
375 yzdvel1(j,k,bi,bj) = ubar(i+ip1,j,k ,bi,bj)
376 & - ubar(i+ip1,j,k+1,bi,bj)
377
378 cgg Compute vertical shear from geostrophy/thermal wind.
379 yzdvel2(j,k,bi,bj) =((yzgrdrho(j,k ,bi,bj)*delz(k)/2.)+
380 & (yzgrdrho(j,k+1,bi,bj)*delz(k+1)/2.))
381 & * gravity / f0 / rhonil
382
383 cgg Make a comparison.
384 fctile = fctile + 100*wcurrent(k,bi,bj) *
385 & maskyzageos(j,k,bi,bj) *
386 & (yzdvel2(j,k,bi,bj) - yzdvel1(j,k,bi,bj))*
387 & (yzdvel2(j,k,bi,bj) - yzdvel1(j,k,bi,bj))
388 enddo
389 enddo
390 c-- End of loop over layers.
391 #endif
392
393 #ifdef ALLOW_OBCSE_CONTROL
394 ip1 = 0
395
396 cgg Make a mask for the velocity shear comparison.
397 do k = 1,nr-1
398 do j = jmin, jmax
399 i = ob_ie(j,bi,bj)
400 if (i.eq.0) then
401 maskyzageos(j,k,bi,bj) =0.
402 else
403 cgg All these points need to be wet.
404 maskyzageos(j,k,bi,bj) =
405 & hfacC(i+ip1,j,k,bi,bj)*hfacC(i+ip1,j+1,k,bi,bj) *
406 & hfacC(i+ip1,j-1,k,bi,bj)*hfacC(i+ip1,j,k+1,bi,bj)*
407 & hfacC(i+ip1,j-1,k+1,bi,bj)*hfacC(i+ip1,j+1,k+1,bi,bj)*
408 & hfacW(i+ip1,j,k,bi,bj)*hfacW(i+ip1,j,k+1,bi,bj)
409 endif
410 enddo
411 enddo
412
413 do k = 1,nr
414
415 call find_rho(bi,bj,imin,imax,jmin,jmax,k,k,
416 & tbar,sbar,rholoc,mythid)
417 _EXCH_XY_R8(rholoc , myThid)
418
419 cgg Compute centered difference horizontal gradient on bdy.
420 do j = jmin, jmax
421 i = ob_ie(j,bi,bj)
422 cgg Negative sign due to geostrophy.
423 yzgrdrho(j,k,bi,bj) =
424 & (rholoc(i+ip1,,j+1,bi,bj)-rholoc(i+ip1,j-1,bi,bj))
425 & /(2.*dyc(i+ip1,j,bi,bj))
426 enddo
427 enddo
428
429 cgg Compute vertical shear from geostrophy/thermal wind.
430 do k = 4,Nr-1
431 do j = jmin,jmax
432 i = ob_ie(j,bi,bj)
433 cgg Retrieve the model vertical shear.
434 yzdvel1(j,k,bi,bj) = ubar(i+ip1,j,k ,bi,bj)
435 & - ubar(i+ip1,j,k+1,bi,bj)
436
437 cgg Compute vertical shear from geostrophy/thermal wind.
438 yzdvel2(j,k,bi,bj) =((yzgrdrho(j,k ,bi,bj)*delz(k)/2.)+
439 & (yzgrdrho(j,k+1,bi,bj)*delz(k+1)/2.))
440 & * gravity / f0 /rhonil
441
442 cgg Make a comparison.
443 fctile = fctile + 100*wcurrent(k,bi,bj) *
444 & maskyzageos(j,k,bi,bj) *
445 & (yzdvel2(j,k,bi,bj) - yzdvel1(j,k,bi,bj))*
446 & (yzdvel2(j,k,bi,bj) - yzdvel1(j,k,bi,bj))
447
448 enddo
449 enddo
450 c-- End of loop over layers.
451 #endif
452
453 fcthread = fcthread + fctile
454 objf_ageos(bi,bj) = objf_ageos(bi,bj) + fctile
455
456 #ifdef ECCO_VERBOSE
457 c-- Print cost function for each tile in each thread.
458 write(msgbuf,'(a)') ' '
459 call print_message( msgbuf, standardmessageunit,
460 & SQUEEZE_RIGHT , mythid)
461 write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
462 & ' cost_Theta: irec,bi,bj = ',irec,bi,bj
463 call print_message( msgbuf, standardmessageunit,
464 & SQUEEZE_RIGHT , mythid)
465 write(msgbuf,'(a,d22.15)')
466 & ' cost function (temperature) = ',
467 & fctile
468 call print_message( msgbuf, standardmessageunit,
469 & SQUEEZE_RIGHT , mythid)
470 write(msgbuf,'(a)') ' '
471 call print_message( msgbuf, standardmessageunit,
472 & SQUEEZE_RIGHT , mythid)
473 #endif
474
475 enddo
476 enddo
477
478 #ifdef ECCO_VERBOSE
479 c-- Print cost function for all tiles.
480 _GLOBAL_SUM_R8( fcthread , myThid )
481 write(msgbuf,'(a)') ' '
482 call print_message( msgbuf, standardmessageunit,
483 & SQUEEZE_RIGHT , mythid)
484 write(msgbuf,'(a,i8.8)')
485 & ' cost_Theta: irec = ',irec
486 call print_message( msgbuf, standardmessageunit,
487 & SQUEEZE_RIGHT , mythid)
488 write(msgbuf,'(a,a,d22.15)')
489 & ' global cost function value',
490 & ' (temperature) = ',fcthread
491 call print_message( msgbuf, standardmessageunit,
492 & SQUEEZE_RIGHT , mythid)
493 write(msgbuf,'(a)') ' '
494 call print_message( msgbuf, standardmessageunit,
495 & SQUEEZE_RIGHT , mythid)
496 #endif
497
498 enddo
499 c-- End of loop over records.
500
501 #endif
502
503 return
504 end
505
506
507

  ViewVC Help
Powered by ViewVC 1.1.22