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

Annotation of /MITgcm/pkg/ecco/cost_gencost_sshv4.F

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

Revision 1.4 - (hide annotations) (download)
Thu Feb 23 23:31:48 2012 UTC (12 years, 2 months ago) by gforget
Branch: MAIN
Changes since 1.3: +40 -3 lines 
- dont create ad files unless adjoint run.
- bug fix : sla_startdate etc. rather than modelstartdate.
- ecco_cost_driver.F : dont do cost_ssh.F when cost_gencost_sshv4.F.
- cost_gencost_sshv4.F : only re-reference MDT if altimeter period is covered.
1 gforget 1.4 C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_gencost_sshv4.F,v 1.3 2011/05/18 03:15:38 gforget Exp $
2 gforget 1.1 C $Name: $
3    
4     #include "COST_CPPOPTIONS.h"
5    
6    
7     subroutine cost_gencost_sshv4(
8     I myiter,
9     I mytime,
10     I mythid
11     & )
12    
13     c ==================================================================
14     c SUBROUTINE cost_gencost_sshv4
15     c ==================================================================
16     c
17     c o Evaluate cost function contributions of sea surface height.
18     c
19     c started: Gael Forget, Oct-2009
20     c
21     c working assumption for the time mean dynamic topography (MDT) constraint:
22     c the various SLA data sets (tp, ers, gfo) have been consistenty cross-referenced,
23     c as done in the RADS data sets. We do not need to know the reference dynamic
24     c topography (or SSH/Geoid). Simply it is assumed that the biases
25     c between instruments have been taken care of. This is only a matter
26     c for the MDT constraint, not for SLA constraints (see below).
27     c
28     cgf 1) there are a few hardcoded numbers that will eventually be put in common
29     cgf blocks/namelists
30     cgf 2) there are a several refinements that should be considered, such as
31     cgf modulating weights with respect to numbers of samples
32     c
33     c ==================================================================
34     c SUBROUTINE cost_gencost_sshv4
35     c ==================================================================
36    
37     implicit none
38    
39     c == global variables ==
40    
41     #include "EEPARAMS.h"
42     #include "SIZE.h"
43     #include "PARAMS.h"
44     #include "GRID.h"
45    
46     #include "ecco_cost.h"
47     #include "ctrl.h"
48     #include "ctrl_dummy.h"
49     #include "optim.h"
50     #include "DYNVARS.h"
51     #ifdef ALLOW_PROFILES
52     #include "profiles.h"
53     #endif
54 gforget 1.4 #include "cal.h"
55 gforget 1.1
56     c == routine arguments ==
57    
58     integer myiter
59     _RL mytime
60     integer mythid
61    
62     #ifdef ALLOW_SSH_COST_CONTRIBUTION
63     #ifdef ALLOW_SMOOTH
64     #ifdef ALLOW_GENCOST_CONTRIBUTION
65     #ifdef ALLOW_GENCOST_FREEFORM
66     c == local variables ==
67    
68     integer bi,bj
69     integer i,j,k
70     integer itlo,ithi
71     integer jtlo,jthi
72     integer jmin,jmax
73     integer imin,imax
74     integer irec,jrec,krec
75     integer ilps
76     integer gwunit
77    
78     logical doglobalread
79     logical ladinit
80    
81     c mapping to gencost
82     integer igen_mdt, igen_lsc
83     integer igen_tp, igen_ers, igen_gfo
84    
85     _RL offset,fac
86     _RL offset_sum
87     _RL psmean ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
88     _RL diagnosfld ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
89    
90     c for PART 1: re-reference MDT (tpmean) to the inferred SLA reference field
91     _RL mean_slaobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
92     _RL mean_slaobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
93    
94     c for PART 2: compute time mean differences over the model period
95     _RL mean_slaobs2(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
96     _RL mean_psMssh_all(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
97     _RL mean_psMssh_all_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
98     _RL mean_psMssh_all_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
99    
100     _RL mean_psMtpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
101     _RL mean_psMtpobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
102     _RL mean_psMtpobs_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
103    
104     _RL mean_psMersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
105     _RL mean_psMersobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
106     _RL mean_psMersobs_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
107    
108     _RL mean_psMgfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
109     _RL mean_psMgfoobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
110     _RL mean_psMgfoobs_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
111    
112     c for PART 4/5: compute smooth/raw anomalies
113     _RL anom_psMslaobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
114     _RL anom_slaobs (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
115     _RL anom_psMslaobs_NUM (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
116    
117     _RL anom_psMtpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
118     _RL anom_psMtpobs_NUM (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
119     _RL anom_tpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
120    
121     _RL anom_psMersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
122     _RL anom_psMersobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
123     _RL anom_ersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
124    
125     _RL anom_psMgfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
126     _RL anom_psMgfoobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
127     _RL anom_gfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
128    
129     integer tpmean_y0,tpmean_y1,year,day
130     integer num_var
131    
132     _RL junk,junkweight
133    
134     integer ndaysave
135     _RL ndaysaveRL
136    
137     character*(80) fname
138     character*(80) fname4test
139     character*(MAX_LEN_MBUF) msgbuf
140    
141 gforget 1.4 LOGICAL doReference
142    
143 gforget 1.1 c == external functions ==
144    
145     integer ilnblnk
146     external ilnblnk
147    
148     c == end of interface ==
149    
150     jtlo = mybylo(mythid)
151     jthi = mybyhi(mythid)
152     itlo = mybxlo(mythid)
153     ithi = mybxhi(mythid)
154     jmin = 1
155     jmax = sny
156     imin = 1
157     imax = snx
158    
159     c-- detect the relevant gencost indices
160     igen_mdt=0
161     igen_tp =0
162     igen_ers=0
163     igen_gfo=0
164     igen_lsc=0
165     do k=1,NGENCOST
166     if (gencost_name(k).EQ.'sshv4-mdt') igen_mdt=k
167     if (gencost_name(k).EQ.'sshv4-tp') igen_tp=k
168     if (gencost_name(k).EQ.'sshv4-ers') igen_ers=k
169     if (gencost_name(k).EQ.'sshv4-gfo') igen_gfo=k
170     if (gencost_name(k).EQ.'sshv4-lsc') igen_lsc=k
171     enddo
172    
173     c-- First, read tiled data.
174     doglobalread = .false.
175     ladinit = .false.
176    
177     write(fname(1:80),'(80a)') ' '
178     ilps=ilnblnk( psbarfile )
179     write(fname(1:80),'(2a,i10.10)')
180     & psbarfile(1:ilps),'.',optimcycle
181    
182    
183     cgf =======================================================
184     cgf PART 1:
185     cgf x Get the MDT (tpmean) ... compute the sample mean
186     cgf (mean_slaobs) of the SLA data (i.e. RADS for tp, ers, and gfo
187     cgf together) over the time interval of the MDT ... subtract
188     cgf mean_slaobs from tpmean.
189     cgf x At this point, tpmean is the inferred SLA reference field.
190     cgf x From there on, tpmean+sla will be directly comparable to
191     cgf the model SSH (psbar).
192     cgf =======================================================
193    
194     c-- Read mean field and mask
195     call cost_ReadTopexMean( mythid )
196    
197     c-- Compute mean_slaobs: sample mean SLA over the time period of tpmean.
198    
199     c pavlis and ecco/rio
200     tpmean_y0=1993
201     tpmean_y1=2004
202     c maximenko
203     c tpmean_y0=1992
204     c tpmean_y1=2002
205     c rio
206     c tpmean_y0=1993
207     c tpmean_y1=1999
208    
209 gforget 1.4 doReference=.FALSE.
210     if ((modelstartdate(1).GT.1992*100).AND.
211     & (modelstartdate(1).LT.2011*100).AND.
212     & (ndaysrec.GE.365)) doReference=.TRUE.
213    
214     write(msgbuf,'(a,l)') ' sshv4:re-reference MDT=',doReference
215     call print_message( msgbuf, standardmessageunit,
216     & SQUEEZE_RIGHT , mythid)
217    
218 gforget 1.1 do bj = jtlo,jthi
219     do bi = itlo,ithi
220     do j = jmin,jmax
221     do i = imin,imax
222     mean_slaobs(i,j,bi,bj) = 0. _d 0
223     mean_slaobs_NUM(i,j,bi,bj) = 0. _d 0
224     enddo
225     enddo
226     enddo
227     enddo
228    
229     do year=tpmean_y0,tpmean_y1
230     do day=1,366
231     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
232     call cost_sla_read_yd( topexfile,
233     & tpobs, tpmask,
234     & year, day, mythid )
235     #endif
236     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
237     call cost_sla_read_yd( ersfile,
238     & ersobs, ersmask,
239     & year, day, mythid )
240     #endif
241     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
242     call cost_sla_read_yd( gfofile,
243     & gfoobs, gfomask,
244     & year, day, mythid )
245     #endif
246     do bj = jtlo,jthi
247     do bi = itlo,ithi
248     do j = jmin,jmax
249     do i = imin,imax
250     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
251     if ( tpmask(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)*
252     & gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then
253     mean_slaobs(i,j,bi,bj)= mean_slaobs(i,j,bi,bj)+
254     & tpobs(i,j,bi,bj)
255     mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
256     endif
257     #endif
258     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
259     if ( ersmask(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)*
260     & gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then
261     mean_slaobs(i,j,bi,bj)= mean_slaobs(i,j,bi,bj)+
262     & ersobs(i,j,bi,bj)
263     mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
264     endif
265     #endif
266     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
267     if ( gfomask(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)*
268     & gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then
269     mean_slaobs(i,j,bi,bj)= mean_slaobs(i,j,bi,bj)+
270     & gfoobs(i,j,bi,bj)
271     mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
272     endif
273     #endif
274     enddo
275     enddo
276     enddo
277     enddo
278    
279     enddo !do day=1,366
280     enddo !do year=tpmean_y0,tpmean_y1
281    
282     do bj = jtlo,jthi
283     do bi = itlo,ithi
284     do j = jmin,jmax
285     do i = imin,imax
286     if ( ( mean_slaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
287     & ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
288     & ( abs(YC(i,j,bi,bj)) .LE. 66. ).AND.
289     & ( tpmeanmask(i,j,bi,bj) .NE. 0. ) ) then
290     mean_slaobs(i,j,bi,bj) = mean_slaobs(i,j,bi,bj) /
291     & mean_slaobs_NUM(i,j,bi,bj)
292     else
293     mean_slaobs(i,j,bi,bj) = 0. _d 0
294     endif
295     enddo
296     enddo
297     enddo
298     enddo
299    
300    
301     c-- smooth mean_slaobs:
302    
303     write(fname4test(1:80),'(1a)') 'sla2mdt_raw'
304     call mdswritefield(fname4test,32,.false.,'RL',
305     & 1,mean_slaobs,1,1,mythid)
306    
307     call smooth_hetero2d(mean_slaobs,maskc,
308     & gencost_scalefile(igen_mdt),300,mythid)
309    
310     write(fname4test(1:80),'(1a)') 'sla2mdt_smooth'
311     call mdswritefield(fname4test,32,.false.,'RL',
312     & 1,mean_slaobs,1,1,mythid)
313    
314     c-- re-reference tpmean:
315     do bj = jtlo,jthi
316     do bi = itlo,ithi
317     do j = jmin,jmax
318     do i = imin,imax
319     if ( ( tpmeanmask(i,j,bi,bj) .NE. 0. ).AND.
320 gforget 1.4 & ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
321     & ( doReference ) ) then
322 gforget 1.1 tpmean(i,j,bi,bj) = tpmean(i,j,bi,bj)
323     & -mean_slaobs(i,j,bi,bj)
324     endif
325     enddo
326     enddo
327     enddo
328     enddo
329    
330    
331     cgf =======================================================
332     cgf PART 2: compute sample means of psbar-slaobs over the
333     cgf period that is covered by the model (i.e. psbar).
334     cgf x for all SLA data sets together: mean_psMssh_all, mean_psMssh_all_MSK,
335     cgf and offset will be used in PART 3 (MDT cost term).
336     cgf x for each SLA data individually. mean_psMtpobs, mean_psMtpobs_MS, etc.
337     cgf will be used in PART 4&5 (SLA cost terms).
338     cgf =======================================================
339    
340     do bj = jtlo,jthi
341     do bi = itlo,ithi
342     do j = jmin,jmax
343     do i = imin,imax
344    
345     psmean(i,j,bi,bj) = 0. _d 0
346     mean_psMtpobs(i,j,bi,bj) = 0. _d 0
347     mean_psMersobs(i,j,bi,bj) = 0. _d 0
348     mean_psMgfoobs(i,j,bi,bj) = 0. _d 0
349     mean_psMssh_all(i,j,bi,bj) = 0. _d 0
350     mean_slaobs2(i,j,bi,bj) = 0. _d 0
351    
352     mean_psMtpobs_NUM(i,j,bi,bj) = 0. _d 0
353     mean_psMersobs_NUM(i,j,bi,bj) = 0. _d 0
354     mean_psMgfoobs_NUM(i,j,bi,bj) = 0. _d 0
355     mean_psMssh_all_NUM(i,j,bi,bj) = 0. _d 0
356    
357     mean_psMtpobs_MSK(i,j,bi,bj) = 0. _d 0
358     mean_psMersobs_MSK(i,j,bi,bj) = 0. _d 0
359     mean_psMgfoobs_MSK(i,j,bi,bj) = 0. _d 0
360    
361     enddo
362     enddo
363     enddo
364     enddo
365     offset = 0. _d 0
366     offset_sum = 0. _d 0
367    
368    
369     do irec = 1, ndaysrec
370    
371     call active_read_xy( fname, psbar, irec, doglobalread,
372     & ladinit, optimcycle, mythid,
373     & xx_psbar_mean_dummy )
374    
375 gforget 1.3 #ifndef ALLOW_PSBAR_MEAN
376     CALL REMOVE_MEAN_RL( 1, psbar, maskInC, maskInC, rA, drF,
377     & 'psbar', myTime, myThid )
378     #endif
379    
380 gforget 1.1 #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
381     call cost_sla_read( topexfile, topexstartdate, topexperiod,
382     & topexintercept, topexslope,
383     & tpobs, tpmask,
384     & irec, mythid )
385     #endif
386     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
387     call cost_sla_read( ersfile, ersstartdate, ersperiod,
388     & ersintercept, ersslope,
389     & ersobs, ersmask,
390     & irec, mythid )
391     #endif
392     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
393     call cost_sla_read( gfofile, gfostartdate, gfoperiod,
394     & gfointercept, gfoslope,
395     & gfoobs, gfomask,
396     & irec, mythid )
397     #endif
398    
399     do bj = jtlo,jthi
400     do bi = itlo,ithi
401     do j = jmin,jmax
402     do i = imin,imax
403     psmean(i,j,bi,bj) = psmean(i,j,bi,bj) +
404     & psbar(i,j,bi,bj) / float(ndaysrec)
405     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
406     if ( tpmask(i,j,bi,bj)*
407     & gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then
408     mean_slaobs2(i,j,bi,bj)=
409     & mean_slaobs2(i,j,bi,bj)+tpobs(i,j,bi,bj)
410     mean_psMtpobs(i,j,bi,bj) =
411     & mean_psMtpobs(i,j,bi,bj) +
412     & psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)
413     mean_psMtpobs_NUM(i,j,bi,bj) =
414     & mean_psMtpobs_NUM(i,j,bi,bj) + 1. _d 0
415     endif
416     #endif
417     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
418     if ( ersmask(i,j,bi,bj)*
419     & gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then
420     mean_slaobs2(i,j,bi,bj)=
421     & mean_slaobs2(i,j,bi,bj)+ersobs(i,j,bi,bj)
422     mean_psMersobs(i,j,bi,bj) =
423     & mean_psMersobs(i,j,bi,bj) +
424     & psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)
425     mean_psMersobs_NUM(i,j,bi,bj) =
426     & mean_psMersobs_NUM(i,j,bi,bj) + 1. _d 0
427     endif
428     #endif
429     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
430     if ( gfomask(i,j,bi,bj)*
431     & gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then
432     mean_slaobs2(i,j,bi,bj)=
433     & mean_slaobs2(i,j,bi,bj)+gfoobs(i,j,bi,bj)
434     mean_psMgfoobs(i,j,bi,bj) =
435     & mean_psMgfoobs(i,j,bi,bj) +
436     & psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)
437     mean_psMgfoobs_NUM(i,j,bi,bj) =
438     & mean_psMgfoobs_NUM(i,j,bi,bj) + 1. _d 0
439     endif
440     #endif
441     enddo
442     enddo
443     enddo
444     enddo
445    
446     c-- END loop over records for the first time.
447     enddo
448    
449     do bj = jtlo,jthi
450     do bi = itlo,ithi
451     do j = jmin,jmax
452     do i = imin,imax
453     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
454 gforget 1.3 if ( ( mean_psMtpobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
455     & ( abs(YC(i,j,bi,bj)) .LE. 66. ) ) then
456 gforget 1.1 mean_psMssh_all(i,j,bi,bj) =
457     & mean_psMssh_all(i,j,bi,bj) +
458     & mean_psMtpobs(i,j,bi,bj)
459     mean_psMssh_all_NUM(i,j,bi,bj) =
460     & mean_psMssh_all_NUM(i,j,bi,bj) +
461     & mean_psMtpobs_NUM(i,j,bi,bj)
462     mean_psMtpobs(i,j,bi,bj) =
463     & mean_psMtpobs(i,j,bi,bj) /
464     & mean_psMtpobs_NUM(i,j,bi,bj)
465     mean_psMtpobs_MSK(i,j,bi,bj) = 1. _d 0
466     endif
467     #endif
468     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
469 gforget 1.3 if ( ( mean_psMersobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
470     & ( abs(YC(i,j,bi,bj)) .LE. 66. ) ) then
471 gforget 1.1 mean_psMssh_all(i,j,bi,bj) =
472     & mean_psMssh_all(i,j,bi,bj) +
473     & mean_psMersobs(i,j,bi,bj)
474     mean_psMssh_all_NUM(i,j,bi,bj) =
475     & mean_psMssh_all_NUM(i,j,bi,bj) +
476     & mean_psMersobs_NUM(i,j,bi,bj)
477     mean_psMersobs(i,j,bi,bj) =
478     & mean_psMersobs(i,j,bi,bj) /
479     & mean_psMersobs_NUM(i,j,bi,bj)
480     mean_psMersobs_MSK(i,j,bi,bj) = 1. _d 0
481     endif
482     #endif
483     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
484 gforget 1.3 if ( ( mean_psMgfoobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
485     & ( abs(YC(i,j,bi,bj)) .LE. 66. ) ) then
486 gforget 1.1 mean_psMssh_all(i,j,bi,bj) =
487     & mean_psMssh_all(i,j,bi,bj) +
488     & mean_psMgfoobs(i,j,bi,bj)
489     mean_psMssh_all_NUM(i,j,bi,bj) =
490     & mean_psMssh_all_NUM(i,j,bi,bj) +
491     & mean_psMgfoobs_NUM(i,j,bi,bj)
492     mean_psMgfoobs(i,j,bi,bj) =
493     & mean_psMgfoobs(i,j,bi,bj) /
494     & mean_psMgfoobs_NUM(i,j,bi,bj)
495     mean_psMgfoobs_MSK(i,j,bi,bj) = 1. _d 0
496     endif
497     #endif
498     if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .NE. 0. ).AND.
499     & ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
500     & ( abs(YC(i,j,bi,bj)) .LE. 66. ).AND.
501 gforget 1.4 & ( tpmeanmask(i,j,bi,bj) .NE. 0. ).AND.
502     & ( doReference ) ) then
503 gforget 1.1 mean_slaobs2(i,j,bi,bj) =
504     & mean_slaobs2(i,j,bi,bj) /
505     & mean_psMssh_all_NUM(i,j,bi,bj)
506     mean_psMssh_all(i,j,bi,bj) =
507     & mean_psMssh_all(i,j,bi,bj) /
508     & mean_psMssh_all_NUM(i,j,bi,bj)-tpmean(i,j,bi,bj)
509     mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
510     offset=offset+mean_psMssh_all(i,j,bi,bj)*
511     & mean_psMssh_all_NUM(i,j,bi,bj)
512     offset_sum=offset_sum+mean_psMssh_all_NUM(i,j,bi,bj)
513 gforget 1.4 elseif ( ( tpmeanmask(i,j,bi,bj) .NE. 0. ) .AND.
514     & ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
515     & ( .NOT.doReference ) ) then
516     mean_slaobs2(i,j,bi,bj) = 0.d0
517     mean_psMssh_all(i,j,bi,bj) = 0. _d 0
518     mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
519     offset=offset+RA(i,j,bi,bj)
520     & *( psmean(i,j,bi,bj) -tpmean(i,j,bi,bj) )
521     offset_sum=offset_sum+RA(i,j,bi,bj)
522 gforget 1.1 else
523     mean_slaobs2(i,j,bi,bj) = 0.d0
524     mean_psMssh_all(i,j,bi,bj) = 0. _d 0
525     mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0
526     endif
527     enddo
528     enddo
529     enddo
530     enddo
531    
532     c-- Do a global summation.
533     _GLOBAL_SUM_RL( offset , mythid )
534     _GLOBAL_SUM_RL( offset_sum , mythid )
535    
536 gforget 1.4 write(msgbuf,'(a,2d22.15)') ' sshv4:offset=',offset,offset_sum
537     call print_message( msgbuf, standardmessageunit,
538     & SQUEEZE_RIGHT , mythid)
539    
540 gforget 1.1 if (offset_sum .eq. 0.0) then
541     _BEGIN_MASTER( mythid )
542     write(msgbuf,'(a)') ' cost_ssh: offset_sum = zero!'
543     call print_message( msgbuf, standardmessageunit,
544     & SQUEEZE_RIGHT , mythid)
545     _END_MASTER( mythid )
546     c stop ' ... stopped in cost_ssh.'
547     else
548     _BEGIN_MASTER( mythid )
549     write(msgbuf,'(a,d22.15)')
550     & ' cost_ssh: offset_sum = ',offset_sum
551     call print_message( msgbuf, standardmessageunit,
552     & SQUEEZE_RIGHT , mythid)
553     _END_MASTER( mythid )
554     endif
555    
556     c-- Compute (average) offset
557     offset = offset / offset_sum
558    
559     c-- subtract offset from mean_psMssh_all and psmean
560     do bj = jtlo,jthi
561     do bi = itlo,ithi
562     do j = jmin,jmax
563     do i = imin,imax
564    
565     if ( (mean_psMssh_all_MSK(i,j,bi,bj) .NE. 0.) .AND.
566     & ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
567     & ( abs(YC(i,j,bi,bj)) .LE. 66. ) ) then
568     c use the re-referencing approach
569     mean_psMssh_all(i,j,bi,bj) =
570     & mean_psMssh_all(i,j,bi,bj) - offset
571     elseif ( ( tpmeanmask(i,j,bi,bj) .NE. 0. ) .AND.
572     & ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
573     & ( abs(YC(i,j,bi,bj)) .GT. 66. ) ) then
574     c use the simpler approach
575     mean_psMssh_all(i,j,bi,bj) =
576     & psmean(i,j,bi,bj) -tpmean(i,j,bi,bj) - offset
577     else
578     mean_psMssh_all(i,j,bi,bj) = 0. _d 0
579     endif
580    
581 gforget 1.4 c use the simpler approach
582     if ( ( .NOT.doReference ).AND.
583     & ( tpmeanmask(i,j,bi,bj) .NE. 0. ) .AND.
584     & ( maskc(i,j,1,bi,bj) .NE. 0. ) ) then
585     mean_psMssh_all(i,j,bi,bj) =
586     & psmean(i,j,bi,bj) -tpmean(i,j,bi,bj) - offset
587     else
588     mean_psMssh_all(i,j,bi,bj) = 0. _d 0
589     endif
590    
591 gforget 1.1 if ( maskc(i,j,1,bi,bj) .NE. 0. )
592     & psmean(i,j,bi,bj)=psmean(i,j,bi,bj)-offset
593     enddo
594     enddo
595     enddo
596     enddo
597    
598     c-- smooth mean_psMssh_all
599     write(fname4test(1:80),'(1a)') 'mdtdiff_raw'
600     call mdswritefield(fname4test,32,.false.,'RL',
601     & 1,mean_psMssh_all,1,1,mythid)
602    
603     call smooth_hetero2d(mean_psMssh_all,maskc,
604     & gencost_scalefile(igen_mdt),300,mythid)
605    
606     write(fname4test(1:80),'(1a)') 'mdtdiff_smooth'
607     call mdswritefield(fname4test,32,.false.,'RL',
608     & 1,mean_psMssh_all,1,1,mythid)
609    
610     write(fname4test(1:80),'(1a)') 'sla2model_raw'
611     call mdswritefield(fname4test,32,.false.,'RL',
612     & 1,mean_slaobs2,1,1,mythid)
613    
614     call smooth_hetero2d(mean_slaobs2,maskc,
615     & gencost_scalefile(igen_mdt),300,mythid)
616    
617     write(fname4test(1:80),'(1a)') 'sla2model_smooth'
618     call mdswritefield(fname4test,32,.false.,'RL',
619     & 1,mean_slaobs2,1,1,mythid)
620    
621     cgf at this point:
622     cgf 1) mean_psMssh_all is the sample mean <psbar-slaobs-tpmean-offset>,
623     cgf to which a smoothing filter has been applied.
624     cgf 2) mean_psMtpobs is the (unsmoothed) sample mean <psbar-tpobs>.
625     cgf And similarly for ers and gfo, each treated separately.
626    
627     #ifdef ALLOW_PROFILES
628     do bj = jtlo,jthi
629     do bi = itlo,ithi
630     do j = 1,sny
631     do i = 1,snx
632     prof_etan_mean(i,j,bi,bj)=psmean(i,j,bi,bj)
633     enddo
634     enddo
635     enddo
636     enddo
637     _EXCH_XY_RL( prof_etan_mean, mythid )
638     #endif
639    
640    
641     cgf =======================================================
642     cgf PART 3: compute MDT cost term
643     cgf =======================================================
644    
645    
646     #ifdef ALLOW_SSH_MEAN_COST_CONTRIBUTION
647    
648     do bj = jtlo,jthi
649     do bi = itlo,ithi
650     do j = jmin,jmax
651     do i = imin,imax
652     junk = mean_psMssh_all(i,j,bi,bj)
653     junkweight = gencost_weight(i,j,bi,bj,igen_mdt)
654     & *tpmeanmask(i,j,bi,bj)
655     objf_gencost(igen_mdt,bi,bj) = objf_gencost(igen_mdt,bi,bj)
656     & + junk*junk*junkweight
657     if ( junkweight .ne. 0. ) num_gencost(igen_mdt,bi,bj) =
658     & num_gencost(igen_mdt,bi,bj) + 1. _d 0
659     diagnosfld(i,j,bi,bj) = junk*junk*junkweight
660     enddo
661     enddo
662     enddo
663     enddo
664    
665     CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean', ' ', diagnosfld,
666     & optimcycle, mythid )
667    
668     #endif /* ALLOW_SSH_MEAN_COST_CONTRIBUTION */
669    
670    
671    
672     cgf =======================================================
673     cgf PART 4: compute smooth SLA cost term
674     cgf =======================================================
675    
676    
677     ndaysave=35
678     ndaysaveRL=ndaysave
679    
680 gforget 1.2 do irec = 1, ndaysrec-ndaysave+1, 7
681 gforget 1.1
682     do bj = jtlo,jthi
683     do bi = itlo,ithi
684     do j = jmin,jmax
685     do i = imin,imax
686     anom_psMslaobs(i,j,bi,bj) = 0. _d 0
687     anom_slaobs(i,j,bi,bj) = 0. _d 0
688     anom_psMtpobs(i,j,bi,bj) = 0. _d 0
689     anom_psMersobs(i,j,bi,bj) = 0. _d 0
690     anom_psMgfoobs(i,j,bi,bj) = 0. _d 0
691     anom_psMslaobs_NUM(i,j,bi,bj) = 0. _d 0
692     anom_psMtpobs_NUM(i,j,bi,bj) = 0. _d 0
693     anom_psMersobs_NUM(i,j,bi,bj) = 0. _d 0
694     anom_psMgfoobs_NUM(i,j,bi,bj) = 0. _d 0
695     enddo
696     enddo
697     enddo
698     enddo
699    
700     c PART 4.1: compute running sample average over ndaysave
701     c ------------------------------------------------------
702    
703     do jrec=1,ndaysave
704    
705     krec=irec+jrec-1
706    
707     call active_read_xy( fname, psbar, krec, doglobalread,
708     & ladinit, optimcycle, mythid,
709     & xx_psbar_mean_dummy )
710    
711 gforget 1.3 #ifndef ALLOW_PSBAR_MEAN
712     CALL REMOVE_MEAN_RL( 1, psbar, maskInC, maskInC, rA, drF,
713     & 'psbar', myTime, myThid )
714     #endif
715    
716 gforget 1.1 #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
717     call cost_sla_read( topexfile, topexstartdate, topexperiod,
718     & topexintercept, topexslope,
719     & tpobs, tpmask,
720     & krec, mythid )
721     #endif
722     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
723     call cost_sla_read( ersfile, ersstartdate, ersperiod,
724     & ersintercept, ersslope,
725     & ersobs, ersmask,
726     & krec, mythid )
727     #endif
728     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
729     call cost_sla_read( gfofile, gfostartdate, gfoperiod,
730     & gfointercept, gfoslope,
731     & gfoobs, gfomask,
732     & krec, mythid )
733     #endif
734    
735     do bj = jtlo,jthi
736     do bi = itlo,ithi
737     do j = jmin,jmax
738     do i = imin,imax
739     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
740     if ( tpmask(i,j,bi,bj)*mean_psMtpobs_MSK(i,j,bi,bj)
741     & .NE.0. ) then
742     anom_psMtpobs(i,j,bi,bj)= anom_psMtpobs(i,j,bi,bj)+
743     & psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)
744     & -mean_psMtpobs(i,j,bi,bj)
745     anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
746     & psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)
747     & -mean_psMtpobs(i,j,bi,bj)
748     anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
749     & tpobs(i,j,bi,bj)
750     anom_psMtpobs_NUM(i,j,bi,bj)=
751     & anom_psMtpobs_NUM(i,j,bi,bj)+1. _d 0
752     anom_psMslaobs_NUM(i,j,bi,bj)=
753     & anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
754     endif
755     #endif
756     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
757     if ( ersmask(i,j,bi,bj)*mean_psMersobs_MSK(i,j,bi,bj)
758     & .NE.0. ) then
759     anom_psMersobs(i,j,bi,bj)= anom_psMersobs(i,j,bi,bj)+
760     & psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)
761     & -mean_psMersobs(i,j,bi,bj)
762     anom_psMersobs_NUM(i,j,bi,bj)=
763     & anom_psMersobs_NUM(i,j,bi,bj)+1. _d 0
764     anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
765     & psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)
766     & -mean_psMersobs(i,j,bi,bj)
767     anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
768     & ersobs(i,j,bi,bj)
769     anom_psMslaobs_NUM(i,j,bi,bj)=
770     & anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
771     endif
772     #endif
773     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
774     if ( gfomask(i,j,bi,bj)*mean_psMgfoobs_MSK(i,j,bi,bj)
775     & .NE.0. ) then
776     anom_psMgfoobs(i,j,bi,bj)= anom_psMgfoobs(i,j,bi,bj)+
777     & psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)
778     & -mean_psMgfoobs(i,j,bi,bj)
779     anom_psMgfoobs_NUM(i,j,bi,bj)=
780     & anom_psMgfoobs_NUM(i,j,bi,bj)+1. _d 0
781     anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
782     & psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)
783     & -mean_psMgfoobs(i,j,bi,bj)
784     anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
785     & gfoobs(i,j,bi,bj)
786     anom_psMslaobs_NUM(i,j,bi,bj)=
787     & anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
788     endif
789     #endif
790     enddo
791     enddo
792     enddo
793     enddo
794    
795     enddo !do jrec=1,ndaysave
796    
797     do bj = jtlo,jthi
798     do bi = itlo,ithi
799     do j = jmin,jmax
800     do i = imin,imax
801     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
802     if ( anom_psMtpobs_NUM(i,j,bi,bj) .NE. 0. ) then
803     anom_psMtpobs(i,j,bi,bj) =
804     & anom_psMtpobs(i,j,bi,bj) /
805     & anom_psMtpobs_NUM(i,j,bi,bj)
806     endif
807     #endif
808     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
809     if ( anom_psMersobs_NUM(i,j,bi,bj) .NE. 0. ) then
810     anom_psMersobs(i,j,bi,bj) =
811     & anom_psMersobs(i,j,bi,bj) /
812     & anom_psMersobs_NUM(i,j,bi,bj)
813     endif
814     #endif
815     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
816     if ( anom_psMgfoobs_NUM(i,j,bi,bj) .NE. 0. ) then
817     anom_psMgfoobs(i,j,bi,bj) =
818     & anom_psMgfoobs(i,j,bi,bj) /
819     & anom_psMgfoobs_NUM(i,j,bi,bj)
820     endif
821     #endif
822     if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
823 gforget 1.3 & ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
824     & ( abs(YC(i,j,bi,bj)) .LE. 66. ) ) then
825 gforget 1.1 anom_psMslaobs(i,j,bi,bj) =
826     & anom_psMslaobs(i,j,bi,bj) /
827     & anom_psMslaobs_NUM(i,j,bi,bj)
828     anom_slaobs(i,j,bi,bj) =
829     & anom_slaobs(i,j,bi,bj) /
830     & anom_psMslaobs_NUM(i,j,bi,bj)
831     else
832     anom_psMslaobs(i,j,bi,bj) = 0. _d 0
833     anom_slaobs(i,j,bi,bj) = 0. _d 0
834     endif
835     enddo
836     enddo
837     enddo
838     enddo
839    
840     c PART 4.2: smooth anom_psMslaobs in space
841     c ----------------------------------------
842    
843     #ifdef ALLOW_GENCOST_SSHV4_OUTPUT
844     write(fname4test(1:80),'(1a)') 'sladiff_raw'
845     call mdswritefield(fname4test,32,.false.,'RL',
846     & 1,anom_psMslaobs,irec,1,mythid)
847    
848     write(fname4test(1:80),'(1a)') 'slaobs_raw'
849     call mdswritefield(fname4test,32,.false.,'RL',
850     & 1,anom_slaobs,irec,1,mythid)
851     #endif
852    
853     call smooth_hetero2d(anom_psMslaobs,maskc,
854     & gencost_scalefile(igen_lsc),300,mythid)
855    
856     #ifdef ALLOW_GENCOST_SSHV4_OUTPUT
857     call smooth_hetero2d(anom_slaobs,maskc,
858     & gencost_scalefile(igen_lsc),300,mythid)
859    
860     write(fname4test(1:80),'(1a)') 'sladiff_smooth'
861     call mdswritefield(fname4test,32,.false.,'RL',
862     & 1,anom_psMslaobs,irec,1,mythid)
863    
864     write(fname4test(1:80),'(1a)') 'slaobs_smooth'
865     call mdswritefield(fname4test,32,.false.,'RL',
866     & 1,anom_slaobs,irec,1,mythid)
867     #endif
868    
869     c PART 4.3: compute cost function term
870     c ------------------------------------
871    
872     do bj = jtlo,jthi
873     do bi = itlo,ithi
874     do j = jmin,jmax
875     do i = imin,imax
876     # if (defined (ALLOW_SSH_GFOANOM_COST_CONTRIBUTION) || \
877     defined (ALLOW_SSH_TPANOM_COST_CONTRIBUTION) || \
878     defined (ALLOW_SSH_ERSANOM_COST_CONTRIBUTION))
879     junk = anom_psMslaobs(i,j,bi,bj)
880     objf_gencost(igen_lsc,bi,bj) = objf_gencost(igen_lsc,bi,bj)
881     & + junk*junk*gencost_weight(i,j,bi,bj,igen_lsc)
882     & *maskc(i,j,1,bi,bj)/ndaysaveRL
883     if ( (gencost_weight(i,j,bi,bj,igen_lsc).GT.0.).AND.
884     & (anom_psMslaobs_NUM(i,j,bi,bj).GT.0.).AND.
885     & (maskc(i,j,1,bi,bj) .ne. 0.) )
886     & num_gencost(igen_lsc,bi,bj) =
887     & num_gencost(igen_lsc,bi,bj) + 1. _d 0 /ndaysaveRL
888     #endif
889     enddo
890     enddo
891     enddo
892     enddo
893    
894     enddo
895    
896    
897     cgf =======================================================
898     cgf PART 5: compute raw SLA cost term
899     cgf =======================================================
900    
901    
902     do irec = 1, ndaysrec
903    
904     call active_read_xy( fname, psbar, irec, doglobalread,
905     & ladinit, optimcycle, mythid,
906     & xx_psbar_mean_dummy )
907    
908 gforget 1.3 #ifndef ALLOW_PSBAR_MEAN
909     CALL REMOVE_MEAN_RL( 1, psbar, maskInC, maskInC, rA, drF,
910     & 'psbar', myTime, myThid )
911     #endif
912    
913 gforget 1.1 #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
914     call cost_readtopex( irec, mythid )
915     #endif
916     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
917     call cost_readers( irec, mythid )
918     #endif
919     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
920     call cost_readgfo( irec, mythid )
921     #endif
922    
923     do bj = jtlo,jthi
924     do bi = itlo,ithi
925     do j = jmin,jmax
926     do i = imin,imax
927     anom_psMtpobs(i,j,bi,bj) = 0. _d 0
928     anom_psMersobs(i,j,bi,bj) = 0. _d 0
929     anom_psMgfoobs(i,j,bi,bj) = 0. _d 0
930     anom_tpobs(i,j,bi,bj) = 0. _d 0
931     anom_ersobs(i,j,bi,bj) = 0. _d 0
932     anom_gfoobs(i,j,bi,bj) = 0. _d 0
933     enddo
934     enddo
935     enddo
936     enddo
937    
938     do bj = jtlo,jthi
939     do bi = itlo,ithi
940     do j = jmin,jmax
941     do i = imin,imax
942     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
943     if ( tpmask(i,j,bi,bj)*mean_psMtpobs_MSK(i,j,bi,bj).NE.0. )
944     & then
945     anom_psMtpobs(i,j,bi,bj)=
946     & psbar(i,j,bi,bj) - tpobs(i,j,bi,bj)
947     & - mean_psMtpobs(i,j,bi,bj)
948     anom_tpobs(i,j,bi,bj)=tpobs(i,j,bi,bj)
949     endif
950     #endif
951     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
952     if ( ersmask(i,j,bi,bj)*mean_psMersobs_MSK(i,j,bi,bj).NE.0. )
953     & then
954     anom_psMersobs(i,j,bi,bj)=
955     & psbar(i,j,bi,bj) - ersobs(i,j,bi,bj)
956     & - mean_psMersobs(i,j,bi,bj)
957     anom_ersobs(i,j,bi,bj)=ersobs(i,j,bi,bj)
958     endif
959     #endif
960     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
961     if ( gfomask(i,j,bi,bj)*mean_psMgfoobs_MSK(i,j,bi,bj).NE.0. )
962     & then
963     anom_psMgfoobs(i,j,bi,bj)=
964     & psbar(i,j,bi,bj) - gfoobs(i,j,bi,bj)
965     & - mean_psMgfoobs(i,j,bi,bj)
966     anom_gfoobs(i,j,bi,bj)=gfoobs(i,j,bi,bj)
967     endif
968     #endif
969     enddo
970     enddo
971     enddo
972     enddo
973    
974     #ifdef ALLOW_GENCOST_SSHV4_OUTPUT
975     write(fname4test(1:80),'(1a)') 'sladiff_tp_raw'
976     call mdswritefield(fname4test,32,.false.,'RL',
977     & 1,anom_psMtpobs,irec,1,mythid)
978     write(fname4test(1:80),'(1a)') 'sladiff_ers_raw'
979     call mdswritefield(fname4test,32,.false.,'RL',
980     & 1,anom_psMersobs,irec,1,mythid)
981     write(fname4test(1:80),'(1a)') 'sladiff_gfo_raw'
982     call mdswritefield(fname4test,32,.false.,'RL',
983     & 1,anom_psMgfoobs,irec,1,mythid)
984    
985     write(fname4test(1:80),'(1a)') 'slaobs_tp_raw'
986     call mdswritefield(fname4test,32,.false.,'RL',
987     & 1,anom_tpobs,irec,1,mythid)
988     write(fname4test(1:80),'(1a)') 'slaobs_ers_raw'
989     call mdswritefield(fname4test,32,.false.,'RL',
990     & 1,anom_ersobs,irec,1,mythid)
991     write(fname4test(1:80),'(1a)') 'slaobs_gfo_raw'
992     call mdswritefield(fname4test,32,.false.,'RL',
993     & 1,anom_gfoobs,irec,1,mythid)
994     #endif
995    
996     do bj = jtlo,jthi
997     do bi = itlo,ithi
998     do j = jmin,jmax
999     do i = imin,imax
1000     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
1001     c-- The array psobs contains SSH anomalies.
1002     junkweight = mean_psMtpobs_MSK(i,j,bi,bj)
1003     & *gencost_weight(i,j,bi,bj,igen_tp)
1004     & *tpmask(i,j,bi,bj)
1005     & *cosphi(i,j,bi,bj)
1006     junk = anom_psMtpobs(i,j,bi,bj)
1007     objf_gencost(igen_tp,bi,bj) =
1008     & objf_gencost(igen_tp,bi,bj)+junk*junk*junkweight
1009     if ( junkweight .ne. 0. )
1010     & num_gencost(igen_tp,bi,bj) =
1011     & num_gencost(igen_tp,bi,bj) + 1. _d 0
1012     #endif
1013     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
1014     c-- The array ersobs contains SSH anomalies.
1015     junkweight = mean_psMersobs_MSK(i,j,bi,bj)
1016     & *gencost_weight(i,j,bi,bj,igen_ers)
1017     & *ersmask(i,j,bi,bj)
1018     & *cosphi(i,j,bi,bj)
1019     junk = anom_psMersobs(i,j,bi,bj)
1020     objf_gencost(igen_ers,bi,bj) =
1021     & objf_gencost(igen_ers,bi,bj)+junk*junk*junkweight
1022     if ( junkweight .ne. 0. )
1023     & num_gencost(igen_ers,bi,bj) =
1024     & num_gencost(igen_ers,bi,bj) + 1. _d 0
1025     #endif
1026     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
1027     c-- The array gfoobs contains SSH anomalies.
1028     junkweight = mean_psMgfoobs_MSK(i,j,bi,bj)
1029     & *gencost_weight(i,j,bi,bj,igen_gfo)
1030     & *gfomask(i,j,bi,bj)
1031     & *cosphi(i,j,bi,bj)
1032     junk = anom_psMgfoobs(i,j,bi,bj)
1033     objf_gencost(igen_gfo,bi,bj) =
1034     & objf_gencost(igen_gfo,bi,bj)+junk*junk*junkweight
1035     if ( junkweight .ne. 0. )
1036     & num_gencost(igen_gfo,bi,bj) =
1037     & num_gencost(igen_gfo,bi,bj) + 1. _d 0
1038     #endif
1039     enddo
1040     enddo
1041     enddo
1042     enddo
1043    
1044     enddo
1045    
1046    
1047     #endif /* ifdef ALLOW_GENCOST_FREEFORM */
1048     #endif /* ifdef ALLOW_GENCOST_CONTRIBUTION */
1049     #endif /* ifdef ALLOW_SMOOTH */
1050     #endif /* ifdef ALLOW_SSH_COST_CONTRIBUTION */
1051    
1052     end
  ViewVC Help
Powered by ViewVC 1.1.22