1 |
C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_gencost_sshv4.F,v 1.35 2016/09/20 15:27:54 gforget Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "ECCO_OPTIONS.h" |
5 |
|
6 |
|
7 |
subroutine cost_gencost_sshv4( |
8 |
I mythid |
9 |
& ) |
10 |
|
11 |
c ================================================================== |
12 |
c SUBROUTINE cost_gencost_sshv4 |
13 |
c ================================================================== |
14 |
c |
15 |
c o Evaluate cost function contributions of sea surface height. |
16 |
c |
17 |
c started: Gael Forget, Oct-2009 |
18 |
c |
19 |
c working assumption for the time mean dynamic topography (MDT) constraint: |
20 |
c the various SLA data sets (tp, ers, gfo) have been consistenty cross-referenced, |
21 |
c as done in the RADS data sets. We do not need to know the reference dynamic |
22 |
c topography (or SSH/Geoid). Simply it is assumed that the biases |
23 |
c between instruments have been taken care of. This is only a matter |
24 |
c for the MDT constraint, not for SLA constraints (see below). |
25 |
c |
26 |
cgf 1) there are a few hardcoded numbers that will eventually be put in common |
27 |
cgf blocks/namelists |
28 |
cgf 2) there are a several refinements that should be considered, such as |
29 |
cgf modulating weights with respect to numbers of samples |
30 |
c |
31 |
c ================================================================== |
32 |
c SUBROUTINE cost_gencost_sshv4 |
33 |
c ================================================================== |
34 |
|
35 |
implicit none |
36 |
|
37 |
c == global variables == |
38 |
|
39 |
#include "EEPARAMS.h" |
40 |
#include "SIZE.h" |
41 |
#include "PARAMS.h" |
42 |
#include "GRID.h" |
43 |
#ifdef ALLOW_CAL |
44 |
# include "cal.h" |
45 |
#endif |
46 |
#ifdef ALLOW_ECCO |
47 |
# include "ecco.h" |
48 |
#endif |
49 |
#ifdef ALLOW_SMOOTH |
50 |
# include "SMOOTH.h" |
51 |
#endif |
52 |
|
53 |
c == routine arguments == |
54 |
|
55 |
integer mythid |
56 |
|
57 |
#ifdef ALLOW_GENCOST_CONTRIBUTION |
58 |
|
59 |
c == functions == |
60 |
LOGICAL MASTER_CPU_THREAD |
61 |
EXTERNAL MASTER_CPU_THREAD |
62 |
|
63 |
c == local variables == |
64 |
|
65 |
integer bi,bj |
66 |
integer i,j,k |
67 |
integer itlo,ithi |
68 |
integer jtlo,jthi |
69 |
integer jmin,jmax |
70 |
integer imin,imax |
71 |
integer irec,jrec,krec |
72 |
integer ilps |
73 |
|
74 |
logical doglobalread |
75 |
logical ladinit |
76 |
|
77 |
c mapping to gencost |
78 |
integer igen_mdt, igen_lsc, igen_gmsl |
79 |
integer igen_tp, igen_ers, igen_gfo |
80 |
integer igen_etaday |
81 |
|
82 |
_RL spval, factor |
83 |
_RL offset |
84 |
_RL offset_sum |
85 |
_RL slaoffset |
86 |
_RL slaoffset_sum |
87 |
_RL slaoffset_weight |
88 |
|
89 |
c local arrays |
90 |
_RL num0array ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) |
91 |
_RL num0total |
92 |
integer tempinteger |
93 |
|
94 |
_RL diagnosfld ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) |
95 |
|
96 |
_RL mdtob(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
97 |
_RL tpob(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
98 |
_RL ersob(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
99 |
_RL gfoob(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
100 |
_RL mdtma(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
101 |
_RL tpma(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
102 |
_RL ersma(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
103 |
_RL gfoma(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
104 |
_RL etaday(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
105 |
|
106 |
character*(MAX_LEN_FNAM) mdtfi, tpfi, ersfi, gfofi |
107 |
integer tpsta(4), erssta(4), gfosta(4) |
108 |
integer mdtsta(4), mdtend(4) |
109 |
_RL tpper, ersper, gfoper |
110 |
|
111 |
c for PART 1: re-reference MDT (mdt) to the inferred SLA reference field |
112 |
_RL psmean ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) |
113 |
_RL mean_slaobs_mdt(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
114 |
_RL mean_slaobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
115 |
|
116 |
c for PART 2: compute time mean differences over the model period |
117 |
_RL mean_slaobs_model(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
118 |
_RL mean_psMssh_all(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
119 |
_RL mean_psMssh_all_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
120 |
_RL mean_psMssh_all_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
121 |
|
122 |
_RL mean_psMslaobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
123 |
_RL mean_psMslaobs_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
124 |
|
125 |
_RL mean_psMtpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
126 |
_RL mean_psMtpobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
127 |
|
128 |
_RL mean_psMersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
129 |
_RL mean_psMersobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
130 |
|
131 |
_RL mean_psMgfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
132 |
_RL mean_psMgfoobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
133 |
|
134 |
c for PART 4/5: compute smooth/raw anomalies |
135 |
_RL anom_psMslaobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
136 |
_RL anom_psMslaobs_NUM (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
137 |
_RL anom_psMslaobs_MSK (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
138 |
_RL anom_psMslaobs_SUB (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
139 |
_RL anom_slaobs (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
140 |
_RL anom_slaobs_SUB (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
141 |
|
142 |
_RL anom_psMtpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
143 |
_RL anom_tpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
144 |
|
145 |
_RL anom_psMersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
146 |
_RL anom_ersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
147 |
|
148 |
_RL anom_psMgfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
149 |
_RL anom_gfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
150 |
|
151 |
integer mdt_y0,mdt_y1,year,day |
152 |
integer num0 |
153 |
|
154 |
_RL junk,junkweight |
155 |
|
156 |
integer ndaysave |
157 |
_RL ndaysaveRL |
158 |
|
159 |
integer k2, k2_mdt, k2_lsc |
160 |
|
161 |
character*(80) fname |
162 |
character*(80) fname4test |
163 |
character*(MAX_LEN_MBUF) msgbuf |
164 |
|
165 |
LOGICAL doReference |
166 |
LOGICAL useEtaMean |
167 |
|
168 |
c == external functions == |
169 |
|
170 |
integer ilnblnk |
171 |
external ilnblnk |
172 |
|
173 |
c == end of interface == |
174 |
|
175 |
jtlo = mybylo(mythid) |
176 |
jthi = mybyhi(mythid) |
177 |
itlo = mybxlo(mythid) |
178 |
ithi = mybxhi(mythid) |
179 |
jmin = 1 |
180 |
jmax = sny |
181 |
imin = 1 |
182 |
imax = snx |
183 |
|
184 |
c-- detect the relevant gencost indices |
185 |
igen_mdt=0 |
186 |
igen_tp =0 |
187 |
igen_ers=0 |
188 |
igen_gfo=0 |
189 |
igen_lsc=0 |
190 |
igen_gmsl=0 |
191 |
igen_etaday=0 |
192 |
do k=1,NGENCOST |
193 |
if (gencost_name(k).EQ.'sshv4-mdt') igen_mdt=k |
194 |
if (gencost_name(k).EQ.'sshv4-tp') igen_tp=k |
195 |
if (gencost_name(k).EQ.'sshv4-ers') igen_ers=k |
196 |
if (gencost_name(k).EQ.'sshv4-gfo') igen_gfo=k |
197 |
if (gencost_name(k).EQ.'sshv4-lsc') igen_lsc=k |
198 |
if (gencost_name(k).EQ.'sshv4-gmsl') igen_gmsl=k |
199 |
enddo |
200 |
|
201 |
k2_mdt=0 |
202 |
k2_lsc=0 |
203 |
do k2 = 1, NGENPPROC |
204 |
if (gencost_posproc(k2,igen_mdt).EQ.'smooth') k2_mdt=k2 |
205 |
if (gencost_posproc(k2,igen_lsc).EQ.'smooth') k2_lsc=k2 |
206 |
enddo |
207 |
|
208 |
call ecco_zero(gencost_weight(:,:,1,1,igen_mdt),1,zeroRL,myThid) |
209 |
call ecco_zero(gencost_weight(:,:,1,1,igen_lsc),1,zeroRL,myThid) |
210 |
call ecco_zero(gencost_weight(:,:,1,1,igen_tp),1,zeroRL,myThid) |
211 |
call ecco_zero(gencost_weight(:,:,1,1,igen_ers),1,zeroRL,myThid) |
212 |
call ecco_zero(gencost_weight(:,:,1,1,igen_gfo),1,zeroRL,myThid) |
213 |
if ( gencost_errfile(igen_mdt) .NE. ' ' ) |
214 |
& call ecco_readwei(gencost_errfile(igen_mdt), |
215 |
& gencost_weight(:,:,1,1,igen_mdt),1,1,mythid) |
216 |
if ( gencost_errfile(igen_lsc) .NE. ' ' ) |
217 |
& call ecco_readwei(gencost_errfile(igen_lsc), |
218 |
& gencost_weight(:,:,1,1,igen_lsc),1,1,mythid) |
219 |
if ( gencost_errfile(igen_tp) .NE. ' ' ) |
220 |
& call ecco_readwei(gencost_errfile(igen_tp), |
221 |
& gencost_weight(:,:,1,1,igen_tp),1,1,mythid) |
222 |
if ( gencost_errfile(igen_ers) .NE. ' ' ) |
223 |
& call ecco_readwei(gencost_errfile(igen_ers), |
224 |
& gencost_weight(:,:,1,1,igen_ers),1,1,mythid) |
225 |
if ( gencost_errfile(igen_gfo) .NE. ' ' ) |
226 |
& call ecco_readwei(gencost_errfile(igen_gfo), |
227 |
& gencost_weight(:,:,1,1,igen_gfo),1,1,mythid) |
228 |
|
229 |
c switch for excluding global mean |
230 |
useEtaMean=.TRUE. |
231 |
|
232 |
write(msgbuf,'(a,l)') |
233 |
& ' sshv4: use global mean of eta useEtaMean=',useEtaMean |
234 |
call print_message( msgbuf, standardmessageunit, |
235 |
& SQUEEZE_RIGHT , mythid) |
236 |
|
237 |
c minimum number of observations to consider re-referenced MDT |
238 |
num0=100 |
239 |
c determine whether or not to re-reference |
240 |
c-- not only model period has to fall within the mdt period |
241 |
c-- but that there has to be at least 365 days of model run |
242 |
c-- so for short model run, this will always get set to FALSE |
243 |
doReference=.FALSE. |
244 |
if ((modelstartdate(1).GT.1992*10000).AND. |
245 |
& (modelstartdate(1).LT.2011*10000).AND. |
246 |
& (ndaysrec.GE.365)) doReference=.TRUE. |
247 |
|
248 |
write(msgbuf,'(a,l)') ' sshv4:re-reference MDT=',doReference |
249 |
call print_message( msgbuf, standardmessageunit, |
250 |
& SQUEEZE_RIGHT , mythid) |
251 |
|
252 |
c-- initialize local arrays |
253 |
do bj = jtlo,jthi |
254 |
do bi = itlo,ithi |
255 |
do j = jmin,jmax |
256 |
do i = imin,imax |
257 |
mdtma(i,j,bi,bj) = 0. _d 0 |
258 |
mdtob(i,j,bi,bj) = 0. _d 0 |
259 |
tpma(i,j,bi,bj) = 0. _d 0 |
260 |
tpob(i,j,bi,bj) = 0. _d 0 |
261 |
ersma(i,j,bi,bj) = 0. _d 0 |
262 |
ersob(i,j,bi,bj) = 0. _d 0 |
263 |
gfoma(i,j,bi,bj) = 0. _d 0 |
264 |
gfoob(i,j,bi,bj) = 0. _d 0 |
265 |
enddo |
266 |
enddo |
267 |
enddo |
268 |
enddo |
269 |
|
270 |
c mdtfi, mdtsta, mdtend |
271 |
if (igen_mdt.GT.0) then |
272 |
ilps=ilnblnk( gencost_datafile(igen_mdt) ) |
273 |
mdtfi=gencost_datafile(igen_mdt)(1:ilps) |
274 |
call cal_CopyDate(gencost_startdate(1,igen_mdt),mdtsta,mythid) |
275 |
call cal_CopyDate(gencost_enddate(1,igen_mdt), mdtend, mythid) |
276 |
endif |
277 |
c tpfi, tpsta, tpper |
278 |
if (igen_tp.GT.0) then |
279 |
ilps=ilnblnk( gencost_datafile(igen_tp) ) |
280 |
tpfi=gencost_datafile(igen_tp)(1:ilps) |
281 |
call cal_CopyDate(gencost_startdate(1,igen_tp),tpsta,mythid) |
282 |
tpper=gencost_period(igen_tp) |
283 |
if (gencost_barfile(igen_tp)(1:5).EQ.'m_eta') |
284 |
& igen_etaday=igen_tp |
285 |
endif |
286 |
c ersfi, erssta, ersper |
287 |
if (igen_ers.GT.0) then |
288 |
ilps=ilnblnk( gencost_datafile(igen_ers) ) |
289 |
ersfi=gencost_datafile(igen_ers)(1:ilps) |
290 |
call cal_CopyDate(gencost_startdate(1,igen_ers),erssta,mythid) |
291 |
ersper=gencost_period(igen_ers) |
292 |
if (gencost_barfile(igen_ers)(1:5).EQ.'m_eta') |
293 |
& igen_etaday=igen_ers |
294 |
endif |
295 |
c gfofi, gfosta, gfoper |
296 |
if (igen_gfo.GT.0) then |
297 |
ilps=ilnblnk( gencost_datafile(igen_gfo) ) |
298 |
gfofi=gencost_datafile(igen_gfo)(1:ilps) |
299 |
call cal_CopyDate(gencost_startdate(1,igen_gfo),gfosta,mythid) |
300 |
gfoper=gencost_period(igen_gfo) |
301 |
if (gencost_barfile(igen_gfo)(1:5).EQ.'m_eta') |
302 |
& igen_etaday=igen_gfo |
303 |
endif |
304 |
|
305 |
c mdt_y0, mdt_y1 |
306 |
mdt_y0=mdtsta(1)/10000 |
307 |
mdt_y1=mdtend(1)/10000 |
308 |
|
309 |
write(msgbuf,'(a,i8,a,i8)') |
310 |
& 'mdt[start,end]date(1): ', mdtsta(1), ',', mdtend(1) |
311 |
call print_message( msgbuf, standardmessageunit, |
312 |
& SQUEEZE_RIGHT , mythid) |
313 |
write(msgbuf,'(a,i4,a,i4)') |
314 |
& ' TP MDT yrrange: ', mdt_y0, ',', mdt_y1 |
315 |
call print_message( msgbuf, standardmessageunit, |
316 |
& SQUEEZE_RIGHT , mythid) |
317 |
|
318 |
c-- First, read tiled data. |
319 |
doglobalread = .false. |
320 |
ladinit = .false. |
321 |
|
322 |
write(fname(1:80),'(80a)') ' ' |
323 |
ilps=ilnblnk( gencost_barfile(igen_etaday) ) |
324 |
write(fname(1:80),'(2a,i10.10)') |
325 |
& gencost_barfile(igen_etaday)(1:ilps),'.',eccoiter |
326 |
|
327 |
cgf ======================================================= |
328 |
cgf PART 1: |
329 |
cgf x Get the MDT (mdt) ... compute the sample mean |
330 |
cgf (mean_slaobs_mdt) of the SLA data (i.e. RADS for tp, ers, and gfo |
331 |
cgf together) over the time interval of the MDT ... subtract |
332 |
cgf mean_slaobs_mdt from mdt. |
333 |
cgf x At this point, mdt is the inferred SLA reference field. |
334 |
cgf x From there on, mdt+sla will be directly comparable to |
335 |
cgf the model SSH (etaday). |
336 |
cgf ======================================================= |
337 |
|
338 |
c-- Read mean field and mask |
339 |
|
340 |
if (using_mdt) |
341 |
&call mdsreadfield( mdtfi, cost_iprec, cost_yftype, 1, |
342 |
& mdtob, 1, mythid ) |
343 |
|
344 |
factor = 0.01 _d 0 |
345 |
spval = -9990. _d 0 |
346 |
|
347 |
do bj = jtlo,jthi |
348 |
do bi = itlo,ithi |
349 |
do j = jmin,jmax |
350 |
do i = imin,imax |
351 |
if ( (maskC(i,j,1,bi,bj) .eq. 0.).OR. |
352 |
#ifndef ALLOW_SHALLOW_ALTIMETRY |
353 |
& (R_low(i,j,bi,bj).GT.-200.).OR. |
354 |
#endif |
355 |
& (mdtob(i,j,bi,bj) .lt. spval ).OR. |
356 |
& (mdtob(i,j,bi,bj) .eq. 0. _d 0) ) then |
357 |
mdtma(i,j,bi,bj) = 0. _d 0 |
358 |
mdtob(i,j,bi,bj) = 0. _d 0 |
359 |
else |
360 |
mdtma(i,j,bi,bj) = 1. _d 0 |
361 |
mdtob(i,j,bi,bj) = mdtob(i,j,bi,bj)*factor |
362 |
endif |
363 |
enddo |
364 |
enddo |
365 |
enddo |
366 |
enddo |
367 |
|
368 |
c-- Compute mean_slaobs_mdt: sample mean SLA over the time period of mdt. |
369 |
|
370 |
do bj = jtlo,jthi |
371 |
do bi = itlo,ithi |
372 |
do j = jmin,jmax |
373 |
do i = imin,imax |
374 |
mean_slaobs_mdt(i,j,bi,bj) = 0. _d 0 |
375 |
mean_slaobs_NUM(i,j,bi,bj) = 0. _d 0 |
376 |
enddo |
377 |
enddo |
378 |
enddo |
379 |
enddo |
380 |
|
381 |
do year=mdt_y0,mdt_y1 |
382 |
do day=1,366 |
383 |
|
384 |
do bj = jtlo,jthi |
385 |
do bi = itlo,ithi |
386 |
do j = jmin,jmax |
387 |
do i = imin,imax |
388 |
tpma(i,j,bi,bj) = 0. _d 0 |
389 |
tpob(i,j,bi,bj) = 0. _d 0 |
390 |
ersma(i,j,bi,bj) = 0. _d 0 |
391 |
ersob(i,j,bi,bj) = 0. _d 0 |
392 |
gfoma(i,j,bi,bj) = 0. _d 0 |
393 |
gfoob(i,j,bi,bj) = 0. _d 0 |
394 |
enddo |
395 |
enddo |
396 |
enddo |
397 |
enddo |
398 |
|
399 |
if (using_tpj) |
400 |
& call cost_sla_read_yd( tpfi, tpsta, |
401 |
& tpob, tpma, |
402 |
& year, day, mythid ) |
403 |
if (using_ers) |
404 |
& call cost_sla_read_yd( ersfi, erssta, |
405 |
& ersob, ersma, |
406 |
& year, day, mythid ) |
407 |
if (using_gfo) |
408 |
& call cost_sla_read_yd( gfofi, gfosta, |
409 |
& gfoob, gfoma, |
410 |
& year, day, mythid ) |
411 |
|
412 |
do bj = jtlo,jthi |
413 |
do bi = itlo,ithi |
414 |
do j = jmin,jmax |
415 |
do i = imin,imax |
416 |
if ( tpma(i,j,bi,bj)*mdtma(i,j,bi,bj)* |
417 |
& gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then |
418 |
mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+ |
419 |
& tpob(i,j,bi,bj) |
420 |
mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0 |
421 |
endif |
422 |
if ( ersma(i,j,bi,bj)*mdtma(i,j,bi,bj)* |
423 |
& gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then |
424 |
mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+ |
425 |
& ersob(i,j,bi,bj) |
426 |
mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0 |
427 |
endif |
428 |
if ( gfoma(i,j,bi,bj)*mdtma(i,j,bi,bj)* |
429 |
& gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then |
430 |
mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+ |
431 |
& gfoob(i,j,bi,bj) |
432 |
mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0 |
433 |
endif |
434 |
enddo |
435 |
enddo |
436 |
enddo |
437 |
enddo |
438 |
|
439 |
enddo !do day=1,366 |
440 |
enddo !do year=mdt_y0,mdt_y1 |
441 |
|
442 |
do bj = jtlo,jthi |
443 |
do bi = itlo,ithi |
444 |
do j = jmin,jmax |
445 |
do i = imin,imax |
446 |
if ( ( mean_slaobs_NUM(i,j,bi,bj) .NE. 0. ).AND. |
447 |
& ( maskc(i,j,1,bi,bj) .NE. 0. ).AND. |
448 |
& ( mdtma(i,j,bi,bj) .NE. 0. ) ) then |
449 |
mean_slaobs_mdt(i,j,bi,bj) = |
450 |
& mean_slaobs_mdt(i,j,bi,bj) / |
451 |
& mean_slaobs_NUM(i,j,bi,bj) |
452 |
else |
453 |
mean_slaobs_mdt(i,j,bi,bj) = 0. _d 0 |
454 |
endif |
455 |
enddo |
456 |
enddo |
457 |
enddo |
458 |
enddo |
459 |
|
460 |
|
461 |
c-- smooth mean_slaobs_mdt: |
462 |
|
463 |
if (gencost_outputlevel(igen_mdt).GT.0) then |
464 |
write(fname4test(1:80),'(1a)') 'sla2mdt_raw' |
465 |
call mdswritefield(fname4test,32,.false.,'RL', |
466 |
& 1,mean_slaobs_mdt,1,1,mythid) |
467 |
endif |
468 |
|
469 |
#ifdef ALLOW_SMOOTH |
470 |
if ( useSMOOTH.AND.(k2_mdt.GT.0) ) |
471 |
& call smooth_hetero2d(mean_slaobs_mdt,maskc, |
472 |
& gencost_posproc_c(k2_mdt,igen_mdt), |
473 |
& gencost_posproc_i(k2_mdt,igen_mdt),mythid) |
474 |
#endif |
475 |
|
476 |
if (gencost_outputlevel(igen_mdt).GT.0) then |
477 |
write(fname4test(1:80),'(1a)') 'sla2mdt_smooth' |
478 |
call mdswritefield(fname4test,32,.false.,'RL', |
479 |
& 1,mean_slaobs_mdt,1,1,mythid) |
480 |
endif |
481 |
|
482 |
c-- re-reference mdt: |
483 |
do bj = jtlo,jthi |
484 |
do bi = itlo,ithi |
485 |
do j = jmin,jmax |
486 |
do i = imin,imax |
487 |
if ( ( mdtma(i,j,bi,bj) .NE. 0. ).AND. |
488 |
& ( maskc(i,j,1,bi,bj) .NE. 0. ).AND. |
489 |
& ( doReference ) ) then |
490 |
mdtob(i,j,bi,bj) = mdtob(i,j,bi,bj) |
491 |
& -mean_slaobs_mdt(i,j,bi,bj) |
492 |
endif |
493 |
enddo |
494 |
enddo |
495 |
enddo |
496 |
enddo |
497 |
|
498 |
|
499 |
cgf ======================================================= |
500 |
cgf PART 2: compute sample means of etaday-slaobs over the |
501 |
cgf period that is covered by the model (i.e. etaday). |
502 |
cgf x for all SLA data sets together: mean_psMssh_all, mean_psMssh_all_MSK, |
503 |
cgf and offset will be used in PART 3 (MDT cost term). |
504 |
cgf x for each SLA data individually. mean_psMtpobs, mean_psMtpobs_MS, etc. |
505 |
cgf will be used in PART 4&5 (SLA cost terms). |
506 |
cgf ======================================================= |
507 |
|
508 |
do bj = jtlo,jthi |
509 |
do bi = itlo,ithi |
510 |
do j = jmin,jmax |
511 |
do i = imin,imax |
512 |
|
513 |
psmean(i,j,bi,bj) = 0. _d 0 |
514 |
mean_psMslaobs(i,j,bi,bj) = 0. _d 0 |
515 |
mean_psMtpobs(i,j,bi,bj) = 0. _d 0 |
516 |
mean_psMersobs(i,j,bi,bj) = 0. _d 0 |
517 |
mean_psMgfoobs(i,j,bi,bj) = 0. _d 0 |
518 |
mean_psMssh_all(i,j,bi,bj) = 0. _d 0 |
519 |
mean_slaobs_model(i,j,bi,bj) = 0. _d 0 |
520 |
|
521 |
mean_psMtpobs_NUM(i,j,bi,bj) = 0. _d 0 |
522 |
mean_psMersobs_NUM(i,j,bi,bj) = 0. _d 0 |
523 |
mean_psMgfoobs_NUM(i,j,bi,bj) = 0. _d 0 |
524 |
mean_psMssh_all_NUM(i,j,bi,bj) = 0. _d 0 |
525 |
mean_psMslaobs_MSK(i,j,bi,bj) = 0. _d 0 |
526 |
|
527 |
enddo |
528 |
enddo |
529 |
enddo |
530 |
enddo |
531 |
offset = 0. _d 0 |
532 |
offset_sum = 0. _d 0 |
533 |
|
534 |
|
535 |
do irec = 1, ndaysrec |
536 |
|
537 |
#ifdef ALLOW_AUTODIFF |
538 |
call active_read_xy( fname, etaday, irec, doglobalread, |
539 |
& ladinit, eccoiter, mythid, |
540 |
& gencost_dummy(igen_etaday) ) |
541 |
#else |
542 |
CALL READ_REC_XY_RL( fname, etaday, iRec, 1, myThid ) |
543 |
#endif |
544 |
|
545 |
|
546 |
if (.NOT.useEtaMean) |
547 |
& CALL REMOVE_MEAN_RL( 1, etaday, maskInC, maskInC, rA, drF, |
548 |
& 'etaday', 0. _d 0, myThid ) |
549 |
|
550 |
do bj = jtlo,jthi |
551 |
do bi = itlo,ithi |
552 |
do j = jmin,jmax |
553 |
do i = imin,imax |
554 |
tpma(i,j,bi,bj) = 0. _d 0 |
555 |
tpob(i,j,bi,bj) = 0. _d 0 |
556 |
ersma(i,j,bi,bj) = 0. _d 0 |
557 |
ersob(i,j,bi,bj) = 0. _d 0 |
558 |
gfoma(i,j,bi,bj) = 0. _d 0 |
559 |
gfoob(i,j,bi,bj) = 0. _d 0 |
560 |
enddo |
561 |
enddo |
562 |
enddo |
563 |
enddo |
564 |
|
565 |
if (using_tpj) |
566 |
& call cost_sla_read( tpfi, tpsta, tpper, |
567 |
& zeroRL, zeroRL, |
568 |
& tpob, tpma, |
569 |
& irec, mythid ) |
570 |
if (using_ers) |
571 |
& call cost_sla_read( ersfi, erssta, ersper, |
572 |
& zeroRL, zeroRL, |
573 |
& ersob, ersma, |
574 |
& irec, mythid ) |
575 |
if (using_gfo) |
576 |
& call cost_sla_read( gfofi, gfosta, gfoper, |
577 |
& zeroRL, zeroRL, |
578 |
& gfoob, gfoma, |
579 |
& irec, mythid ) |
580 |
|
581 |
do bj = jtlo,jthi |
582 |
do bi = itlo,ithi |
583 |
do j = jmin,jmax |
584 |
do i = imin,imax |
585 |
psmean(i,j,bi,bj) = psmean(i,j,bi,bj) + |
586 |
& etaday(i,j,bi,bj) / float(ndaysrec) |
587 |
if ( tpma(i,j,bi,bj)* |
588 |
& gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then |
589 |
mean_slaobs_model(i,j,bi,bj)= |
590 |
& mean_slaobs_model(i,j,bi,bj)+tpob(i,j,bi,bj) |
591 |
mean_psMtpobs(i,j,bi,bj) = |
592 |
& mean_psMtpobs(i,j,bi,bj) + |
593 |
& etaday(i,j,bi,bj)-tpob(i,j,bi,bj) |
594 |
mean_psMtpobs_NUM(i,j,bi,bj) = |
595 |
& mean_psMtpobs_NUM(i,j,bi,bj) + 1. _d 0 |
596 |
endif |
597 |
if ( ersma(i,j,bi,bj)* |
598 |
& gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then |
599 |
mean_slaobs_model(i,j,bi,bj)= |
600 |
& mean_slaobs_model(i,j,bi,bj)+ersob(i,j,bi,bj) |
601 |
mean_psMersobs(i,j,bi,bj) = |
602 |
& mean_psMersobs(i,j,bi,bj) + |
603 |
& etaday(i,j,bi,bj)-ersob(i,j,bi,bj) |
604 |
mean_psMersobs_NUM(i,j,bi,bj) = |
605 |
& mean_psMersobs_NUM(i,j,bi,bj) + 1. _d 0 |
606 |
endif |
607 |
if ( gfoma(i,j,bi,bj)* |
608 |
& gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then |
609 |
mean_slaobs_model(i,j,bi,bj)= |
610 |
& mean_slaobs_model(i,j,bi,bj)+gfoob(i,j,bi,bj) |
611 |
mean_psMgfoobs(i,j,bi,bj) = |
612 |
& mean_psMgfoobs(i,j,bi,bj) + |
613 |
& etaday(i,j,bi,bj)-gfoob(i,j,bi,bj) |
614 |
mean_psMgfoobs_NUM(i,j,bi,bj) = |
615 |
& mean_psMgfoobs_NUM(i,j,bi,bj) + 1. _d 0 |
616 |
endif |
617 |
enddo |
618 |
enddo |
619 |
enddo |
620 |
enddo |
621 |
|
622 |
c-- END loop over records for the first time. |
623 |
enddo |
624 |
|
625 |
call ecco_zero(num0array,1,oneRL,mythid) |
626 |
|
627 |
do bj = jtlo,jthi |
628 |
do bi = itlo,ithi |
629 |
do j = jmin,jmax |
630 |
do i = imin,imax |
631 |
if ( ( mean_psMtpobs_NUM(i,j,bi,bj) .NE. 0. ) |
632 |
& ) then |
633 |
mean_psMssh_all(i,j,bi,bj) = |
634 |
& mean_psMssh_all(i,j,bi,bj) + |
635 |
& mean_psMtpobs(i,j,bi,bj) |
636 |
mean_psMssh_all_NUM(i,j,bi,bj) = |
637 |
& mean_psMssh_all_NUM(i,j,bi,bj) + |
638 |
& mean_psMtpobs_NUM(i,j,bi,bj) |
639 |
mean_psMtpobs(i,j,bi,bj) = |
640 |
& mean_psMtpobs(i,j,bi,bj) / |
641 |
& mean_psMtpobs_NUM(i,j,bi,bj) |
642 |
endif |
643 |
if ( ( mean_psMersobs_NUM(i,j,bi,bj) .NE. 0. ) |
644 |
& ) then |
645 |
mean_psMssh_all(i,j,bi,bj) = |
646 |
& mean_psMssh_all(i,j,bi,bj) + |
647 |
& mean_psMersobs(i,j,bi,bj) |
648 |
mean_psMssh_all_NUM(i,j,bi,bj) = |
649 |
& mean_psMssh_all_NUM(i,j,bi,bj) + |
650 |
& mean_psMersobs_NUM(i,j,bi,bj) |
651 |
mean_psMersobs(i,j,bi,bj) = |
652 |
& mean_psMersobs(i,j,bi,bj) / |
653 |
& mean_psMersobs_NUM(i,j,bi,bj) |
654 |
endif |
655 |
if ( ( mean_psMgfoobs_NUM(i,j,bi,bj) .NE. 0. ) |
656 |
& ) then |
657 |
mean_psMssh_all(i,j,bi,bj) = |
658 |
& mean_psMssh_all(i,j,bi,bj) + |
659 |
& mean_psMgfoobs(i,j,bi,bj) |
660 |
mean_psMssh_all_NUM(i,j,bi,bj) = |
661 |
& mean_psMssh_all_NUM(i,j,bi,bj) + |
662 |
& mean_psMgfoobs_NUM(i,j,bi,bj) |
663 |
mean_psMgfoobs(i,j,bi,bj) = |
664 |
& mean_psMgfoobs(i,j,bi,bj) / |
665 |
& mean_psMgfoobs_NUM(i,j,bi,bj) |
666 |
endif |
667 |
|
668 |
if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .GT. 0 ).AND. |
669 |
& ( maskc(i,j,1,bi,bj) .NE. 0. ) |
670 |
& ) then |
671 |
mean_psMslaobs(i,j,bi,bj) = |
672 |
& mean_psMssh_all(i,j,bi,bj) / |
673 |
& mean_psMssh_all_NUM(i,j,bi,bj) |
674 |
mean_psMslaobs_MSK(i,j,bi,bj) = 1. _d 0 |
675 |
else |
676 |
mean_psMslaobs(i,j,bi,bj) = 0. _d 0 |
677 |
mean_psMslaobs_MSK(i,j,bi,bj) = 0. _d 0 |
678 |
endif |
679 |
|
680 |
if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .GT. num0 ).AND. |
681 |
& ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND. |
682 |
& ( mdtma(i,j,bi,bj) .NE. 0. ).AND. |
683 |
& ( doReference ) ) then |
684 |
mean_slaobs_model(i,j,bi,bj) = |
685 |
& mean_slaobs_model(i,j,bi,bj) / |
686 |
& mean_psMssh_all_NUM(i,j,bi,bj) |
687 |
mean_psMssh_all(i,j,bi,bj) = |
688 |
& mean_psMssh_all(i,j,bi,bj) / |
689 |
& mean_psMssh_all_NUM(i,j,bi,bj)-mdtob(i,j,bi,bj) |
690 |
mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0 |
691 |
offset=offset+RA(i,j,bi,bj)*mean_psMssh_all(i,j,bi,bj) |
692 |
offset_sum=offset_sum+RA(i,j,bi,bj) |
693 |
elseif ( ( mdtma(i,j,bi,bj) .NE. 0. ) .AND. |
694 |
& ( maskc(i,j,1,bi,bj) .NE. 0. ).AND. |
695 |
& ( .NOT.doReference ) ) then |
696 |
mean_slaobs_model(i,j,bi,bj) = 0.d0 |
697 |
mean_psMssh_all(i,j,bi,bj) = 0. _d 0 |
698 |
mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0 |
699 |
offset=offset+RA(i,j,bi,bj) |
700 |
& *( psmean(i,j,bi,bj) -mdtob(i,j,bi,bj) ) |
701 |
offset_sum=offset_sum+RA(i,j,bi,bj) |
702 |
num0array(i,j,bi,bj)=0. _d 0 |
703 |
else |
704 |
mean_slaobs_model(i,j,bi,bj) = 0.d0 |
705 |
mean_psMssh_all(i,j,bi,bj) = 0. _d 0 |
706 |
mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0 |
707 |
num0array(i,j,bi,bj)=0. _d 0 |
708 |
endif |
709 |
enddo |
710 |
enddo |
711 |
enddo |
712 |
enddo |
713 |
|
714 |
|
715 |
c-- Do a global summation. |
716 |
_GLOBAL_SUM_RL( offset , mythid ) |
717 |
_GLOBAL_SUM_RL( offset_sum , mythid ) |
718 |
|
719 |
catn--- add warning that written mean_slaobs_model is all zero |
720 |
c due to having less data points that hardcoded num0 |
721 |
num0total=0. _d 0 |
722 |
do bj = jtlo,jthi |
723 |
do bi = itlo,ithi |
724 |
do j = jmin,jmax |
725 |
do i = imin,imax |
726 |
num0total=num0total+num0array(i,j,bi,bj) |
727 |
enddo |
728 |
enddo |
729 |
enddo |
730 |
enddo |
731 |
if(num0total.lt.1. _d 0) then |
732 |
write(msgbuf,'(A,i5,A)') |
733 |
& '** WARNING ** S/R COST_GENCOST_SSHV4: There are <',num0,' pts' |
734 |
call print_message( msgbuf, errormessageunit, |
735 |
& SQUEEZE_RIGHT , mythid) |
736 |
write(msgbuf,'(A)') |
737 |
& ' at all grid points for combined tp+ers+gfo.' |
738 |
call print_message( msgbuf, errormessageunit, |
739 |
& SQUEEZE_RIGHT , mythid) |
740 |
write(msgbuf,'(A)') |
741 |
& 'So, all model slaob minus model sla2model_raw are set to 0.' |
742 |
call print_message( msgbuf, errormessageunit, |
743 |
& SQUEEZE_RIGHT , mythid) |
744 |
endif |
745 |
catn------- |
746 |
|
747 |
if (offset_sum .eq. 0.0) then |
748 |
if (gencost_outputlevel(igen_gmsl).GT.0) then |
749 |
_BEGIN_MASTER( mythid ) |
750 |
write(msgbuf,'(a)') ' sshv4: offset_sum = zero!' |
751 |
call print_message( msgbuf, standardmessageunit, |
752 |
& SQUEEZE_RIGHT , mythid) |
753 |
_END_MASTER( mythid ) |
754 |
endif |
755 |
c stop ' ... stopped in cost_ssh.' |
756 |
else |
757 |
if (gencost_outputlevel(igen_gmsl).GT.0) then |
758 |
_BEGIN_MASTER( mythid ) |
759 |
write(msgbuf,'(a,2d22.15)') ' sshv4:offset,sum=', |
760 |
& offset,offset_sum |
761 |
call print_message( msgbuf, standardmessageunit, |
762 |
& SQUEEZE_RIGHT , mythid) |
763 |
_END_MASTER( mythid ) |
764 |
endif |
765 |
endif |
766 |
|
767 |
c-- Compute (average) offset |
768 |
offset = offset / offset_sum |
769 |
|
770 |
c-- subtract offset from mean_psMssh_all and psmean |
771 |
do bj = jtlo,jthi |
772 |
do bi = itlo,ithi |
773 |
do j = jmin,jmax |
774 |
do i = imin,imax |
775 |
|
776 |
if ( (mean_psMssh_all_NUM(i,j,bi,bj) .GT. num0) .AND. |
777 |
& ( mdtma(i,j,bi,bj) .NE. 0. ) .AND. |
778 |
& ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND. |
779 |
& ( doReference ) ) then |
780 |
c use the re-referencing approach |
781 |
mean_psMssh_all(i,j,bi,bj) = |
782 |
& mean_psMssh_all(i,j,bi,bj) - offset |
783 |
mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0 |
784 |
|
785 |
elseif ( ( mdtma(i,j,bi,bj) .NE. 0. ) .AND. |
786 |
& ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND. |
787 |
& ( .NOT.doReference ) ) then |
788 |
c use the simpler approach |
789 |
mean_psMssh_all(i,j,bi,bj) = |
790 |
& psmean(i,j,bi,bj) -mdtob(i,j,bi,bj) - offset |
791 |
mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0 |
792 |
|
793 |
else |
794 |
mean_psMssh_all(i,j,bi,bj) = 0. _d 0 |
795 |
mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0 |
796 |
endif |
797 |
|
798 |
if ( maskc(i,j,1,bi,bj) .NE. 0. ) |
799 |
& psmean(i,j,bi,bj)=psmean(i,j,bi,bj)-offset |
800 |
|
801 |
enddo |
802 |
enddo |
803 |
enddo |
804 |
enddo |
805 |
|
806 |
c-- smooth mean_psMssh_all |
807 |
if (gencost_outputlevel(igen_mdt).GT.0) then |
808 |
write(fname4test(1:80),'(1a)') 'mdtdiff_raw' |
809 |
call mdswritefield(fname4test,32,.false.,'RL', |
810 |
& 1,mean_psMssh_all,1,1,mythid) |
811 |
endif |
812 |
|
813 |
#ifdef ALLOW_SMOOTH |
814 |
if ( useSMOOTH.AND.(k2_mdt.GT.0) ) |
815 |
& call smooth_hetero2d(mean_psMssh_all,maskc, |
816 |
& gencost_posproc_c(k2_mdt,igen_mdt), |
817 |
& gencost_posproc_i(k2_mdt,igen_mdt),mythid) |
818 |
#endif |
819 |
|
820 |
if (gencost_outputlevel(igen_mdt).GT.0) then |
821 |
write(fname4test(1:80),'(1a)') 'mdtdiff_smooth' |
822 |
call mdswritefield(fname4test,32,.false.,'RL', |
823 |
& 1,mean_psMssh_all,1,1,mythid) |
824 |
endif |
825 |
|
826 |
if (gencost_outputlevel(igen_mdt).GT.0) then |
827 |
write(fname4test(1:80),'(1a)') 'sla2model_raw' |
828 |
call mdswritefield(fname4test,32,.false.,'RL', |
829 |
& 1,mean_slaobs_model,1,1,mythid) |
830 |
endif |
831 |
|
832 |
#ifdef ALLOW_SMOOTH |
833 |
if ( useSMOOTH.AND.(k2_mdt.GT.0) ) |
834 |
& call smooth_hetero2d(mean_slaobs_model,maskc, |
835 |
& gencost_posproc_c(k2_mdt,igen_mdt), |
836 |
& gencost_posproc_i(k2_mdt,igen_mdt),mythid) |
837 |
#endif |
838 |
|
839 |
if (gencost_outputlevel(igen_mdt).GT.0) then |
840 |
write(fname4test(1:80),'(1a)') 'sla2model_smooth' |
841 |
call mdswritefield(fname4test,32,.false.,'RL', |
842 |
& 1,mean_slaobs_model,1,1,mythid) |
843 |
endif |
844 |
|
845 |
cgf at this point: |
846 |
cgf 1) mean_psMssh_all is the sample mean <etaday-slaobs-mdt-offset>, |
847 |
cgf to which a smoothing filter has been applied. |
848 |
cgf 2) mean_psMtpobs is the (unsmoothed) sample mean <etaday-tpobs>. |
849 |
cgf And similarly for ers and gfo, each treated separately. |
850 |
|
851 |
cgf ======================================================= |
852 |
cgf PART 3: compute MDT cost term |
853 |
cgf ======================================================= |
854 |
|
855 |
|
856 |
do bj = jtlo,jthi |
857 |
do bi = itlo,ithi |
858 |
do j = jmin,jmax |
859 |
do i = imin,imax |
860 |
if (mean_psMssh_all_MSK(i,j,bi,bj).NE.0. _d 0) then |
861 |
junk = mean_psMssh_all(i,j,bi,bj) |
862 |
junkweight = gencost_weight(i,j,bi,bj,igen_mdt) |
863 |
& * mdtma(i,j,bi,bj) |
864 |
else |
865 |
junk = 0. _d 0 |
866 |
junkweight = 0. _d 0 |
867 |
endif |
868 |
objf_gencost(bi,bj,igen_mdt) = objf_gencost(bi,bj,igen_mdt) |
869 |
& + junk*junk*junkweight |
870 |
if ( junkweight .ne. 0. ) num_gencost(bi,bj,igen_mdt) = |
871 |
& num_gencost(bi,bj,igen_mdt) + 1. _d 0 |
872 |
diagnosfld(i,j,bi,bj) = junk*junk*junkweight |
873 |
enddo |
874 |
enddo |
875 |
enddo |
876 |
enddo |
877 |
|
878 |
if (gencost_outputlevel(igen_mdt).GT.0) then |
879 |
write(fname4test(1:80),'(1a)') 'misfits_mdt' |
880 |
call mdswritefield(fname4test,32,.false.,'RL', |
881 |
& 1,diagnosfld,1,1,mythid) |
882 |
endif |
883 |
|
884 |
cgf ======================================================= |
885 |
cgf PART 4: compute smooth SLA cost term |
886 |
cgf ======================================================= |
887 |
|
888 |
|
889 |
ndaysave=35 |
890 |
ndaysaveRL=ndaysave |
891 |
|
892 |
catn add a warning if have less nrecday than the hardcoded 35days |
893 |
tempinteger=ndaysrec-ndaysave+1 |
894 |
if(tempinteger.lt.1) then |
895 |
write(msgbuf,'(A,i5)') |
896 |
& '** WARNING ** S/R COST_GENCOST_SSHV4: There are < ',ndaysave |
897 |
call print_message( msgbuf, errormessageunit, |
898 |
& SQUEEZE_RIGHT , mythid) |
899 |
write(msgbuf,'(A)') |
900 |
& ' days required to calculate running mean tp+ers+gfo.' |
901 |
call print_message( msgbuf, errormessageunit, |
902 |
& SQUEEZE_RIGHT , mythid) |
903 |
write(msgbuf,'(A)') |
904 |
& 'PART 4 in cost_gencost_sshv4 is skipped entirely, and there' |
905 |
call print_message( msgbuf, errormessageunit, |
906 |
& SQUEEZE_RIGHT , mythid) |
907 |
write(msgbuf,'(A)') |
908 |
& 'will be NO report of [tp,ers,gfo,lsc,gmsl] in costfunction.' |
909 |
call print_message( msgbuf, errormessageunit, |
910 |
& SQUEEZE_RIGHT , mythid) |
911 |
endif |
912 |
catn ----------- |
913 |
|
914 |
do irec = 1, ndaysrec-ndaysave+1 |
915 |
|
916 |
do bj = jtlo,jthi |
917 |
do bi = itlo,ithi |
918 |
do j = jmin,jmax |
919 |
do i = imin,imax |
920 |
anom_psMslaobs(i,j,bi,bj) = 0. _d 0 |
921 |
anom_psMslaobs_MSK(i,j,bi,bj) = 0. _d 0 |
922 |
anom_psMslaobs_SUB(i,j,bi,bj) = 0. _d 0 |
923 |
anom_slaobs(i,j,bi,bj) = 0. _d 0 |
924 |
anom_slaobs_SUB(i,j,bi,bj) = 0. _d 0 |
925 |
anom_psMslaobs_NUM(i,j,bi,bj) = 0. _d 0 |
926 |
enddo |
927 |
enddo |
928 |
enddo |
929 |
enddo |
930 |
|
931 |
c PART 4.1: compute running sample average over ndaysave |
932 |
c ------------------------------------------------------ |
933 |
|
934 |
do jrec=1,ndaysave |
935 |
|
936 |
krec=irec+jrec-1 |
937 |
|
938 |
#ifdef ALLOW_AUTODIFF |
939 |
call active_read_xy( fname, etaday, krec, doglobalread, |
940 |
& ladinit, eccoiter, mythid, |
941 |
& gencost_dummy(igen_etaday) ) |
942 |
#else |
943 |
CALL READ_REC_XY_RL( fname, etaday, kRec, 1, myThid ) |
944 |
#endif |
945 |
|
946 |
if (.NOT.useEtaMean) |
947 |
& CALL REMOVE_MEAN_RL( 1, etaday, maskInC, maskInC, rA, drF, |
948 |
& 'etaday', 0. _d 0, myThid ) |
949 |
|
950 |
do bj = jtlo,jthi |
951 |
do bi = itlo,ithi |
952 |
do j = jmin,jmax |
953 |
do i = imin,imax |
954 |
tpma(i,j,bi,bj) = 0. _d 0 |
955 |
tpob(i,j,bi,bj) = 0. _d 0 |
956 |
ersma(i,j,bi,bj) = 0. _d 0 |
957 |
ersob(i,j,bi,bj) = 0. _d 0 |
958 |
gfoma(i,j,bi,bj) = 0. _d 0 |
959 |
gfoob(i,j,bi,bj) = 0. _d 0 |
960 |
enddo |
961 |
enddo |
962 |
enddo |
963 |
enddo |
964 |
|
965 |
if (using_tpj) |
966 |
& call cost_sla_read( tpfi, tpsta, tpper, |
967 |
& zeroRL, zeroRL, |
968 |
& tpob, tpma, |
969 |
& krec, mythid ) |
970 |
if (using_ers) |
971 |
& call cost_sla_read( ersfi, erssta, ersper, |
972 |
& zeroRL, zeroRL, |
973 |
& ersob, ersma, |
974 |
& krec, mythid ) |
975 |
if (using_gfo) |
976 |
& call cost_sla_read( gfofi, gfosta, gfoper, |
977 |
& zeroRL, zeroRL, |
978 |
& gfoob, gfoma, |
979 |
& krec, mythid ) |
980 |
|
981 |
do bj = jtlo,jthi |
982 |
do bi = itlo,ithi |
983 |
do j = jmin,jmax |
984 |
do i = imin,imax |
985 |
if ( tpma(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj) |
986 |
& .NE.0. ) then |
987 |
anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+ |
988 |
& etaday(i,j,bi,bj)-tpob(i,j,bi,bj) |
989 |
& -mean_psMslaobs(i,j,bi,bj) |
990 |
anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+ |
991 |
& tpob(i,j,bi,bj) |
992 |
anom_psMslaobs_NUM(i,j,bi,bj)= |
993 |
& anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0 |
994 |
if ( (2*jrec.EQ.ndaysave).OR.(2*jrec-1.EQ.ndaysave) ) |
995 |
& anom_psMslaobs_MSK(i,j,bi,bj) = 1. _d 0 |
996 |
endif |
997 |
if ( ersma(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj) |
998 |
& .NE.0. ) then |
999 |
anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+ |
1000 |
& etaday(i,j,bi,bj)-ersob(i,j,bi,bj) |
1001 |
& -mean_psMslaobs(i,j,bi,bj) |
1002 |
anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+ |
1003 |
& ersob(i,j,bi,bj) |
1004 |
anom_psMslaobs_NUM(i,j,bi,bj)= |
1005 |
& anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0 |
1006 |
if ( (2*jrec.EQ.ndaysave).OR.(2*jrec-1.EQ.ndaysave) ) |
1007 |
& anom_psMslaobs_MSK(i,j,bi,bj) = 1. _d 0 |
1008 |
endif |
1009 |
if ( gfoma(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj) |
1010 |
& .NE.0. ) then |
1011 |
anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+ |
1012 |
& etaday(i,j,bi,bj)-gfoob(i,j,bi,bj) |
1013 |
& -mean_psMslaobs(i,j,bi,bj) |
1014 |
anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+ |
1015 |
& gfoob(i,j,bi,bj) |
1016 |
anom_psMslaobs_NUM(i,j,bi,bj)= |
1017 |
& anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0 |
1018 |
if ( (2*jrec.EQ.ndaysave).OR.(2*jrec-1.EQ.ndaysave) ) |
1019 |
& anom_psMslaobs_MSK(i,j,bi,bj) = 1. _d 0 |
1020 |
endif |
1021 |
enddo |
1022 |
enddo |
1023 |
enddo |
1024 |
enddo |
1025 |
|
1026 |
enddo !do jrec=1,ndaysave |
1027 |
|
1028 |
do bj = jtlo,jthi |
1029 |
do bi = itlo,ithi |
1030 |
do j = jmin,jmax |
1031 |
do i = imin,imax |
1032 |
if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND. |
1033 |
& ( maskc(i,j,1,bi,bj) .NE. 0. ) |
1034 |
& ) then |
1035 |
anom_psMslaobs(i,j,bi,bj) = |
1036 |
& anom_psMslaobs(i,j,bi,bj) / |
1037 |
& anom_psMslaobs_NUM(i,j,bi,bj) |
1038 |
anom_slaobs(i,j,bi,bj) = |
1039 |
& anom_slaobs(i,j,bi,bj) / |
1040 |
& anom_psMslaobs_NUM(i,j,bi,bj) |
1041 |
else |
1042 |
anom_psMslaobs(i,j,bi,bj) = 0. _d 0 |
1043 |
anom_slaobs(i,j,bi,bj) = 0. _d 0 |
1044 |
endif |
1045 |
enddo |
1046 |
enddo |
1047 |
enddo |
1048 |
enddo |
1049 |
|
1050 |
c PART 4.11: separate time variable offset |
1051 |
c ---------------------------------------- |
1052 |
|
1053 |
slaoffset = 0. _d 0 |
1054 |
slaoffset_sum = 0. _d 0 |
1055 |
slaoffset_weight = 0. _d 0 |
1056 |
|
1057 |
do bj = jtlo,jthi |
1058 |
do bi = itlo,ithi |
1059 |
do j = jmin,jmax |
1060 |
do i = imin,imax |
1061 |
if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND. |
1062 |
& ( maskc(i,j,1,bi,bj) .NE. 0. ) |
1063 |
& ) then |
1064 |
slaoffset=slaoffset+RA(i,j,bi,bj)* |
1065 |
& anom_psMslaobs(i,j,bi,bj) |
1066 |
slaoffset_sum=slaoffset_sum+RA(i,j,bi,bj) |
1067 |
c Of interest is the total weight for an average of N indep sample (not the area weighted average of weight) |
1068 |
c Since slaoffset is itself area weighted, however, the total weight is only approx the simple weight sum : |
1069 |
slaoffset_weight=slaoffset_weight+ |
1070 |
& gencost_weight(i,j,bi,bj,igen_lsc) |
1071 |
endif |
1072 |
enddo |
1073 |
enddo |
1074 |
enddo |
1075 |
enddo |
1076 |
|
1077 |
_GLOBAL_SUM_RL( slaoffset , mythid ) |
1078 |
_GLOBAL_SUM_RL( slaoffset_weight , mythid ) |
1079 |
_GLOBAL_SUM_RL( slaoffset_sum , mythid ) |
1080 |
|
1081 |
if (slaoffset_sum .eq. 0.0) then |
1082 |
if (gencost_outputlevel(igen_gmsl).GT.0) then |
1083 |
_BEGIN_MASTER( mythid ) |
1084 |
write(msgbuf,'(a)') ' sshv4: slaoffset_sum = zero!' |
1085 |
call print_message( msgbuf, standardmessageunit, |
1086 |
& SQUEEZE_RIGHT , mythid) |
1087 |
_END_MASTER( mythid ) |
1088 |
endif |
1089 |
c stop ' ... stopped in cost_ssh.' |
1090 |
else |
1091 |
if (gencost_outputlevel(igen_gmsl).GT.0) then |
1092 |
_BEGIN_MASTER( mythid ) |
1093 |
write(msgbuf,'(a,3d22.15)') ' sshv4:slaoffset,sum,weight=', |
1094 |
& slaoffset,slaoffset_sum,slaoffset_weight |
1095 |
call print_message( msgbuf, standardmessageunit, |
1096 |
& SQUEEZE_RIGHT , mythid) |
1097 |
_END_MASTER( mythid ) |
1098 |
endif |
1099 |
endif |
1100 |
|
1101 |
c-- Compute (average) slaoffset |
1102 |
slaoffset = slaoffset / slaoffset_sum |
1103 |
|
1104 |
c-- Subtract slaoffset from anom_psMslaobs |
1105 |
do bj = jtlo,jthi |
1106 |
do bi = itlo,ithi |
1107 |
do j = jmin,jmax |
1108 |
do i = imin,imax |
1109 |
if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND. |
1110 |
& ( maskc(i,j,1,bi,bj) .NE. 0. ) |
1111 |
& ) then |
1112 |
anom_psMslaobs(i,j,bi,bj) = |
1113 |
& anom_psMslaobs(i,j,bi,bj) - slaoffset |
1114 |
anom_slaobs(i,j,bi,bj) = |
1115 |
& anom_slaobs(i,j,bi,bj) - slaoffset |
1116 |
else |
1117 |
anom_psMslaobs(i,j,bi,bj) = 0. _d 0 |
1118 |
anom_slaobs(i,j,bi,bj) = 0. _d 0 |
1119 |
endif |
1120 |
enddo |
1121 |
enddo |
1122 |
enddo |
1123 |
enddo |
1124 |
|
1125 |
c PART 4.2: smooth anom_psMslaobs in space |
1126 |
c ---------------------------------------- |
1127 |
|
1128 |
if (gencost_outputlevel(igen_lsc).GT.0) then |
1129 |
write(fname4test(1:80),'(1a)') 'sladiff_raw' |
1130 |
call mdswritefield(fname4test,32,.false.,'RL', |
1131 |
& 1,anom_psMslaobs,irec,1,mythid) |
1132 |
|
1133 |
write(fname4test(1:80),'(1a)') 'slaobs_raw' |
1134 |
call mdswritefield(fname4test,32,.false.,'RL', |
1135 |
& 1,anom_slaobs,irec,1,mythid) |
1136 |
endif |
1137 |
|
1138 |
#ifdef ALLOW_SMOOTH |
1139 |
if ( useSMOOTH.AND.(k2_lsc.GT.0) ) |
1140 |
& call smooth_hetero2d(anom_psMslaobs,maskc, |
1141 |
& gencost_posproc_c(k2_lsc,igen_lsc), |
1142 |
& gencost_posproc_i(k2_lsc,igen_lsc),mythid) |
1143 |
#endif |
1144 |
|
1145 |
if (gencost_outputlevel(igen_lsc).GT.0) then |
1146 |
#ifdef ALLOW_SMOOTH |
1147 |
if ( useSMOOTH.AND.(k2_lsc.GT.0) ) |
1148 |
& call smooth_hetero2d(anom_slaobs,maskc, |
1149 |
& gencost_posproc_c(k2_lsc,igen_lsc), |
1150 |
& gencost_posproc_i(k2_lsc,igen_lsc),mythid) |
1151 |
#endif |
1152 |
c |
1153 |
write(fname4test(1:80),'(1a)') 'sladiff_smooth' |
1154 |
call mdswritefield(fname4test,32,.false.,'RL', |
1155 |
& 1,anom_psMslaobs,irec,1,mythid) |
1156 |
c |
1157 |
write(fname4test(1:80),'(1a)') 'slaobs_smooth' |
1158 |
call mdswritefield(fname4test,32,.false.,'RL', |
1159 |
& 1,anom_slaobs,irec,1,mythid) |
1160 |
endif |
1161 |
|
1162 |
c PART 4.3: compute cost function term |
1163 |
c ------------------------------------ |
1164 |
|
1165 |
do bj = jtlo,jthi |
1166 |
do bi = itlo,ithi |
1167 |
do j = jmin,jmax |
1168 |
do i = imin,imax |
1169 |
if ( (gencost_weight(i,j,bi,bj,igen_lsc).GT.0.).AND. |
1170 |
& (anom_psMslaobs_MSK(i,j,bi,bj).GT.0.) ) then |
1171 |
junk = anom_psMslaobs(i,j,bi,bj) |
1172 |
anom_slaobs_SUB(i,j,bi,bj) = anom_slaobs(i,j,bi,bj) |
1173 |
anom_psMslaobs_SUB(i,j,bi,bj) = anom_psMslaobs(i,j,bi,bj) |
1174 |
else |
1175 |
junk = 0. _d 0 |
1176 |
anom_slaobs_SUB(i,j,bi,bj)= 0. _d 0 |
1177 |
anom_psMslaobs_SUB(i,j,bi,bj)= 0. _d 0 |
1178 |
endif |
1179 |
objf_gencost(bi,bj,igen_lsc) = objf_gencost(bi,bj,igen_lsc) |
1180 |
& + junk*junk*gencost_weight(i,j,bi,bj,igen_lsc) |
1181 |
if ( (gencost_weight(i,j,bi,bj,igen_lsc).GT.0.).AND. |
1182 |
& (anom_psMslaobs_MSK(i,j,bi,bj).GT.0.) ) |
1183 |
& num_gencost(bi,bj,igen_lsc) = |
1184 |
& num_gencost(bi,bj,igen_lsc) + 1. _d 0 |
1185 |
enddo |
1186 |
enddo |
1187 |
enddo |
1188 |
enddo |
1189 |
|
1190 |
if (gencost_outputlevel(igen_lsc).GT.0) then |
1191 |
write(fname4test(1:80),'(1a)') 'sladiff_subsample' |
1192 |
call mdswritefield(fname4test,32,.false.,'RL', |
1193 |
& 1,anom_psMslaobs_SUB,irec,1,mythid) |
1194 |
c |
1195 |
write(fname4test(1:80),'(1a)') 'slaobs_subsample' |
1196 |
call mdswritefield(fname4test,32,.false.,'RL', |
1197 |
& 1,anom_slaobs_SUB,irec,1,mythid) |
1198 |
endif |
1199 |
|
1200 |
c PART 4.4: compute cost function term for global mean sea level |
1201 |
c -------------------------------------------------------------- |
1202 |
|
1203 |
IF ( ( MASTER_CPU_THREAD(myThid).AND.(igen_gmsl.NE.0) ).AND. |
1204 |
& ( slaoffset_sum.GT.0.25 _d 0 * globalArea ) ) THEN |
1205 |
junk=slaoffset |
1206 |
junkweight=1. _d 0 |
1207 |
objf_gencost(1,1,igen_gmsl) = objf_gencost(1,1,igen_gmsl) |
1208 |
& + junk*junk*junkweight |
1209 |
num_gencost(1,1,igen_gmsl) = num_gencost(1,1,igen_gmsl) |
1210 |
& + 1. _d 0 |
1211 |
c print*,'igen_gmsl',igen_gmsl, |
1212 |
ENDIF |
1213 |
|
1214 |
cgf ======================================================= |
1215 |
cgf PART 5: compute raw SLA cost term |
1216 |
cgf ======================================================= |
1217 |
|
1218 |
c time associated with this ndaysrec mean |
1219 |
krec = irec + (ndaysave/2) |
1220 |
|
1221 |
#ifdef ALLOW_AUTODIFF |
1222 |
call active_read_xy( fname, etaday, krec, doglobalread, |
1223 |
& ladinit, eccoiter, mythid, |
1224 |
& gencost_dummy(igen_etaday) ) |
1225 |
#else |
1226 |
CALL READ_REC_XY_RL( fname, etaday, kRec, 1, myThid ) |
1227 |
#endif |
1228 |
|
1229 |
if (.NOT.useEtaMean) |
1230 |
& CALL REMOVE_MEAN_RL( 1, etaday, maskInC, maskInC, rA, drF, |
1231 |
& 'etaday', 0. _d 0, myThid ) |
1232 |
|
1233 |
do bj = jtlo,jthi |
1234 |
do bi = itlo,ithi |
1235 |
do j = jmin,jmax |
1236 |
do i = imin,imax |
1237 |
tpma(i,j,bi,bj) = 0. _d 0 |
1238 |
tpob(i,j,bi,bj) = 0. _d 0 |
1239 |
ersma(i,j,bi,bj) = 0. _d 0 |
1240 |
ersob(i,j,bi,bj) = 0. _d 0 |
1241 |
gfoma(i,j,bi,bj) = 0. _d 0 |
1242 |
gfoob(i,j,bi,bj) = 0. _d 0 |
1243 |
enddo |
1244 |
enddo |
1245 |
enddo |
1246 |
enddo |
1247 |
|
1248 |
if (using_tpj) |
1249 |
& call cost_sla_read( tpfi, tpsta, tpper, |
1250 |
& zeroRL, zeroRL, |
1251 |
& tpob, tpma, |
1252 |
& krec, mythid ) |
1253 |
if (using_ers) |
1254 |
& call cost_sla_read( ersfi, erssta, ersper, |
1255 |
& zeroRL, zeroRL, |
1256 |
& ersob, ersma, |
1257 |
& krec, mythid ) |
1258 |
if (using_gfo) |
1259 |
& call cost_sla_read( gfofi, gfosta, gfoper, |
1260 |
& zeroRL, zeroRL, |
1261 |
& gfoob, gfoma, |
1262 |
& krec, mythid ) |
1263 |
|
1264 |
do bj = jtlo,jthi |
1265 |
do bi = itlo,ithi |
1266 |
do j = jmin,jmax |
1267 |
do i = imin,imax |
1268 |
anom_psMtpobs(i,j,bi,bj) = 0. _d 0 |
1269 |
anom_psMersobs(i,j,bi,bj) = 0. _d 0 |
1270 |
anom_psMgfoobs(i,j,bi,bj) = 0. _d 0 |
1271 |
anom_tpobs(i,j,bi,bj) = 0. _d 0 |
1272 |
anom_ersobs(i,j,bi,bj) = 0. _d 0 |
1273 |
anom_gfoobs(i,j,bi,bj) = 0. _d 0 |
1274 |
enddo |
1275 |
enddo |
1276 |
enddo |
1277 |
enddo |
1278 |
|
1279 |
do bj = jtlo,jthi |
1280 |
do bi = itlo,ithi |
1281 |
do j = jmin,jmax |
1282 |
do i = imin,imax |
1283 |
if ( tpma(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj).NE.0. ) |
1284 |
& then |
1285 |
anom_psMtpobs(i,j,bi,bj)= |
1286 |
& etaday(i,j,bi,bj) - tpob(i,j,bi,bj) |
1287 |
& - mean_psMslaobs(i,j,bi,bj) - slaoffset |
1288 |
& - anom_psMslaobs(i,j,bi,bj) |
1289 |
anom_tpobs(i,j,bi,bj)=tpob(i,j,bi,bj) - slaoffset |
1290 |
& - anom_slaobs(i,j,bi,bj) |
1291 |
endif |
1292 |
if ( ersma(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj).NE.0. ) |
1293 |
& then |
1294 |
anom_psMersobs(i,j,bi,bj)= |
1295 |
& etaday(i,j,bi,bj) - ersob(i,j,bi,bj) |
1296 |
& - mean_psMslaobs(i,j,bi,bj) - slaoffset |
1297 |
& - anom_psMslaobs(i,j,bi,bj) |
1298 |
anom_ersobs(i,j,bi,bj)=ersob(i,j,bi,bj) - slaoffset |
1299 |
& - anom_slaobs(i,j,bi,bj) |
1300 |
endif |
1301 |
if ( gfoma(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj).NE.0. ) |
1302 |
& then |
1303 |
anom_psMgfoobs(i,j,bi,bj)= |
1304 |
& etaday(i,j,bi,bj) - gfoob(i,j,bi,bj) |
1305 |
& - mean_psMslaobs(i,j,bi,bj) - slaoffset |
1306 |
& - anom_psMslaobs(i,j,bi,bj) |
1307 |
anom_gfoobs(i,j,bi,bj)=gfoob(i,j,bi,bj) - slaoffset |
1308 |
& - anom_slaobs(i,j,bi,bj) |
1309 |
endif |
1310 |
enddo |
1311 |
enddo |
1312 |
enddo |
1313 |
enddo |
1314 |
|
1315 |
if (gencost_outputlevel(igen_tp).GT.0) then |
1316 |
write(fname4test(1:80),'(1a)') 'sladiff_tp_raw' |
1317 |
call mdswritefield(fname4test,32,.false.,'RL', |
1318 |
& 1,anom_psMtpobs,irec,1,mythid) |
1319 |
write(fname4test(1:80),'(1a)') 'slaobs_tp_raw' |
1320 |
call mdswritefield(fname4test,32,.false.,'RL', |
1321 |
& 1,anom_tpobs,irec,1,mythid) |
1322 |
endif |
1323 |
|
1324 |
if (gencost_outputlevel(igen_ers).GT.0) then |
1325 |
write(fname4test(1:80),'(1a)') 'sladiff_ers_raw' |
1326 |
call mdswritefield(fname4test,32,.false.,'RL', |
1327 |
& 1,anom_psMersobs,irec,1,mythid) |
1328 |
write(fname4test(1:80),'(1a)') 'slaobs_ers_raw' |
1329 |
call mdswritefield(fname4test,32,.false.,'RL', |
1330 |
& 1,anom_ersobs,irec,1,mythid) |
1331 |
endif |
1332 |
|
1333 |
if (gencost_outputlevel(igen_gfo).GT.0) then |
1334 |
write(fname4test(1:80),'(1a)') 'sladiff_gfo_raw' |
1335 |
call mdswritefield(fname4test,32,.false.,'RL', |
1336 |
& 1,anom_psMgfoobs,irec,1,mythid) |
1337 |
write(fname4test(1:80),'(1a)') 'slaobs_gfo_raw' |
1338 |
call mdswritefield(fname4test,32,.false.,'RL', |
1339 |
& 1,anom_gfoobs,irec,1,mythid) |
1340 |
endif |
1341 |
|
1342 |
do bj = jtlo,jthi |
1343 |
do bi = itlo,ithi |
1344 |
do j = jmin,jmax |
1345 |
do i = imin,imax |
1346 |
c-- The array psobs contains SSH anomalies. |
1347 |
junkweight = mean_psMslaobs_MSK(i,j,bi,bj) |
1348 |
& *gencost_weight(i,j,bi,bj,igen_tp) |
1349 |
& *tpma(i,j,bi,bj) |
1350 |
junk = anom_psMtpobs(i,j,bi,bj) |
1351 |
objf_gencost(bi,bj,igen_tp) = |
1352 |
& objf_gencost(bi,bj,igen_tp)+junk*junk*junkweight |
1353 |
if ( junkweight .ne. 0. ) |
1354 |
& num_gencost(bi,bj,igen_tp) = |
1355 |
& num_gencost(bi,bj,igen_tp) + 1. _d 0 |
1356 |
c-- The array ersobs contains SSH anomalies. |
1357 |
junkweight = mean_psMslaobs_MSK(i,j,bi,bj) |
1358 |
& *gencost_weight(i,j,bi,bj,igen_ers) |
1359 |
& *ersma(i,j,bi,bj) |
1360 |
junk = anom_psMersobs(i,j,bi,bj) |
1361 |
objf_gencost(bi,bj,igen_ers) = |
1362 |
& objf_gencost(bi,bj,igen_ers)+junk*junk*junkweight |
1363 |
if ( junkweight .ne. 0. ) |
1364 |
& num_gencost(bi,bj,igen_ers) = |
1365 |
& num_gencost(bi,bj,igen_ers) + 1. _d 0 |
1366 |
c-- The array gfoobs contains SSH anomalies. |
1367 |
junkweight = mean_psMslaobs_MSK(i,j,bi,bj) |
1368 |
& *gencost_weight(i,j,bi,bj,igen_gfo) |
1369 |
& *gfoma(i,j,bi,bj) |
1370 |
junk = anom_psMgfoobs(i,j,bi,bj) |
1371 |
objf_gencost(bi,bj,igen_gfo) = |
1372 |
& objf_gencost(bi,bj,igen_gfo)+junk*junk*junkweight |
1373 |
if ( junkweight .ne. 0. ) |
1374 |
& num_gencost(bi,bj,igen_gfo) = |
1375 |
& num_gencost(bi,bj,igen_gfo) + 1. _d 0 |
1376 |
enddo |
1377 |
enddo |
1378 |
enddo |
1379 |
enddo |
1380 |
|
1381 |
enddo |
1382 |
|
1383 |
#endif /* ifdef ALLOW_GENCOST_CONTRIBUTION */ |
1384 |
|
1385 |
end |