/[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.2 - (hide annotations) (download)
Mon Oct 11 16:38:53 2004 UTC (19 years, 7 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint57d_post, checkpoint57b_post, checkpoint57c_pre, checkpoint55j_post, checkpoint56b_post, checkpoint55h_post, checkpoint57e_post, checkpoint56c_post, checkpoint57a_post, checkpoint55g_post, checkpoint55f_post, checkpoint57a_pre, checkpoint55i_post, checkpoint57, checkpoint56, eckpoint57e_pre, checkpoint57c_post, checkpoint55e_post, checkpoint56a_post, checkpoint55d_post
Changes since 1.1: +2 -2 lines
o ECCO specific cost function terms (up-to-date with 1x1 runs)
o ecco_cost_weights is modified to 1x1 runs
o modifs to allow observations to be read in as
  single file or yearly files

1 heimbach 1.2 C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_drifter.F,v 1.1 2003/11/06 22:10:07 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     www(i,j,bi,bj) = cosphi(i,j,bi,bj)
206     enddo
207     enddo
208    
209     i6min=1
210    
211     do j = jmin,jmax-1,2
212     do i = i6min,imax-5,6
213     masktmp(i,j,bi,bj) =
214     & (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    
241     do j = jmin,jmax-1,2
242     do i = i6min,imax-5,6
243     masktmp(i,j,bi,bj) =
244     & (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     & +mask13c(i+4,j+1,bi,bj)+mask13c(i+5,j+1,bi,bj))
250     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    
271     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     CADJ STORE wud(:,:,bi,bj)
290     CADJ & = tapelev_ini_bibj_k, key=ikey, byte=isbyte
291     CADJ STORE wvd(:,:,bi,bj)
292     CADJ & = tapelev_ini_bibj_k, key=ikey, byte=isbyte
293     CADJ STORE u6bar(:,:,bi,bj)
294     CADJ & = tapelev_ini_bibj_k, key=ikey, byte=isbyte
295     CADJ STORE v6bar(:,:,bi,bj)
296     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     fctile_drift = fctile_drift
303     & + (wud(i,j,bi,bj)*www(i,j,bi,bj)*mask6c(i,j,bi,bj)*
304     & (u6bar(i,j,bi,bj) - uddat(i,j,bi,bj))*
305     & (u6bar(i,j,bi,bj) - uddat(i,j,bi,bj)) )
306     & + (wvd(i,j,bi,bj)*www(i,j,bi,bj)*mask6c(i,j,bi,bj)*
307     & (v6bar(i,j,bi,bj) - vddat(i,j,bi,bj))*
308     & (v6bar(i,j,bi,bj) - vddat(i,j,bi,bj)) )
309     enddo
310     enddo
311    
312     fcthread_drift = fcthread_drift + fctile_drift
313     objf_drift(bi,bj) = objf_drift(bi,bj) + fctile_drift
314    
315     #ifdef ECCO_VERBOSE
316     c-- Print cost function for each tile in each thread.
317     write(msgbuf,'(a)') ' '
318     call print_message( msgbuf, standardmessageunit,
319     & SQUEEZE_RIGHT , mythid)
320     write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
321     & ' cost_drifter: irec,bi,bj = ',irec,bi,bj
322     call print_message( msgbuf, standardmessageunit,
323     & SQUEEZE_RIGHT , mythid)
324     write(msgbuf,'(a,d22.15)')
325     & ' cost function (temperature) = ',
326 heimbach 1.2 & fctile_drift
327 heimbach 1.1 call print_message( msgbuf, standardmessageunit,
328     & SQUEEZE_RIGHT , mythid)
329     #endif
330    
331     enddo
332     enddo
333    
334     #ifdef ECCO_VERBOSE
335     c-- Print cost function for all tiles.
336     _GLOBAL_SUM_R8( fcthread_drift , myThid )
337     write(msgbuf,'(a)') ' '
338     call print_message( msgbuf, standardmessageunit,
339     & SQUEEZE_RIGHT , mythid)
340     write(msgbuf,'(a,i8.8)')
341     & ' cost_drift: irec = ',irec
342     call print_message( msgbuf, standardmessageunit,
343     & SQUEEZE_RIGHT , mythid)
344     write(msgbuf,'(a,a,d22.15)')
345     & ' global cost function value',
346     & ' (drifters) = ',fcthread_drift
347     call print_message( msgbuf, standardmessageunit,
348     & SQUEEZE_RIGHT , mythid)
349     write(msgbuf,'(a)') ' '
350     call print_message( msgbuf, standardmessageunit,
351     & SQUEEZE_RIGHT , mythid)
352     #endif
353    
354     #endif
355    
356     return
357     end
358    

  ViewVC Help
Powered by ViewVC 1.1.22