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

Annotation of /MITgcm/pkg/ecco/cost_obcss.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: checkpoint58l_post, checkpoint57t_post, checkpoint57o_post, checkpoint58e_post, checkpoint57v_post, checkpoint58u_post, checkpoint58w_post, checkpoint57m_post, checkpoint57s_post, checkpoint57k_post, checkpoint57d_post, checkpoint57g_post, checkpoint57b_post, checkpoint57c_pre, checkpoint58r_post, checkpoint55j_post, checkpoint56b_post, checkpoint57i_post, checkpoint57y_post, checkpoint57e_post, checkpoint55h_post, checkpoint58n_post, checkpoint58x_post, checkpoint57g_pre, checkpoint58t_post, checkpoint58h_post, checkpoint56c_post, checkpoint57y_pre, checkpoint57f_pre, checkpoint57a_post, checkpoint58q_post, checkpoint55g_post, checkpoint58j_post, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint55f_post, checkpoint59c, checkpoint59b, checkpoint59h, checkpoint57r_post, checkpoint59, checkpoint58, checkpoint57a_pre, checkpoint55i_post, checkpoint57, checkpoint56, eckpoint57e_pre, checkpoint57h_done, checkpoint58f_post, checkpoint57x_post, checkpoint57n_post, checkpoint58d_post, checkpoint58c_post, checkpoint57w_post, checkpoint57p_post, checkpint57u_post, checkpoint57f_post, checkpoint58a_post, checkpoint58i_post, checkpoint57q_post, checkpoint58g_post, checkpoint58o_post, checkpoint57z_post, checkpoint57c_post, checkpoint58y_post, checkpoint55e_post, checkpoint58k_post, checkpoint58v_post, checkpoint58s_post, checkpoint58p_post, checkpoint57j_post, checkpoint58b_post, checkpoint57h_pre, checkpoint58m_post, checkpoint57l_post, checkpoint57h_post, checkpoint56a_post, checkpoint55d_post
Changes since 1.1: +136 -125 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.1
2     #include "COST_CPPOPTIONS.h"
3     #ifdef ALLOW_OBCS
4     # include "OBCS_OPTIONS.h"
5     #endif
6    
7     subroutine cost_obcss(
8     I myiter,
9     I mytime,
10     I mythid
11     & )
12    
13     c ==================================================================
14     c SUBROUTINE cost_obcss
15     c ==================================================================
16     c
17     c o cost function contribution obc
18     c
19     c ==================================================================
20     c SUBROUTINE cost_obcss
21     c ==================================================================
22    
23     implicit none
24    
25     c == global variables ==
26    
27     #include "EEPARAMS.h"
28     #include "SIZE.h"
29     #include "PARAMS.h"
30     #include "GRID.h"
31     #include "DYNVARS.h"
32     #ifdef ALLOW_OBCS
33     # include "OBCS.h"
34     #endif
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     integer bi,bj
51     integer i,j,k
52     integer itlo,ithi
53     integer jtlo,jthi
54     integer jmin,jmax
55     integer imin,imax
56     integer irec
57     integer il
58     integer iobcs
59     integer jp1
60    
61     _RL fctile
62     _RL fcthread
63     _RL dummy
64    
65 heimbach 1.2 character*(80) fnametheta
66     character*(80) fnamesalt
67     character*(80) fnameuvel
68     character*(80) fnamevvel
69 heimbach 1.1
70     logical doglobalread
71     logical ladinit
72    
73     #ifdef ECCO_VERBOSE
74     character*(MAX_LEN_MBUF) msgbuf
75     #endif
76    
77     c == external functions ==
78    
79     integer ilnblnk
80     external ilnblnk
81    
82     c == end of interface ==
83    
84     jtlo = mybylo(mythid)
85     jthi = mybyhi(mythid)
86     itlo = mybxlo(mythid)
87     ithi = mybxhi(mythid)
88     jmin = 1
89     jmax = sny
90     imin = 1
91     imax = snx
92    
93     c-- Read tiled data.
94     doglobalread = .false.
95     ladinit = .false.
96    
97     #ifdef ALLOW_OBCSS_COST_CONTRIBUTION
98    
99     jp1 = 1
100     fcthread = 0. _d 0
101    
102 heimbach 1.2 c-- Loop over records.
103     do irec = 1,nmonsrec
104 heimbach 1.1
105 heimbach 1.2 c-- temperature
106     iobcs = 1
107     c-- Read time averages and the monthly mean data.
108     il = ilnblnk( tbarfile )
109     write(fnametheta(1:80),'(2a,i10.10)')
110     & tbarfile(1:il),'.',optimcycle
111     call active_read_xyz( fnametheta, tbar, irec,
112     & doglobalread, ladinit,
113     & optimcycle, mythid,
114     & xx_tbar_mean_dummy )
115 heimbach 1.1
116 heimbach 1.2 call mdsreadfieldxz( OBStFile, readBinaryPrec, 'RS',
117     & nr, OBSt, irec, mythid)
118 heimbach 1.1
119 heimbach 1.2 do bj = jtlo,jthi
120     do bi = itlo,ithi
121     c-- Loop over the model layers
122     fctile = 0. _d 0
123     do k = 1,nr
124     c-- Compute model data misfit and cost function term for
125     c the temperature field.
126     do i = imin,imax
127     j = OB_Js(I,bi,bj)
128     if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
129     fctile = fctile +
130     & wobcss(k,iobcs)*cosphi(i,j,bi,bj)*
131     & (tbar(i,j,k,bi,bj) - OBSt(i,k,bi,bj))*
132     & (tbar(i,j,k,bi,bj) - OBSt(i,k,bi,bj))
133     endif
134     enddo
135     enddo
136     c-- End of loop over layers.
137     fcthread = fcthread + fctile
138     objf_obcss(bi,bj) = objf_obcss(bi,bj) + fctile
139     enddo
140 heimbach 1.1 enddo
141    
142 heimbach 1.2 c-- salt
143     iobcs = 2
144     c-- Read time averages and the monthly mean data.
145     il = ilnblnk( sbarfile )
146     write(fnamesalt(1:80),'(2a,i10.10)')
147     & sbarfile(1:il),'.',optimcycle
148     call active_read_xyz( fnamesalt, sbar, irec,
149     & doglobalread, ladinit,
150     & optimcycle, mythid,
151     & xx_sbar_mean_dummy )
152    
153     call mdsreadfieldxz( OBSsFile, readBinaryPrec, 'RS',
154     & nr, OBSs, irec, mythid)
155    
156     do bj = jtlo,jthi
157     do bi = itlo,ithi
158     c-- Loop over the model layers
159     fctile = 0. _d 0
160     do k = 1,nr
161     c-- Compute model data misfit and cost function term for
162     c the temperature field.
163     do i = imin,imax
164     j = OB_Js(I,bi,bj)
165     if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
166     fctile = fctile +
167     & wobcss(k,iobcs)*cosphi(i,j,bi,bj)*
168     & (sbar(i,j,k,bi,bj) - OBSs(i,k,bi,bj))*
169     & (sbar(i,j,k,bi,bj) - OBSs(i,k,bi,bj))
170     endif
171     enddo
172 heimbach 1.1 enddo
173 heimbach 1.2 c-- End of loop over layers.
174     fcthread = fcthread + fctile
175     objf_obcss(bi,bj) = objf_obcss(bi,bj) + fctile
176 heimbach 1.1 enddo
177 heimbach 1.2 enddo
178 heimbach 1.1
179 heimbach 1.2 c-- uvel
180     iobcs = 3
181     c-- Read time averages and the monthly mean data.
182     il = ilnblnk( ubarfile )
183     write(fnameuvel(1:80),'(2a,i10.10)')
184     & ubarfile(1:il),'.',optimcycle
185     call active_read_xyz( fnameuvel, ubar, irec,
186     & doglobalread, ladinit,
187     & optimcycle, mythid,
188     & dummy )
189 heimbach 1.1
190 heimbach 1.2 call mdsreadfieldxz( OBSuFile, readBinaryPrec, 'RS',
191     & nr, OBSu, irec, mythid)
192 heimbach 1.1 do bj = jtlo,jthi
193     do bi = itlo,ithi
194 heimbach 1.2 c-- Loop over the model layers
195     fctile = 0. _d 0
196     do k = 1,nr
197     c-- Compute model data misfit and cost function term for
198     c the temperature field.
199     do i = imin,imax
200     j = OB_Js(I,bi,bj)
201     if (maskW(i,j,k,bi,bj) .ne. 0.) then
202     fctile = fctile +
203     & wobcss(k,iobcs)*cosphi(i,j,bi,bj)*
204     & (ubar(i,j,k,bi,bj) - OBSu(i,k,bi,bj))*
205     & (ubar(i,j,k,bi,bj) - OBSu(i,k,bi,bj))
206     endif
207     enddo
208     enddo
209     c-- End of loop over layers.
210     fcthread = fcthread + fctile
211     objf_obcss(bi,bj) = objf_obcss(bi,bj) + fctile
212     enddo
213     enddo
214    
215     c-- vvel
216     iobcs = 4
217     c-- Read time averages and the monthly mean data.
218     il = ilnblnk( vbarfile )
219     write(fnamevvel(1:80),'(2a,i10.10)')
220     & vbarfile(1:il),'.',optimcycle
221     call active_read_xyz( fnamevvel, vbar, irec,
222     & doglobalread, ladinit,
223     & optimcycle, mythid,
224     & dummy )
225 heimbach 1.1
226 heimbach 1.2 call mdsreadfieldxz( OBSvFile, readBinaryPrec, 'RS',
227     & nr, OBSv, irec, mythid)
228    
229     do bj = jtlo,jthi
230     do bi = itlo,ithi
231     c-- Loop over the model layers
232 heimbach 1.1 fctile = 0. _d 0
233 heimbach 1.2 do k = 1,nr
234     c-- Compute model data misfit and cost function term for
235     c the temperature field.
236     do i = imin,imax
237     j = OB_Js(I,bi,bj)
238     if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
239     fctile = fctile +
240     & wobcss(k,iobcs)*cosphi(i,j,bi,bj)*
241     & (vbar(i,j,k,bi,bj) - OBSv(i,k,bi,bj))*
242     & (vbar(i,j,k,bi,bj) - OBSv(i,k,bi,bj))
243     endif
244     enddo
245 heimbach 1.1 enddo
246 heimbach 1.2 c-- End of loop over layers.
247     fcthread = fcthread + fctile
248 heimbach 1.1 objf_obcss(bi,bj) = objf_obcss(bi,bj) + fctile
249     enddo
250     enddo
251    
252     #ifdef ECCO_VERBOSE
253     c-- Print cost function for all tiles.
254     _GLOBAL_SUM_R8( fcthread , myThid )
255     write(msgbuf,'(a)') ' '
256     call print_message( msgbuf, standardmessageunit,
257     & SQUEEZE_RIGHT , mythid)
258     write(msgbuf,'(a,i8.8)')
259     & ' cost_obcss: irec = ',irec
260     call print_message( msgbuf, standardmessageunit,
261     & SQUEEZE_RIGHT , mythid)
262     write(msgbuf,'(a,a,d22.15)')
263     & ' global cost function value',
264     & ' (obcss) = ',fcthread
265     call print_message( msgbuf, standardmessageunit,
266     & SQUEEZE_RIGHT , mythid)
267     write(msgbuf,'(a)') ' '
268     call print_message( msgbuf, standardmessageunit,
269     & SQUEEZE_RIGHT , mythid)
270     #endif
271    
272     enddo
273     c-- End of loop over records.
274    
275     #endif
276    
277     return
278     end
279    
280    
281    
282    
283    
284    
285    

  ViewVC Help
Powered by ViewVC 1.1.22