/[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.6 - (hide annotations) (download)
Fri Sep 24 23:24:09 2010 UTC (13 years, 8 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint63p, checkpoint63q, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint63, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x
Changes since 1.5: +5 -1 lines
Changes to accomodate ALLOW_AUTODIFF_WHTAPEIO.

1 gforget 1.6 C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_drifter.F,v 1.5 2009/04/28 18:13:28 jmc Exp $
2 jmc 1.4 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 gforget 1.6 #ifndef ALLOW_AUTODIFF_WHTAPEIO
60    
61 heimbach 1.1 c == local variables ==
62    
63     integer bi,bj
64     integer i,j,k
65     integer itlo,ithi
66     integer jtlo,jthi
67     integer jmin,jmax
68     integer imin,imax
69     integer i6min,i6max
70     integer iglomin
71     integer irec
72     integer ilu
73    
74     _RL fctile_drift
75     _RL fcthread_drift
76     _RL www (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
77     _RL wud (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
78     _RL wvd (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
79     _RL uddat (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
80     _RL u6bar (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
81     _RL vddat (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
82     _RL v6bar (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
83     _RL udmod (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
84     _RL vdmod (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
85     _RL mask13c(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
86     _RL mask6c (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
87     _RL masktmp(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
88    
89     character*(80) fnameud
90     character*(80) fnamevd
91    
92     logical doglobalread
93     logical ladinit
94    
95     character*(MAX_LEN_MBUF) msgbuf
96    
97     c == external functions ==
98    
99     integer ilnblnk
100     external ilnblnk
101    
102     c == end of interface ==
103    
104     jtlo = mybylo(mythid)
105     jthi = mybyhi(mythid)
106     itlo = mybxlo(mythid)
107     ithi = mybxhi(mythid)
108     jmin = 1
109     jmax = sny
110     imin = 1
111     imax = snx
112    
113     c-- Read tiled data.
114     doglobalread = .false.
115     ladinit = .false.
116    
117     #ifdef ALLOW_DRIFTER_COST_CONTRIBUTION
118    
119     if (optimcycle .ge. 0) then
120     ilu = ilnblnk( ubarfile )
121     write(fnameud(1:80),'(2a,i10.10)')
122     & ubarfile(1:ilu),'.',optimcycle
123     ilu = ilnblnk( vbarfile )
124     write(fnamevd(1:80),'(2a,i10.10)')
125     & vbarfile(1:ilu),'.',optimcycle
126     endif
127    
128     fcthread_drift = 0. _d 0
129    
130     do bj = jtlo,jthi
131     do bi = itlo,ithi
132    
133     #ifdef ALLOW_AUTODIFF_TAMC
134     act1 = bi - myBxLo(myThid)
135     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
136     act2 = bj - myByLo(myThid)
137     max2 = myByHi(myThid) - myByLo(myThid) + 1
138     act3 = myThid - 1
139     max3 = nTx*nTy
140     act4 = 0
141     ikey = (act1 + 1) + act2*max1
142     & + act3*max1*max2
143     & + act4*max1*max2*max3
144     #endif /* ALLOW_AUTODIFF_TAMC */
145    
146     k = 2
147     do irec = 1,nmonsrec
148    
149     c-- Read time averages and the monthly mean data.
150     call active_read_xyz( fnameud, ubar, irec,
151     & doglobalread, ladinit,
152     & optimcycle, mythid,
153     & xx_ubar_mean_dummy )
154    
155     call active_read_xyz( fnamevd, vbar, irec,
156     & doglobalread, ladinit,
157     & optimcycle, mythid,
158     & xx_vbar_mean_dummy )
159    
160     do j = jmin,jmax
161     do i = imin,imax
162     if(irec.eq.1)then
163     udmod(i,j,bi,bj)=ubar(i,j,k,bi,bj)
164     vdmod(i,j,bi,bj)=vbar(i,j,k,bi,bj)
165     elseif(irec.eq.nmonsrec)then
166     udmod(i,j,bi,bj)=udmod(i,j,bi,bj)/float(nmonsrec)
167     vdmod(i,j,bi,bj)=vdmod(i,j,bi,bj)/float(nmonsrec)
168     else
169     udmod(i,j,bi,bj)=udmod(i,j,bi,bj)+ubar(i,j,k,bi,bj)
170     vdmod(i,j,bi,bj)=vdmod(i,j,bi,bj)+vbar(i,j,k,bi,bj)
171     endif
172     enddo
173     enddo
174     enddo
175    
176     c-- Read drifter data
177     call mdsreadfield( udriftfile, 32, 'RL', 1, udriftdat, 1,
178     & mythid)
179     call mdsreadfield( vdriftfile, 32, 'RL', 1, vdriftdat, 1,
180     & mythid)
181     c-- Read error data
182     call mdsreadfield( udrifterrfile, 32, 'RL', 1, wudrift, 1,
183     & mythid)
184     call mdsreadfield( vdrifterrfile, 32, 'RL', 1, wvdrift, 1,
185     & mythid)
186    
187     fctile_drift = 0. _d 0
188    
189     c-- Calculate mask for tracer cells (0 => land, 1 => water)
190     do j = jmin,jmax
191     do i = imin,imax
192     mask13c(i,j,bi,bj) = 1. _d 0
193     if (_hFacC(i,j,k,bi,bj) .eq. 0.)
194     & mask13c(i,j,bi,bj) = 0. _d 0
195    
196     cph(
197     cph print *, 'WARNING: SPECIFIC SETUP FOR ECCO'
198     cph below statement could be replaced by following
199     cph to make it independnet of Nr:
200     cph
201     cph if ( rC(K) .GT. -1000. ) then
202     cph)
203     c set mask13c=0 in areas shallower than 1000m
204     if (_hFacC(i,j,13,bi,bj) .eq. 0.) then
205     mask13c(i,j,bi,bj) = 0. _d 0
206     endif
207    
208     enddo
209     enddo
210    
211     i6min=1
212    
213     do j = jmin,jmax-1,2
214     do i = i6min,imax-5,6
215 jmc 1.4 masktmp(i,j,bi,bj) =
216 heimbach 1.1 & (mask13c(i,j,bi,bj)+mask13c(i+1,j,bi,bj)
217     & +mask13c(i+2,j,bi,bj)+mask13c(i+3,j,bi,bj)
218     & +mask13c(i+4,j,bi,bj)+mask13c(i+5,j,bi,bj)
219     & +mask13c(i,j+1,bi,bj)+mask13c(i+1,j+1,bi,bj)
220     & +mask13c(i+2,j+1,bi,bj)+mask13c(i+3,j+1,bi,bj)
221     & +mask13c(i+4,j+1,bi,bj)+mask13c(i+5,j+1,bi,bj))
222     if ( masktmp(i,j,bi,bj) .eq. 0.0 ) then
223     u6bar(i,j,bi,bj) = 0.0
224     else
225     u6bar(i,j,bi,bj) = (
226     & udmod(i,j,bi,bj)*mask13c(i,j,bi,bj)
227     & + udmod(i+1,j,bi,bj)*mask13c(i+1,j,bi,bj)
228     & + udmod(i+2,j,bi,bj)*mask13c(i+2,j,bi,bj)
229     & + udmod(i+3,j,bi,bj)*mask13c(i+3,j,bi,bj)
230     & + udmod(i+4,j,bi,bj)*mask13c(i+4,j,bi,bj)
231     & + udmod(i+5,j,bi,bj)*mask13c(i+5,j,bi,bj)
232     & + udmod(i,j+1,bi,bj)*mask13c(i,j+1,bi,bj)
233     & + udmod(i+1,j+1,bi,bj)*mask13c(i+1,j+1,bi,bj)
234     & + udmod(i+2,j+1,bi,bj)*mask13c(i+2,j+1,bi,bj)
235     & + udmod(i+3,j+1,bi,bj)*mask13c(i+3,j+1,bi,bj)
236     & + udmod(i+4,j+1,bi,bj)*mask13c(i+4,j+1,bi,bj)
237     & + udmod(i+5,j+1,bi,bj)*mask13c(i+5,j+1,bi,bj) )
238     & / ( masktmp(i,j,bi,bj) )
239     endif
240     enddo
241     enddo
242 jmc 1.4
243 heimbach 1.1 do j = jmin,jmax-1,2
244     do i = i6min,imax-5,6
245 jmc 1.4 masktmp(i,j,bi,bj) =
246 heimbach 1.1 & (mask13c(i,j,bi,bj)+mask13c(i+1,j,bi,bj)
247     & +mask13c(i+2,j,bi,bj)+mask13c(i+3,j,bi,bj)
248     & +mask13c(i+4,j,bi,bj)+mask13c(i+5,j,bi,bj)
249     & +mask13c(i,j+1,bi,bj)+mask13c(i+1,j+1,bi,bj)
250     & +mask13c(i+2,j+1,bi,bj)+mask13c(i+3,j+1,bi,bj)
251 jmc 1.4 & +mask13c(i+4,j+1,bi,bj)+mask13c(i+5,j+1,bi,bj))
252 heimbach 1.1 if ( masktmp(i,j,bi,bj) .eq.0.0 ) then
253     v6bar(i,j,bi,bj) = 0.0
254     else
255     v6bar(i,j,bi,bj) = (
256     & vdmod(i,j,bi,bj)*mask13c(i,j,bi,bj)
257     & + vdmod(i+1,j,bi,bj)*mask13c(i+1,j,bi,bj)
258     & + vdmod(i+2,j,bi,bj)*mask13c(i+2,j,bi,bj)
259     & + vdmod(i+3,j,bi,bj)*mask13c(i+3,j,bi,bj)
260     & + vdmod(i+4,j,bi,bj)*mask13c(i+4,j,bi,bj)
261     & + vdmod(i+5,j,bi,bj)*mask13c(i+5,j,bi,bj)
262     & + vdmod(i,j+1,bi,bj)*mask13c(i,j+1,bi,bj)
263     & + vdmod(i+1,j+1,bi,bj)*mask13c(i+1,j+1,bi,bj)
264     & + vdmod(i+2,j+1,bi,bj)*mask13c(i+2,j+1,bi,bj)
265     & + vdmod(i+3,j+1,bi,bj)*mask13c(i+3,j+1,bi,bj)
266     & + vdmod(i+4,j+1,bi,bj)*mask13c(i+4,j+1,bi,bj)
267     & + vdmod(i+5,j+1,bi,bj)*mask13c(i+5,j+1,bi,bj) )
268     & / ( masktmp(i,j,bi,bj) )
269     endif
270     enddo
271     enddo
272 jmc 1.4
273 heimbach 1.1 do j = jmin,jmax-1,2
274     do i = i6min,imax-5, 6
275     c-- change unit from cm/s to m/s
276     uddat(i,j,bi,bj) = 0.01*udriftdat(i,j,bi,bj)
277     vddat(i,j,bi,bj) = 0.01*vdriftdat(i,j,bi,bj)
278     c-- 5 cm/s lower limit
279     wud(i,j,bi,bj) = 1e4*max(wudrift(i,j,bi,bj),5.D0)**(-2)
280     wvd(i,j,bi,bj) = 1e4*max(wvdrift(i,j,bi,bj),5.D0)**(-2)
281     c wud(i,j,bi,bj) = 1.0
282     c wvd(i,j,bi,bj) = 1.0
283     mask6c(i,j,bi,bj) = 1.0
284     if ( uddat(i,j,bi,bj).eq.0.0) mask6c(i,j,bi,bj)=0.0
285     if ( abs(uddat(i,j,bi,bj)).gt.900) mask6c(i,j,bi,bj)=0.0
286     if ( vddat(i,j,bi,bj).eq.0.0) mask6c(i,j,bi,bj)=0.0
287     if ( abs(vddat(i,j,bi,bj)).gt.900) mask6c(i,j,bi,bj)=0.0
288     enddo
289     enddo
290    
291 jmc 1.4 CADJ STORE wud(:,:,bi,bj)
292 heimbach 1.1 CADJ & = tapelev_ini_bibj_k, key=ikey, byte=isbyte
293 jmc 1.4 CADJ STORE wvd(:,:,bi,bj)
294 heimbach 1.1 CADJ & = tapelev_ini_bibj_k, key=ikey, byte=isbyte
295 jmc 1.4 CADJ STORE u6bar(:,:,bi,bj)
296 heimbach 1.1 CADJ & = tapelev_ini_bibj_k, key=ikey, byte=isbyte
297 jmc 1.4 CADJ STORE v6bar(:,:,bi,bj)
298 heimbach 1.1 CADJ & = tapelev_ini_bibj_k, key=ikey, byte=isbyte
299    
300     c-- Compute model data misfit and cost function term for
301     c drifters.
302     do j = jmin,jmax-1,2
303     do i = i6min,imax-5, 6
304 jmc 1.4 fctile_drift = fctile_drift
305 heimbach 1.3 & + (wud(i,j,bi,bj)*cosphi(i,j,bi,bj)*
306     & mask6c(i,j,bi,bj)*
307 heimbach 1.1 & (u6bar(i,j,bi,bj) - uddat(i,j,bi,bj))*
308 jmc 1.4 & (u6bar(i,j,bi,bj) - uddat(i,j,bi,bj)) )
309 heimbach 1.3 & + (wvd(i,j,bi,bj)*cosphi(i,j,bi,bj)*
310     & mask6c(i,j,bi,bj)*
311 heimbach 1.1 & (v6bar(i,j,bi,bj) - vddat(i,j,bi,bj))*
312     & (v6bar(i,j,bi,bj) - vddat(i,j,bi,bj)) )
313 heimbach 1.3 if ( cosphi(i,j,bi,bj)*mask6c(i,j,bi,bj) .ne. 0. ) then
314     if ( wud(i,j,bi,bj) .ne. 0. )
315     & num_drift(bi,bj) = num_drift(bi,bj) + 1. _d 0
316     if ( wvd(i,j,bi,bj) .ne. 0. )
317     & num_drift(bi,bj) = num_drift(bi,bj) + 1. _d 0
318     endif
319 heimbach 1.1 enddo
320     enddo
321    
322     fcthread_drift = fcthread_drift + fctile_drift
323     objf_drift(bi,bj) = objf_drift(bi,bj) + fctile_drift
324    
325     #ifdef ECCO_VERBOSE
326     c-- Print cost function for each tile in each thread.
327     write(msgbuf,'(a)') ' '
328     call print_message( msgbuf, standardmessageunit,
329     & SQUEEZE_RIGHT , mythid)
330     write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
331     & ' cost_drifter: irec,bi,bj = ',irec,bi,bj
332     call print_message( msgbuf, standardmessageunit,
333     & SQUEEZE_RIGHT , mythid)
334     write(msgbuf,'(a,d22.15)')
335     & ' cost function (temperature) = ',
336 heimbach 1.2 & fctile_drift
337 heimbach 1.1 call print_message( msgbuf, standardmessageunit,
338     & SQUEEZE_RIGHT , mythid)
339     #endif
340    
341     enddo
342     enddo
343    
344     #ifdef ECCO_VERBOSE
345     c-- Print cost function for all tiles.
346 jmc 1.5 _GLOBAL_SUM_RL( fcthread_drift , myThid )
347 heimbach 1.1 write(msgbuf,'(a)') ' '
348     call print_message( msgbuf, standardmessageunit,
349     & SQUEEZE_RIGHT , mythid)
350     write(msgbuf,'(a,i8.8)')
351     & ' cost_drift: irec = ',irec
352     call print_message( msgbuf, standardmessageunit,
353     & SQUEEZE_RIGHT , mythid)
354     write(msgbuf,'(a,a,d22.15)')
355     & ' global cost function value',
356     & ' (drifters) = ',fcthread_drift
357     call print_message( msgbuf, standardmessageunit,
358     & SQUEEZE_RIGHT , mythid)
359     write(msgbuf,'(a)') ' '
360     call print_message( msgbuf, standardmessageunit,
361     & SQUEEZE_RIGHT , mythid)
362     #endif
363    
364 gforget 1.6 #endif /* ALLOW_AUTODIFF_WHTAPEIO */
365    
366 heimbach 1.1 #endif
367    
368     return
369     end
370    

  ViewVC Help
Powered by ViewVC 1.1.22