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

Contents of /MITgcm/pkg/ecco/cost_obcsw.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (show annotations) (download)
Mon Oct 11 16:38:53 2004 UTC (19 years, 7 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: +132 -117 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
2 #include "COST_CPPOPTIONS.h"
3 #ifdef ALLOW_OBCS
4 # include "OBCS_OPTIONS.h"
5 #endif
6
7 subroutine cost_obcsw(
8 I myiter,
9 I mytime,
10 I mythid
11 & )
12
13 c ==================================================================
14 c SUBROUTINE cost_obcsw
15 c ==================================================================
16 c
17 c o cost function contribution obc
18 c
19 c ==================================================================
20 c SUBROUTINE cost_obcsw
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 ip1
60
61 _RL fctile
62 _RL fcthread
63 _RL dummy
64
65 character*(80) fnametheta
66 character*(80) fnamesalt
67 character*(80) fnameuvel
68 character*(80) fnamevvel
69
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_OBCSW_COST_CONTRIBUTION
98
99 ip1 = 1
100 fcthread = 0. _d 0
101
102 c-- Loop over records.
103 do irec = 1,nmonsrec
104
105 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
116 call mdsreadfieldyz( OBWtFile, readBinaryPrec, 'RS',
117 & nr, OBWt, irec, mythid)
118
119 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 j = jmin,jmax
127 i = OB_Iw(J,bi,bj)
128 if (maskS(i+ip1,j,k,bi,bj) .ne. 0.) then
129 fctile = fctile +
130 & wobcsw(k,iobcs)*cosphi(i,j,bi,bj)*
131 & (tbar(i,j,k,bi,bj) - OBWt(j,k,bi,bj))*
132 & (tbar(i,j,k,bi,bj) - OBWt(j,k,bi,bj))
133 endif
134 enddo
135 enddo
136 c-- End of loop over layers.
137 fcthread = fcthread + fctile
138 objf_obcsw(bi,bj) = objf_obcsw(bi,bj) + fctile
139 enddo
140 enddo
141
142 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 mdsreadfieldyz( OBWsFile, readBinaryPrec, 'RS',
154 & nr, OBWs, 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 j = jmin,jmax
164 i = OB_Iw(J,bi,bj)
165 if (maskS(i+ip1,j,k,bi,bj) .ne. 0.) then
166 fctile = fctile +
167 & wobcsw(k,iobcs)*cosphi(i,j,bi,bj)*
168 & (sbar(i,j,k,bi,bj) - OBWs(j,k,bi,bj))*
169 & (sbar(i,j,k,bi,bj) - OBWs(j,k,bi,bj))
170 endif
171 enddo
172 enddo
173 c-- End of loop over layers.
174 fcthread = fcthread + fctile
175 objf_obcsw(bi,bj) = objf_obcsw(bi,bj) + fctile
176 enddo
177 enddo
178
179 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
190 call mdsreadfieldyz( OBWuFile, readBinaryPrec, 'RS',
191 & nr, OBWu, irec, mythid)
192 do bj = jtlo,jthi
193 do bi = itlo,ithi
194 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 j = jmin,jmax
200 i = OB_Iw(J,bi,bj)
201 if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
202 fctile = fctile +
203 & wobcsw(k,iobcs)*cosphi(i,j,bi,bj)*
204 & (ubar(i,j,k,bi,bj) - OBWu(j,k,bi,bj))*
205 & (ubar(i,j,k,bi,bj) - OBWu(j,k,bi,bj))
206 endif
207 enddo
208 enddo
209 c-- End of loop over layers.
210 fcthread = fcthread + fctile
211 objf_obcsw(bi,bj) = objf_obcsw(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
226 call mdsreadfieldyz( OBWvFile, readBinaryPrec, 'RS',
227 & nr, OBWv, irec, mythid)
228
229 do bj = jtlo,jthi
230 do bi = itlo,ithi
231 c-- Loop over the model layers
232 fctile = 0. _d 0
233 do k = 1,nr
234 c-- Compute model data misfit and cost function term for
235 c the temperature field.
236 do j = jmin,jmax
237 i = OB_Iw(J,bi,bj)
238 if (maskS(i,j,k,bi,bj) .ne. 0.) then
239 fctile = fctile +
240 & wobcsw(k,iobcs)*cosphi(i,j,bi,bj)*
241 & (vbar(i,j,k,bi,bj) - OBWv(j,k,bi,bj))*
242 & (vbar(i,j,k,bi,bj) - OBWv(j,k,bi,bj))
243 endif
244 enddo
245 enddo
246 c-- End of loop over layers.
247 fcthread = fcthread + fctile
248 objf_obcsw(bi,bj) = objf_obcsw(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_obcsw: 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 & ' (obcsw) = ',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