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

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

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


Revision 1.7 - (show annotations) (download)
Tue Oct 9 00:02:51 2007 UTC (16 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint60, checkpoint61, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59k, checkpoint59j, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61l, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i
Changes since 1.6: +17 -16 lines
add missing cvs $Header:$ or $Name:$

1 C $Header: $
2 C $Name: $
3
4 #include "COST_CPPOPTIONS.h"
5
6
7 subroutine cost_XBT(
8 I myiter,
9 I mytime,
10 I mythid
11 & )
12
13 c ==================================================================
14 c SUBROUTINE cost_XBT
15 c ==================================================================
16 c
17 c o Evaluate cost function contribution of XBT temperature data.
18 c
19 c started: Elisabeth Remy eremy@ucsd.edu 30-Aug-2000
20 c
21 c
22 c ==================================================================
23 c SUBROUTINE cost_XBT
24 c ==================================================================
25
26 implicit none
27
28 c == global variables ==
29
30 #include "EEPARAMS.h"
31 #include "SIZE.h"
32 #include "GRID.h"
33 #include "DYNVARS.h"
34
35 #include "cal.h"
36 #include "ecco_cost.h"
37 #include "ctrl.h"
38 #include "ctrl_dummy.h"
39 #include "optim.h"
40
41 c == routine arguments ==
42
43 integer myiter
44 _RL mytime
45 integer mythid
46
47 c == local variables ==
48
49 integer bi,bj
50 integer i,j,k
51 integer itlo,ithi
52 integer jtlo,jthi
53 integer jmin,jmax
54 integer imin,imax
55 integer nrec
56 integer irec
57 integer ilu
58
59 _RL fctile_xbt
60 _RL fcthread_xbt
61 _RL www (1-olx:snx+olx,1-oly:sny+oly)
62 _RL tmpobs (1-olx:snx+olx,1-oly:sny+oly)
63 _RL cmask (1-olx:snx+olx,1-oly:sny+oly)
64 _RL spval
65 _RL spmax
66 _RL ztop,rl_35,rl_0
67
68 character*(80) fnametheta
69
70 logical doglobalread
71 logical ladinit
72
73 character*(MAX_LEN_MBUF) msgbuf
74
75 cnew(
76 integer il
77 integer mody, modm
78 integer iyear, imonth
79 character*(80) fnametmp
80 logical exst
81 cnew)
82
83 c == external functions ==
84
85 integer ilnblnk
86 external ilnblnk
87 _RL SW_PTMP
88 external SW_PTMP
89
90 c == end of interface ==
91
92 jtlo = mybylo(mythid)
93 jthi = mybyhi(mythid)
94 itlo = mybxlo(mythid)
95 ithi = mybxhi(mythid)
96 jmin = 1
97 jmax = sny
98 imin = 1
99 imax = snx
100
101 spval = -1.8
102 spmax = 40.
103 ztop = -.981*1.027
104 rl_35 = 35.0
105 rl_0 = 0.0
106
107 c-- Read state record from global file.
108 doglobalread = .false.
109 ladinit = .false.
110
111 #ifdef ALLOW_XBT_COST_CONTRIBUTION
112
113 if (optimcycle .ge. 0) then
114 ilu=ilnblnk( tbarfile )
115 write(fnametheta(1:80),'(2a,i10.10)')
116 & tbarfile(1:ilu),'.',optimcycle
117 endif
118
119 fcthread_xbt = 0. _d 0
120
121 cnew(
122 mody = modelstartdate(1)/10000
123 modm = modelstartdate(1)/100 - mody*100
124 cnew)
125
126 c-- Loop over records.
127 do irec = 1,nmonsrec
128
129 c-- Read time averages and the monthly mean data.
130 call active_read_xyz( fnametheta, tbar, irec,
131 & doglobalread, ladinit,
132 & optimcycle, mythid
133 & , xx_tbar_mean_dummy )
134
135 cnew(
136 iyear = mody + INT((modm-1+irec-1)/12)
137 imonth = 1 + MOD(modm-1+irec-1,12)
138 il=ilnblnk(xbtfile)
139 write(fnametmp(1:80),'(2a,i4)')
140 & xbtfile(1:il), '_', iyear
141 inquire( file=fnametmp, exist=exst )
142 if (.NOT. exst) then
143 write(fnametmp(1:80),'(a)') xbtfile(1:il)
144 imonth = irec
145 endif
146
147 call mdsreadfield( fnametmp, cost_iprec, 'RL', nr, xbtobs,
148 & imonth, mythid)
149 cnew)
150
151 c-- Loop over this thread's tiles.
152 do bj = jtlo,jthi
153 do bi = itlo,ithi
154 c-- Loop over the model layers
155
156 fctile_xbt = 0. _d 0
157
158 do k = 1,nr
159
160 c-- Determine the weights to be used.
161 do j = jmin,jmax
162 do i = imin,imax
163 cmask(i,j) = cosphi(i,j,bi,bj)
164 if ( xbtobs(i,j,k,bi,bj) .lt. spval .or.
165 & xbtobs(i,j,k,bi,bj) .gt. spmax .or.
166 & xbtobs(i,j,k,bi,bj) .eq. 0. ) then
167 cmask(i,j) = 0. _d 0
168 endif
169
170 cph(
171 cph print *, 'WARNING: SPECIFIC SETUP FOR ECCO'
172 cph below statement could be replaced by following
173 cph to make it independnet of Nr:
174 cph
175 cph if ( rC(K) .GT. -1000. ) then
176 cph)
177
178 enddo
179 enddo
180
181 do j = jmin,jmax
182 do i = imin,imax
183 if ( _hFacC(i,j,k,bi,bj) .ne. 0. .AND.
184 & _hFacC(i,j,13,bi,bj) .ne. 0. ) then
185 tmpobs(i,j) = SW_PTMP(rl_35,
186 $ xbtobs(i,j,k,bi,bj),ztop*rc(k),rl_0)
187 c-- The array ctdtobs contains CTD temperature.
188
189 fctile_xbt = fctile_xbt +
190 & (wtheta2(i,j,k,bi,bj)*cmask(i,j))*
191 & (tbar(i,j,k,bi,bj)-tmpobs(i,j))*
192 & (tbar(i,j,k,bi,bj)-tmpobs(i,j))
193 if ( wtheta2(i,j,k,bi,bj)*cmask(i,j) .ne. 0. )
194 & num_xbt(bi,bj) = num_xbt(bi,bj) + 1. _d 0
195 endif
196 enddo
197 enddo
198 c-- End of loop over layers.
199 enddo
200
201 fcthread_xbt = fcthread_xbt + fctile_xbt
202 objf_xbt(bi,bj) = objf_xbt(bi,bj) + fctile_xbt
203
204 #ifdef ECCO_VERBOSE
205 write(msgbuf,'(a)') ' '
206 call print_message( msgbuf, standardmessageunit,
207 & SQUEEZE_RIGHT , mythid)
208 write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
209 & ' COST_XBT: irec,bi,bj = ',irec,bi,bj
210 call print_message( msgbuf, standardmessageunit,
211 & SQUEEZE_RIGHT , mythid)
212 write(msgbuf,'(a,d22.15)')
213 & ' COST_XBT: cost function = ', fctile_xbt
214 call print_message( msgbuf, standardmessageunit,
215 & SQUEEZE_RIGHT , mythid)
216 write(msgbuf,'(a)') ' '
217 call print_message( msgbuf, standardmessageunit,
218 & SQUEEZE_RIGHT , mythid)
219 #endif
220
221 enddo
222 enddo
223
224 #ifdef ECCO_VERBOSE
225 c-- Print cost function for all tiles.
226 _GLOBAL_SUM_R8( fcthread_xbt , myThid )
227 write(msgbuf,'(a)') ' '
228 call print_message( msgbuf, standardmessageunit,
229 & SQUEEZE_RIGHT , mythid)
230 write(msgbuf,'(a,i8.8)')
231 & ' cost_XBT: irec = ',irec
232 call print_message( msgbuf, standardmessageunit,
233 & SQUEEZE_RIGHT , mythid)
234 write(msgbuf,'(a,a,d22.15)')
235 & ' global cost function value',
236 & ' ( XBT temp. ) = ',fcthread_xbt
237 call print_message( msgbuf, standardmessageunit,
238 & SQUEEZE_RIGHT , mythid)
239 write(msgbuf,'(a)') ' '
240 call print_message( msgbuf, standardmessageunit,
241 & SQUEEZE_RIGHT , mythid)
242 #endif
243
244 enddo
245 c-- End of second loop over records.
246
247 #else
248 c-- Do not enter the calculation of the XBT temperature contribution
249 c-- to the final cost function.
250
251 fctile_xbt = 0. _d 0
252 fcthread_xbt = 0. _d 0
253
254 crg
255 nrec = 1
256 crg
257
258 _BEGIN_MASTER( mythid )
259 write(msgbuf,'(a)') ' '
260 call print_message( msgbuf, standardmessageunit,
261 & SQUEEZE_RIGHT , mythid)
262 write(msgbuf,'(a,a)')
263 & ' cost_XBT: no contribution of XBT temperature ',
264 & ' to cost function.'
265 call print_message( msgbuf, standardmessageunit,
266 & SQUEEZE_RIGHT , mythid)
267 write(msgbuf,'(a,a,i9.8)')
268 & ' cost_XBT: number of records that would have',
269 & ' been processed: ',nrec
270 call print_message( msgbuf, standardmessageunit,
271 & SQUEEZE_RIGHT , mythid)
272 write(msgbuf,'(a)') ' '
273 call print_message( msgbuf, standardmessageunit,
274 & SQUEEZE_RIGHT , mythid)
275 _END_MASTER( mythid )
276 #endif
277
278 return
279 end

  ViewVC Help
Powered by ViewVC 1.1.22