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

Annotation of /MITgcm/pkg/ecco/cost_drifter.F

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


Revision 1.4 - (hide 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, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61l, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i
Changes since 1.3: +13 -12 lines
add missing cvs $Header:$ or $Name:$

1 jmc 1.4 C $Header: $
2     C $Name: $
3 heimbach 1.1
4     #include "COST_CPPOPTIONS.h"
5    
6    
7     subroutine cost_drifter(
8     I myiter,
9     I mytime,
10     I mythid
11     & )
12    
13     c ==================================================================
14     c SUBROUTINE cost_drifter
15     c ==================================================================
16     c
17     c o Evaluate cost function contribution of temperature.
18     c
19     c started: Christian Eckert eckert@mit.edu 30-Jun-1999
20     c
21     c changed: Christian Eckert eckert@mit.edu 25-Feb-2000
22     c
23     c - Restructured the code in order to create a package
24     c for the MITgcmUV.
25     c
26     c changed: Patrick Heimbach heimbach@mit.edu 27-May-2000
27     c
28     c - set ladinit to .true. to initialise adubar file
29     c
30     c ==================================================================
31     c SUBROUTINE cost_drifter
32     c ==================================================================
33    
34     implicit none
35    
36     c == global variables ==
37    
38     #include "EEPARAMS.h"
39     #include "SIZE.h"
40     #include "GRID.h"
41     #include "DYNVARS.h"
42    
43     #include "cal.h"
44     #include "ecco_cost.h"
45     #include "ctrl.h"
46     #include "ctrl_dummy.h"
47     #include "optim.h"
48     #ifdef ALLOW_AUTODIFF_TAMC
49     # include "tamc.h"
50     # include "tamc_keys.h"
51     #endif /* ALLOW_AUTODIFF_TAMC */
52    
53     c == routine arguments ==
54    
55     integer myiter
56     _RL mytime
57     integer mythid
58    
59     c == local variables ==
60    
61     integer bi,bj
62     integer i,j,k
63     integer itlo,ithi
64     integer jtlo,jthi
65     integer jmin,jmax
66     integer imin,imax
67     integer i6min,i6max
68     integer iglomin
69     integer irec
70     integer ilu
71    
72     _RL fctile_drift
73     _RL fcthread_drift
74     _RL www (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
75     _RL wud (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
76     _RL wvd (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
77     _RL uddat (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
78     _RL u6bar (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
79     _RL vddat (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
80     _RL v6bar (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
81     _RL udmod (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
82     _RL vdmod (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
83     _RL mask13c(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
84     _RL mask6c (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
85     _RL masktmp(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
86    
87     character*(80) fnameud
88     character*(80) fnamevd
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 ALLOW_DRIFTER_COST_CONTRIBUTION
116    
117     if (optimcycle .ge. 0) then
118     ilu = ilnblnk( ubarfile )
119     write(fnameud(1:80),'(2a,i10.10)')
120     & ubarfile(1:ilu),'.',optimcycle
121     ilu = ilnblnk( vbarfile )
122     write(fnamevd(1:80),'(2a,i10.10)')
123     & vbarfile(1:ilu),'.',optimcycle
124     endif
125    
126     fcthread_drift = 0. _d 0
127    
128     do bj = jtlo,jthi
129     do bi = itlo,ithi
130    
131     #ifdef ALLOW_AUTODIFF_TAMC
132     act1 = bi - myBxLo(myThid)
133     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
134     act2 = bj - myByLo(myThid)
135     max2 = myByHi(myThid) - myByLo(myThid) + 1
136     act3 = myThid - 1
137     max3 = nTx*nTy
138     act4 = 0
139     ikey = (act1 + 1) + act2*max1
140     & + act3*max1*max2
141     & + act4*max1*max2*max3
142     #endif /* ALLOW_AUTODIFF_TAMC */
143    
144     k = 2
145     do irec = 1,nmonsrec
146    
147     c-- Read time averages and the monthly mean data.
148     call active_read_xyz( fnameud, ubar, irec,
149     & doglobalread, ladinit,
150     & optimcycle, mythid,
151     & xx_ubar_mean_dummy )
152    
153     call active_read_xyz( fnamevd, vbar, irec,
154     & doglobalread, ladinit,
155     & optimcycle, mythid,
156     & xx_vbar_mean_dummy )
157    
158     do j = jmin,jmax
159     do i = imin,imax
160     if(irec.eq.1)then
161     udmod(i,j,bi,bj)=ubar(i,j,k,bi,bj)
162     vdmod(i,j,bi,bj)=vbar(i,j,k,bi,bj)
163     elseif(irec.eq.nmonsrec)then
164     udmod(i,j,bi,bj)=udmod(i,j,bi,bj)/float(nmonsrec)
165     vdmod(i,j,bi,bj)=vdmod(i,j,bi,bj)/float(nmonsrec)
166     else
167     udmod(i,j,bi,bj)=udmod(i,j,bi,bj)+ubar(i,j,k,bi,bj)
168     vdmod(i,j,bi,bj)=vdmod(i,j,bi,bj)+vbar(i,j,k,bi,bj)
169     endif
170     enddo
171     enddo
172     enddo
173    
174     c-- Read drifter data
175     call mdsreadfield( udriftfile, 32, 'RL', 1, udriftdat, 1,
176     & mythid)
177     call mdsreadfield( vdriftfile, 32, 'RL', 1, vdriftdat, 1,
178     & mythid)
179     c-- Read error data
180     call mdsreadfield( udrifterrfile, 32, 'RL', 1, wudrift, 1,
181     & mythid)
182     call mdsreadfield( vdrifterrfile, 32, 'RL', 1, wvdrift, 1,
183     & mythid)
184    
185     fctile_drift = 0. _d 0
186    
187     c-- Calculate mask for tracer cells (0 => land, 1 => water)
188     do j = jmin,jmax
189     do i = imin,imax
190     mask13c(i,j,bi,bj) = 1. _d 0
191     if (_hFacC(i,j,k,bi,bj) .eq. 0.)
192     & mask13c(i,j,bi,bj) = 0. _d 0
193    
194     cph(
195     cph print *, 'WARNING: SPECIFIC SETUP FOR ECCO'
196     cph below statement could be replaced by following
197     cph to make it independnet of Nr:
198     cph
199     cph if ( rC(K) .GT. -1000. ) then
200     cph)
201     c set mask13c=0 in areas shallower than 1000m
202     if (_hFacC(i,j,13,bi,bj) .eq. 0.) then
203     mask13c(i,j,bi,bj) = 0. _d 0
204     endif
205    
206     enddo
207     enddo
208    
209     i6min=1
210    
211     do j = jmin,jmax-1,2
212     do i = i6min,imax-5,6
213 jmc 1.4 masktmp(i,j,bi,bj) =
214 heimbach 1.1 & (mask13c(i,j,bi,bj)+mask13c(i+1,j,bi,bj)
215     & +mask13c(i+2,j,bi,bj)+mask13c(i+3,j,bi,bj)
216     & +mask13c(i+4,j,bi,bj)+mask13c(i+5,j,bi,bj)
217     & +mask13c(i,j+1,bi,bj)+mask13c(i+1,j+1,bi,bj)
218     & +mask13c(i+2,j+1,bi,bj)+mask13c(i+3,j+1,bi,bj)
219     & +mask13c(i+4,j+1,bi,bj)+mask13c(i+5,j+1,bi,bj))
220     if ( masktmp(i,j,bi,bj) .eq. 0.0 ) then
221     u6bar(i,j,bi,bj) = 0.0
222     else
223     u6bar(i,j,bi,bj) = (
224     & udmod(i,j,bi,bj)*mask13c(i,j,bi,bj)
225     & + udmod(i+1,j,bi,bj)*mask13c(i+1,j,bi,bj)
226     & + udmod(i+2,j,bi,bj)*mask13c(i+2,j,bi,bj)
227     & + udmod(i+3,j,bi,bj)*mask13c(i+3,j,bi,bj)
228     & + udmod(i+4,j,bi,bj)*mask13c(i+4,j,bi,bj)
229     & + udmod(i+5,j,bi,bj)*mask13c(i+5,j,bi,bj)
230     & + udmod(i,j+1,bi,bj)*mask13c(i,j+1,bi,bj)
231     & + udmod(i+1,j+1,bi,bj)*mask13c(i+1,j+1,bi,bj)
232     & + udmod(i+2,j+1,bi,bj)*mask13c(i+2,j+1,bi,bj)
233     & + udmod(i+3,j+1,bi,bj)*mask13c(i+3,j+1,bi,bj)
234     & + udmod(i+4,j+1,bi,bj)*mask13c(i+4,j+1,bi,bj)
235     & + udmod(i+5,j+1,bi,bj)*mask13c(i+5,j+1,bi,bj) )
236     & / ( masktmp(i,j,bi,bj) )
237     endif
238     enddo
239     enddo
240 jmc 1.4
241 heimbach 1.1 do j = jmin,jmax-1,2
242     do i = i6min,imax-5,6
243 jmc 1.4 masktmp(i,j,bi,bj) =
244 heimbach 1.1 & (mask13c(i,j,bi,bj)+mask13c(i+1,j,bi,bj)
245     & +mask13c(i+2,j,bi,bj)+mask13c(i+3,j,bi,bj)
246     & +mask13c(i+4,j,bi,bj)+mask13c(i+5,j,bi,bj)
247     & +mask13c(i,j+1,bi,bj)+mask13c(i+1,j+1,bi,bj)
248     & +mask13c(i+2,j+1,bi,bj)+mask13c(i+3,j+1,bi,bj)
249 jmc 1.4 & +mask13c(i+4,j+1,bi,bj)+mask13c(i+5,j+1,bi,bj))
250 heimbach 1.1 if ( masktmp(i,j,bi,bj) .eq.0.0 ) then
251     v6bar(i,j,bi,bj) = 0.0
252     else
253     v6bar(i,j,bi,bj) = (
254     & vdmod(i,j,bi,bj)*mask13c(i,j,bi,bj)
255     & + vdmod(i+1,j,bi,bj)*mask13c(i+1,j,bi,bj)
256     & + vdmod(i+2,j,bi,bj)*mask13c(i+2,j,bi,bj)
257     & + vdmod(i+3,j,bi,bj)*mask13c(i+3,j,bi,bj)
258     & + vdmod(i+4,j,bi,bj)*mask13c(i+4,j,bi,bj)
259     & + vdmod(i+5,j,bi,bj)*mask13c(i+5,j,bi,bj)
260     & + vdmod(i,j+1,bi,bj)*mask13c(i,j+1,bi,bj)
261     & + vdmod(i+1,j+1,bi,bj)*mask13c(i+1,j+1,bi,bj)
262     & + vdmod(i+2,j+1,bi,bj)*mask13c(i+2,j+1,bi,bj)
263     & + vdmod(i+3,j+1,bi,bj)*mask13c(i+3,j+1,bi,bj)
264     & + vdmod(i+4,j+1,bi,bj)*mask13c(i+4,j+1,bi,bj)
265     & + vdmod(i+5,j+1,bi,bj)*mask13c(i+5,j+1,bi,bj) )
266     & / ( masktmp(i,j,bi,bj) )
267     endif
268     enddo
269     enddo
270 jmc 1.4
271 heimbach 1.1 do j = jmin,jmax-1,2
272     do i = i6min,imax-5, 6
273     c-- change unit from cm/s to m/s
274     uddat(i,j,bi,bj) = 0.01*udriftdat(i,j,bi,bj)
275     vddat(i,j,bi,bj) = 0.01*vdriftdat(i,j,bi,bj)
276     c-- 5 cm/s lower limit
277     wud(i,j,bi,bj) = 1e4*max(wudrift(i,j,bi,bj),5.D0)**(-2)
278     wvd(i,j,bi,bj) = 1e4*max(wvdrift(i,j,bi,bj),5.D0)**(-2)
279     c wud(i,j,bi,bj) = 1.0
280     c wvd(i,j,bi,bj) = 1.0
281     mask6c(i,j,bi,bj) = 1.0
282     if ( uddat(i,j,bi,bj).eq.0.0) mask6c(i,j,bi,bj)=0.0
283     if ( abs(uddat(i,j,bi,bj)).gt.900) mask6c(i,j,bi,bj)=0.0
284     if ( vddat(i,j,bi,bj).eq.0.0) mask6c(i,j,bi,bj)=0.0
285     if ( abs(vddat(i,j,bi,bj)).gt.900) mask6c(i,j,bi,bj)=0.0
286     enddo
287     enddo
288    
289 jmc 1.4 CADJ STORE wud(:,:,bi,bj)
290 heimbach 1.1 CADJ & = tapelev_ini_bibj_k, key=ikey, byte=isbyte
291 jmc 1.4 CADJ STORE wvd(:,:,bi,bj)
292 heimbach 1.1 CADJ & = tapelev_ini_bibj_k, key=ikey, byte=isbyte
293 jmc 1.4 CADJ STORE u6bar(:,:,bi,bj)
294 heimbach 1.1 CADJ & = tapelev_ini_bibj_k, key=ikey, byte=isbyte
295 jmc 1.4 CADJ STORE v6bar(:,:,bi,bj)
296 heimbach 1.1 CADJ & = tapelev_ini_bibj_k, key=ikey, byte=isbyte
297    
298     c-- Compute model data misfit and cost function term for
299     c drifters.
300     do j = jmin,jmax-1,2
301     do i = i6min,imax-5, 6
302 jmc 1.4 fctile_drift = fctile_drift
303 heimbach 1.3 & + (wud(i,j,bi,bj)*cosphi(i,j,bi,bj)*
304     & mask6c(i,j,bi,bj)*
305 heimbach 1.1 & (u6bar(i,j,bi,bj) - uddat(i,j,bi,bj))*
306 jmc 1.4 & (u6bar(i,j,bi,bj) - uddat(i,j,bi,bj)) )
307 heimbach 1.3 & + (wvd(i,j,bi,bj)*cosphi(i,j,bi,bj)*
308     & mask6c(i,j,bi,bj)*
309 heimbach 1.1 & (v6bar(i,j,bi,bj) - vddat(i,j,bi,bj))*
310     & (v6bar(i,j,bi,bj) - vddat(i,j,bi,bj)) )
311 heimbach 1.3 if ( cosphi(i,j,bi,bj)*mask6c(i,j,bi,bj) .ne. 0. ) then
312     if ( wud(i,j,bi,bj) .ne. 0. )
313     & num_drift(bi,bj) = num_drift(bi,bj) + 1. _d 0
314     if ( wvd(i,j,bi,bj) .ne. 0. )
315     & num_drift(bi,bj) = num_drift(bi,bj) + 1. _d 0
316     endif
317 heimbach 1.1 enddo
318     enddo
319    
320     fcthread_drift = fcthread_drift + fctile_drift
321     objf_drift(bi,bj) = objf_drift(bi,bj) + fctile_drift
322    
323     #ifdef ECCO_VERBOSE
324     c-- Print cost function for each tile in each thread.
325     write(msgbuf,'(a)') ' '
326     call print_message( msgbuf, standardmessageunit,
327     & SQUEEZE_RIGHT , mythid)
328     write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
329     & ' cost_drifter: irec,bi,bj = ',irec,bi,bj
330     call print_message( msgbuf, standardmessageunit,
331     & SQUEEZE_RIGHT , mythid)
332     write(msgbuf,'(a,d22.15)')
333     & ' cost function (temperature) = ',
334 heimbach 1.2 & fctile_drift
335 heimbach 1.1 call print_message( msgbuf, standardmessageunit,
336     & SQUEEZE_RIGHT , mythid)
337     #endif
338    
339     enddo
340     enddo
341    
342     #ifdef ECCO_VERBOSE
343     c-- Print cost function for all tiles.
344     _GLOBAL_SUM_R8( fcthread_drift , myThid )
345     write(msgbuf,'(a)') ' '
346     call print_message( msgbuf, standardmessageunit,
347     & SQUEEZE_RIGHT , mythid)
348     write(msgbuf,'(a,i8.8)')
349     & ' cost_drift: irec = ',irec
350     call print_message( msgbuf, standardmessageunit,
351     & SQUEEZE_RIGHT , mythid)
352     write(msgbuf,'(a,a,d22.15)')
353     & ' global cost function value',
354     & ' (drifters) = ',fcthread_drift
355     call print_message( msgbuf, standardmessageunit,
356     & SQUEEZE_RIGHT , mythid)
357     write(msgbuf,'(a)') ' '
358     call print_message( msgbuf, standardmessageunit,
359     & SQUEEZE_RIGHT , mythid)
360     #endif
361    
362     #endif
363    
364     return
365     end
366    

  ViewVC Help
Powered by ViewVC 1.1.22