/[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.3 - (hide annotations) (download)
Mon Mar 28 23:49:49 2005 UTC (19 years, 2 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint58l_post, checkpoint57t_post, checkpoint57o_post, checkpoint58e_post, checkpoint57v_post, checkpoint58u_post, checkpoint57f_post, checkpoint57s_post, checkpoint57k_post, checkpoint57g_post, checkpoint58r_post, checkpoint57i_post, checkpoint57y_post, checkpoint58g_post, checkpoint57x_post, checkpoint57m_post, checkpoint58n_post, checkpoint58x_post, checkpoint57g_pre, checkpoint58t_post, checkpoint58h_post, checkpoint58w_post, checkpoint57y_pre, checkpoint57f_pre, checkpoint58q_post, checkpoint58j_post, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59h, checkpoint57r_post, checkpoint59, checkpoint58, checkpoint57h_done, checkpoint58f_post, checkpoint57n_post, checkpoint58d_post, checkpoint58c_post, checkpoint57w_post, checkpoint57p_post, checkpint57u_post, checkpoint58a_post, checkpoint58i_post, checkpoint57q_post, checkpoint58o_post, checkpoint57z_post, checkpoint58y_post, checkpoint58k_post, checkpoint58v_post, checkpoint58s_post, checkpoint58p_post, checkpoint57j_post, checkpoint58b_post, checkpoint57h_pre, checkpoint58m_post, checkpoint57l_post, checkpoint57h_post
Changes since 1.2: +11 -4 lines
Adding counters to cost terms.

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

  ViewVC Help
Powered by ViewVC 1.1.22