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

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

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


Revision 1.6 - (show annotations) (download)
Tue Oct 9 00:02:50 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.5: +18 -17 lines
add missing cvs $Header:$ or $Name:$

1 C $Header: $
2 C $Name: $
3
4 #include "COST_CPPOPTIONS.h"
5
6
7 subroutine cost_CTDT(
8 I myiter,
9 I mytime,
10 I mythid
11 & )
12
13 c ==================================================================
14 c SUBROUTINE cost_CTDT
15 c ==================================================================
16 c
17 c o Evaluate cost function contribution of CTD 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_CTDT
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_ctdt
60 _RL fcthread_ctdt
61 _RL www (1-olx:snx+olx,1-oly:sny+oly)
62 _RL wtmp (1-olx:snx+olx,1-oly:sny+oly)
63 _RL tmpobs (1-olx:snx+olx,1-oly:sny+oly)
64 _RL tmpbar (1-olx:snx+olx,1-oly:sny+oly)
65 _RL cmask (1-olx:snx+olx,1-oly:sny+oly)
66 _RL spval
67 _RL spmax
68
69 character*(80) fnametheta
70
71 logical doglobalread
72 logical ladinit
73
74 character*(MAX_LEN_MBUF) msgbuf
75
76 cnew(
77 integer il
78 integer mody, modm
79 integer iyear, imonth
80 character*(80) fnametmp
81 logical exst
82 cnew)
83
84 c == external functions ==
85
86 integer ilnblnk
87 external ilnblnk
88
89 c == end of interface ==
90
91 jtlo = mybylo(mythid)
92 jthi = mybyhi(mythid)
93 itlo = mybxlo(mythid)
94 ithi = mybxhi(mythid)
95 jmin = 1
96 jmax = sny
97 imin = 1
98 imax = snx
99
100 spval = -1.8
101 spmax = 40.
102
103 c-- Read state record from global file.
104 doglobalread = .false.
105 ladinit = .false.
106
107 #ifdef ALLOW_CTDT_COST_CONTRIBUTION
108
109 if (optimcycle .ge. 0) then
110 ilu=ilnblnk( tbarfile )
111 write(fnametheta(1:80),'(2a,i10.10)')
112 & tbarfile(1:ilu),'.',optimcycle
113 endif
114
115 fcthread_ctdt = 0. _d 0
116
117 cnew(
118 mody = modelstartdate(1)/10000
119 modm = modelstartdate(1)/100 - mody*100
120 cnew)
121
122 c-- Loop over records.
123 do irec = 1,nmonsrec
124
125 c-- Read time averages and the monthly mean data.
126 call active_read_xyz( fnametheta, tbar, irec,
127 & doglobalread, ladinit,
128 & optimcycle, mythid, xx_tbar_mean_dummy )
129
130 cnew(
131 iyear = mody + INT((modm-1+irec-1)/12)
132 imonth = 1 + MOD(modm-1+irec-1,12)
133 il=ilnblnk(ctdtfile)
134 write(fnametmp(1:80),'(2a,i4)')
135 & ctdtfile(1:il), '_', iyear
136 inquire( file=fnametmp, exist=exst )
137 if (.NOT. exst) then
138 write(fnametmp(1:80),'(a)') ctdtfile(1:il)
139 imonth = irec
140 endif
141
142 call mdsreadfield( fnametmp, cost_iprec, 'RL', nr, ctdtobs,
143 & imonth, mythid)
144 cnew)
145
146 c-- Loop over this thread's tiles.
147 do bj = jtlo,jthi
148 do bi = itlo,ithi
149 c-- Loop over the model layers
150
151 fctile_ctdt = 0. _d 0
152
153 do k = 1,nr
154
155 c-- Determine the weights to be used.
156 do j = jmin,jmax
157 do i = imin,imax
158 cmask(i,j) = 1. _d 0
159 if (ctdtobs(i,j,k,bi,bj) .lt. spval .or.
160 & ctdtobs(i,j,k,bi,bj) .gt. spmax .or.
161 & ctdtobs(i,j,k,bi,bj) .eq. 0. ) then
162 cmask(i,j) = 0. _d 0
163 endif
164
165 cph(
166 cph print *, 'WARNING: SPECIFIC SETUP FOR ECCO'
167 cph below statement could be replaced by following
168 cph to make it independnet of Nr:
169 cph
170 cph if ( rC(K) .GT. -1000. ) then
171 cph)
172 c set cmask=0 in areas shallower than 1000m
173
174 if (
175 & (_hFacC(i,j,13,bi,bj) .ne. 0.).and.
176 & (_hFacC(i,j,k,bi,bj) .ne. 0.)) then
177
178 www(i,j) = cosphi(i,j,bi,bj)*cmask(i,j)
179 tmpobs(i,j) = ctdtobs(i,j,k,bi,bj)
180 tmpbar(i,j) = tbar(i,j,k,bi,bj)
181 wtmp(i,j) = wtheta2(i,j,k,bi,bj)
182
183 c-- The array ctdtobs contains CTD temperature.
184 fctile_ctdt = fctile_ctdt +
185 & (wtmp(i,j)*www(i,j))*
186 & (tmpbar(i,j)-tmpobs(i,j))*
187 & (tmpbar(i,j)-tmpobs(i,j))
188 if ( wtmp(i,j)*www(i,j) .ne. 0. )
189 & num_ctdt(bi,bj) = num_ctdt(bi,bj) + 1. _d 0
190 endif
191 enddo
192 enddo
193 enddo
194 c-- End of loop over layers.
195
196 fcthread_ctdt = fcthread_ctdt + fctile_ctdt
197 objf_ctdt(bi,bj) = objf_ctdt(bi,bj) + fctile_ctdt
198
199 #ifdef ECCO_VERBOSE
200 write(msgbuf,'(a)') ' '
201 call print_message( msgbuf, standardmessageunit,
202 & SQUEEZE_RIGHT , mythid)
203 write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
204 & ' COST_CTDT: irec,bi,bj = ',irec,bi,bj
205 call print_message( msgbuf, standardmessageunit,
206 & SQUEEZE_RIGHT , mythid)
207 write(msgbuf,'(a,d22.15)')
208 & ' COST_CTDT: cost function = ', fctile_ctdt
209 call print_message( msgbuf, standardmessageunit,
210 & SQUEEZE_RIGHT , mythid)
211 write(msgbuf,'(a)') ' '
212 call print_message( msgbuf, standardmessageunit,
213 & SQUEEZE_RIGHT , mythid)
214 #endif
215
216 enddo
217 enddo
218
219 #ifdef ECCO_VERBOSE
220 c-- Print cost function for all tiles.
221 _GLOBAL_SUM_R8( fcthread_ctdt , myThid )
222 write(msgbuf,'(a)') ' '
223 call print_message( msgbuf, standardmessageunit,
224 & SQUEEZE_RIGHT , mythid)
225 write(msgbuf,'(a,i8.8)')
226 & ' cost_CTDT: irec = ',irec
227 call print_message( msgbuf, standardmessageunit,
228 & SQUEEZE_RIGHT , mythid)
229 write(msgbuf,'(a,a,d22.15)')
230 & ' global cost function value',
231 & ' ( CTD temp. ) = ',fcthread_ctdt
232 call print_message( msgbuf, standardmessageunit,
233 & SQUEEZE_RIGHT , mythid)
234 write(msgbuf,'(a)') ' '
235 call print_message( msgbuf, standardmessageunit,
236 & SQUEEZE_RIGHT , mythid)
237 #endif
238
239 enddo
240 c-- End of second loop over records.
241
242 #else
243 c-- Do not enter the calculation of the CTD temperature contribution
244 c-- to the final cost function.
245
246 fctile_ctdt = 0. _d 0
247 fcthread_ctdt = 0. _d 0
248
249 crg
250 nrec = 1
251 crg
252
253 _BEGIN_MASTER( mythid )
254 write(msgbuf,'(a)') ' '
255 call print_message( msgbuf, standardmessageunit,
256 & SQUEEZE_RIGHT , mythid)
257 write(msgbuf,'(a,a)')
258 & ' cost_CTDT: no contribution of CTD temperature ',
259 & ' to cost function.'
260 call print_message( msgbuf, standardmessageunit,
261 & SQUEEZE_RIGHT , mythid)
262 write(msgbuf,'(a,a,i9.8)')
263 & ' cost_CTDT: number of records that would have',
264 & ' been processed: ',nrec
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 _END_MASTER( mythid )
271 #endif
272
273 return
274 end

  ViewVC Help
Powered by ViewVC 1.1.22