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

Annotation of /MITgcm/pkg/ecco/cost_driftw.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, 8 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: +1 -20 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/cost/Attic/cost_driftw.F,v 1.1.2.1 2002/04/04 10:58:58 heimbach Exp $
2 heimbach 1.1
3     #include "COST_CPPOPTIONS.h"
4    
5    
6     subroutine cost_driftw(
7     I myiter,
8     I mytime,
9     I mythid
10     & )
11    
12     c ==================================================================
13     c SUBROUTINE cost_Driftw
14     c ==================================================================
15     c
16     c o Evaluate cost function contribution of the w difference
17     c between the first and the last year.
18     c
19     c started: from cost_drift
20     c
21     c Armin Koehl akoehl@ucsd.edu 26-Feb-2002
22     c
23     c ==================================================================
24     c SUBROUTINE cost_Driftw
25     c ==================================================================
26    
27     implicit none
28    
29     c == global variables ==
30    
31     #include "EEPARAMS.h"
32     #include "SIZE.h"
33     #include "GRID.h"
34     #include "DYNVARS.h"
35    
36     #include "cal.h"
37     #include "ecco_cost.h"
38     #include "ctrl.h"
39     #include "ctrl_dummy.h"
40     #include "optim.h"
41    
42     c == routine arguments ==
43    
44     integer myiter
45     _RL mytime
46     integer mythid
47    
48     c == local variables ==
49    
50     _RS one_rs
51     parameter( one_rs = 1. )
52    
53     integer bi,bj
54     integer i,j,k
55     integer itlo,ithi
56     integer jtlo,jthi
57     integer jmin,jmax
58     integer imin,imax
59     integer irec
60     integer ilw
61     integer nf, nl
62    
63     _RL fctilew
64     _RL fcthread_wdrift
65    
66     character*(80) fnamew
67    
68     logical doglobalread
69     logical ladinit
70    
71     character*(MAX_LEN_MBUF) msgbuf
72    
73     c == external functions ==
74    
75     integer ilnblnk
76     external ilnblnk
77    
78     c == end of interface ==
79    
80     jtlo = mybylo(mythid)
81     jthi = mybyhi(mythid)
82     itlo = mybxlo(mythid)
83     ithi = mybxhi(mythid)
84     jmin = 1
85     jmax = sny
86     imin = 1
87     imax = snx
88    
89     c-- Read tiled data.
90     doglobalread = .false.
91     ladinit = .false.
92    
93     #ifdef ALLOW_DRIFTW_COST_CONTRIBUTION
94    
95     if (optimcycle .ge. 0) then
96     ilw = ilnblnk( wbarfile )
97     write(fnamew(1:80),'(2a,i10.10)')
98     & wbarfile(1:ilw),'.',optimcycle
99     endif
100    
101     if (nmonsrec.lt.12) then
102     print*
103     print*,' cost_Driftw: assimilation period smaller'
104     print*,' than 1 year, no drift computed.'
105     print*
106     stop ' ... stopped in cost_Driftw.'
107     endif
108    
109     fcthread_wdrift = 0. _d 0
110    
111     do bj = jtlo,jthi
112     do bi = itlo,ithi
113     do k = 1,nr
114     do j = jmin,jmax
115     do i = imin,imax
116     wfmean(i,j,k,bi,bj) = 0.0
117     wlmean(i,j,k,bi,bj) = 0.0
118     enddo
119     enddo
120     enddo
121     enddo
122     enddo
123    
124     nf = 0
125     nl = 0
126    
127     c-- Loop over records.
128     do irec = 1,12
129    
130     c-- Read time averages and the monthly mean data.
131     call active_read_xyz( fnamew, wbar, irec,
132     & doglobalread, ladinit,
133     & optimcycle, mythid,
134     & xx_wbar_mean_dummy )
135    
136     nf = nf + 1
137     do bj = jtlo,jthi
138     do bi = itlo,ithi
139     do k = 1,nr
140     do j = jmin,jmax
141     do i = imin,imax
142     wfmean(i,j,k,bi,bj) = wfmean(i,j,k,bi,bj) +
143     & wbar(i,j,k,bi,bj)
144     enddo
145     enddo
146     enddo
147     enddo
148     enddo
149    
150     enddo
151    
152     do irec = nmonsrec-12+1, nmonsrec
153    
154     c-- Read time averages and the monthly mean data.
155     call active_read_xyz( fnamew, wbar, irec,
156     & doglobalread, ladinit,
157     & optimcycle, mythid,
158     & xx_wbar_mean_dummy )
159    
160     nl = nl + 1
161    
162     do bj = jtlo,jthi
163     do bi = itlo,ithi
164     do k = 1,nr
165     do j = jmin,jmax
166     do i = imin,imax
167     wlmean(i,j,k,bi,bj) = wlmean(i,j,k,bi,bj) +
168     & wbar(i,j,k,bi,bj)
169     enddo
170     enddo
171     enddo
172     enddo
173     enddo
174    
175     enddo
176    
177    
178     do bj = jtlo,jthi
179     do bi = itlo,ithi
180    
181     c-- Loop over the model layers
182     fctilew = 0. _d 0
183    
184     do k = 1,nr
185    
186     c-- Compute model misfit and cost function term for
187     c the vertical velovity field. The error is 1e-4 m/s.
188     do j = jmin,jmax
189     do i = imin,imax
190     if (_hFacC(i,j,k,bi,bj) .ne. 0.) then
191     fctilew = fctilew +
192     & (2.5e11*cosphi(i,j,bi,bj)*
193     & (wlmean(i,j,k,bi,bj)/nl - wfmean(i,j,k,bi,bj)/nf)*
194     & (wlmean(i,j,k,bi,bj)/nl - wfmean(i,j,k,bi,bj)/nf))
195     endif
196     enddo
197     enddo
198    
199     enddo
200     c-- End of loop over layers.
201    
202     fcthread_wdrift = fcthread_wdrift + fctilew
203     objf_wdrift(bi,bj) = objf_wdrift(bi,bj) + fctilew
204    
205     #ifdef ECCO_VERBOSE
206     c-- Print cost function for each tile in each thread.
207     write(msgbuf,'(a)') ' '
208     call print_message( msgbuf, standardmessageunit,
209     & SQUEEZE_RIGHT , mythid)
210     write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
211     & ' cost_Driftw: irec,bi,bj = ',irec,bi,bj
212     call print_message( msgbuf, standardmessageunit,
213     & SQUEEZE_RIGHT , mythid)
214     write(msgbuf,'(a,d22.15)')
215     & ' cost function (wvel) = ',
216     & fctilet+fctiles
217     call print_message( msgbuf, standardmessageunit,
218     & SQUEEZE_RIGHT , mythid)
219     write(msgbuf,'(a)') ' '
220     call print_message( msgbuf, standardmessageunit,
221     & SQUEEZE_RIGHT , mythid)
222     #endif
223    
224     enddo
225     enddo
226    
227     #ifdef ECCO_VERBOSE
228     c-- Print cost function for all tiles.
229     _GLOBAL_SUM_R8( fcthread_tdrift , myThid )
230     _GLOBAL_SUM_R8( fcthread_tdrifs , myThid )
231     write(msgbuf,'(a)') ' '
232     call print_message( msgbuf, standardmessageunit,
233     & SQUEEZE_RIGHT , mythid)
234     write(msgbuf,'(a,i8.8)')
235     & ' cost_Drift: irec = ',irec
236     call print_message( msgbuf, standardmessageunit,
237     & SQUEEZE_RIGHT , mythid)
238     write(msgbuf,'(a,a,d22.15)')
239     & ' cost function value',
240     & ' (wvel) = ',fcthread_tdrifw
241     call print_message( msgbuf, standardmessageunit,
242     & SQUEEZE_RIGHT , mythid)
243     write(msgbuf,'(a)') ' '
244     call print_message( msgbuf, standardmessageunit,
245     & SQUEEZE_RIGHT , mythid)
246     #endif
247    
248     #else
249     c-- Do not enter the calculation of the drift contribution to
250     c-- the final cost function.
251    
252     _BEGIN_MASTER( mythid )
253     write(msgbuf,'(a)') ' '
254     call print_message( msgbuf, standardmessageunit,
255     & SQUEEZE_RIGHT , mythid)
256     write(msgbuf,'(a,a)')
257     & ' cost_Driftw: no contribution of drift between the first ',
258     & 'and the last year to cost function.'
259     call print_message( msgbuf, standardmessageunit,
260     & SQUEEZE_RIGHT , mythid)
261     write(msgbuf,'(a,a,i9.8)')
262     & ' cost_Driftw: number of records that would have',
263     & ' been processed: ',nmonsrec
264     call print_message( msgbuf, standardmessageunit,
265     & SQUEEZE_RIGHT , mythid)
266     write(msgbuf,'(a)') ' '
267     call print_message( msgbuf, standardmessageunit,
268     & SQUEEZE_RIGHT , mythid)
269     _END_MASTER( mythid )
270     #endif
271    
272     return
273     end
274    

  ViewVC Help
Powered by ViewVC 1.1.22