/[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.12 - (hide annotations) (download)
Tue Nov 19 19:56:51 2013 UTC (10 years, 6 months ago) by gforget
Branch: MAIN
Changes since 1.11: +7 -6 lines
- fix bug, reported by O. Wang, in MDT re-referencing, which was always off.

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

  ViewVC Help
Powered by ViewVC 1.1.22