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

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

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

revision 1.2 by gforget, Mon Oct 25 20:42:31 2010 UTC revision 1.21 by gforget, Mon Jun 9 17:52:28 2014 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "COST_CPPOPTIONS.h"  #include "ECCO_OPTIONS.h"
5    
6    
7        subroutine cost_gencost_sshv4(        subroutine cost_gencost_sshv4(
# Line 44  c     == global variables == Line 44  c     == global variables ==
44  #include "GRID.h"  #include "GRID.h"
45    
46  #include "ecco_cost.h"  #include "ecco_cost.h"
47    #include "CTRL_SIZE.h"
48  #include "ctrl.h"  #include "ctrl.h"
49  #include "ctrl_dummy.h"  #include "ctrl_dummy.h"
50  #include "optim.h"  #include "optim.h"
# Line 51  c     == global variables == Line 52  c     == global variables ==
52  #ifdef ALLOW_PROFILES  #ifdef ALLOW_PROFILES
53  #include "profiles.h"  #include "profiles.h"
54  #endif  #endif
55    #include "cal.h"
56    #ifdef ALLOW_SMOOTH
57    #include "SMOOTH.h"
58    #endif
59    
60  c     == routine arguments ==  c     == routine arguments ==
61    
# Line 59  c     == routine arguments == Line 64  c     == routine arguments ==
64        integer mythid        integer mythid
65    
66  #ifdef ALLOW_SSH_COST_CONTRIBUTION  #ifdef ALLOW_SSH_COST_CONTRIBUTION
 #ifdef ALLOW_SMOOTH  
67  #ifdef ALLOW_GENCOST_CONTRIBUTION  #ifdef ALLOW_GENCOST_CONTRIBUTION
 #ifdef ALLOW_GENCOST_FREEFORM  
68  c     == local variables ==  c     == local variables ==
69    
70        integer bi,bj        integer bi,bj
# Line 79  c     == local variables == Line 82  c     == local variables ==
82    
83  c mapping to gencost  c mapping to gencost
84        integer igen_mdt, igen_lsc        integer igen_mdt, igen_lsc
85        integer igen_tp, igen_ers, igen_gfo              integer igen_tp, igen_ers, igen_gfo
86    
87        _RL offset,fac        _RL offset
88        _RL offset_sum        _RL offset_sum
89        _RL psmean    ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )        _RL slaoffset
90          _RL slaoffset_sum
91    
92        _RL diagnosfld ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )        _RL diagnosfld ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
93    
94  c for PART 1: re-reference MDT (tpmean) to the inferred SLA reference field  c for PART 1: re-reference MDT (mdt) to the inferred SLA reference field
95        _RL mean_slaobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL psmean    ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
96          _RL mean_slaobs_mdt(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
97        _RL mean_slaobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL mean_slaobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
98    
99  c for PART 2: compute time mean differences over the model period  c for PART 2: compute time mean differences over the model period
100        _RL mean_slaobs2(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL mean_slaobs_model(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
101        _RL mean_psMssh_all(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL mean_psMssh_all(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
102        _RL mean_psMssh_all_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL mean_psMssh_all_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
103        _RL mean_psMssh_all_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL mean_psMssh_all_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
104    
105          _RL mean_psMslaobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
106          _RL mean_psMslaobs_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
107    
108        _RL mean_psMtpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL mean_psMtpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
109        _RL mean_psMtpobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL mean_psMtpobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
       _RL mean_psMtpobs_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)  
110    
111        _RL mean_psMersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL mean_psMersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
112        _RL mean_psMersobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL mean_psMersobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
       _RL mean_psMersobs_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)  
113    
114        _RL mean_psMgfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL mean_psMgfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
115        _RL mean_psMgfoobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL mean_psMgfoobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
       _RL mean_psMgfoobs_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)  
116    
117  c for PART 4/5: compute smooth/raw anomalies  c for PART 4/5: compute smooth/raw anomalies
118        _RL anom_psMslaobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL anom_psMslaobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
       _RL anom_slaobs (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)  
119        _RL anom_psMslaobs_NUM (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL anom_psMslaobs_NUM (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
120          _RL anom_psMslaobs_MSK (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
121          _RL anom_psMslaobs_SUB (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
122          _RL anom_slaobs (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
123          _RL anom_slaobs_SUB  (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
124    
125        _RL anom_psMtpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL anom_psMtpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
       _RL anom_psMtpobs_NUM (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)  
126        _RL anom_tpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL anom_tpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
127    
128        _RL anom_psMersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL anom_psMersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
       _RL anom_psMersobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)  
129        _RL anom_ersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL anom_ersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
130    
131        _RL anom_psMgfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL anom_psMgfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
       _RL anom_psMgfoobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)  
132        _RL anom_gfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL anom_gfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
133    
134        integer tpmean_y0,tpmean_y1,year,day        integer mdt_y0,mdt_y1,year,day
135        integer num_var        integer num_var, num0
136    
137        _RL junk,junkweight        _RL junk,junkweight
138    
# Line 137  c for PART 4/5: compute smooth/raw anoma Line 143  c for PART 4/5: compute smooth/raw anoma
143        character*(80) fname4test        character*(80) fname4test
144        character*(MAX_LEN_MBUF) msgbuf        character*(MAX_LEN_MBUF) msgbuf
145    
146          LOGICAL doReference
147    
148  c     == external functions ==  c     == external functions ==
149    
150        integer  ilnblnk        integer  ilnblnk
# Line 179  c--   First, read tiled data. Line 187  c--   First, read tiled data.
187    
188  cgf =======================================================  cgf =======================================================
189  cgf PART 1:  cgf PART 1:
190  cgf        x Get the MDT (tpmean) ... compute the sample mean  cgf        x Get the MDT (mdt) ... compute the sample mean
191  cgf        (mean_slaobs) of the SLA data (i.e. RADS for tp, ers, and gfo  cgf        (mean_slaobs_mdt) of the SLA data (i.e. RADS for tp, ers, and gfo
192  cgf        together) over the time interval of the MDT ... subtract  cgf        together) over the time interval of the MDT ... subtract
193  cgf        mean_slaobs from tpmean.  cgf        mean_slaobs_mdt from mdt.
194  cgf        x At this point, tpmean is the inferred SLA reference field.  cgf        x At this point, mdt is the inferred SLA reference field.
195  cgf        x From there on, tpmean+sla will be directly comparable to  cgf        x From there on, mdt+sla will be directly comparable to
196  cgf        the model SSH (psbar).  cgf        the model SSH (psbar).
197  cgf =======================================================  cgf =======================================================
198    
199  c--   Read mean field and mask  c--   Read mean field and mask
200        call cost_ReadTopexMean( mythid )        call cost_ReadTopexMean( mythid )
201    
202  c--   Compute mean_slaobs: sample mean SLA over the time period of tpmean.  c--   Compute mean_slaobs_mdt: sample mean SLA over the time period of mdt.
203    
204  c pavlis and ecco/rio  c pavlis and ecco/rio
205        tpmean_y0=1993  c      mdt_y0=1993
206        tpmean_y1=2004  c      mdt_y1=2004
207  c maximenko  c maximenko
208  c      tpmean_y0=1992  c      mdt_y0=1992
209  c      tpmean_y1=2002  c      mdt_y1=2002
210  c rio  c rio
211  c      tpmean_y0=1993  c      mdt_y0=1993
212  c      tpmean_y1=1999  c      mdt_y1=1999
213    c generalized:
214           mdt_y0=mdtstartdate(1)/10000
215           mdt_y1=mdtenddate(1)/10000
216    
217           write(msgbuf,'(a,i8,a,i8)')
218         &  'mdt[start,end]date(1): ',mdtstartdate(1),
219         &  ',',mdtenddate(1)
220            call print_message( msgbuf, standardmessageunit,
221         &                      SQUEEZE_RIGHT , mythid)
222    
223            write(msgbuf,'(a,i4,a,i4)')
224         &  '  TP MDT yrrange:          ',
225         &                        mdt_y0,',',
226         &                        mdt_y1
227            call print_message( msgbuf, standardmessageunit,
228         &                      SQUEEZE_RIGHT , mythid)
229    
230    c minimum number of observations to consider re-referenced MDT
231          num0=100
232    
233          doReference=.FALSE.
234          if ((modelstartdate(1).GT.1992*10000).AND.
235         &    (modelstartdate(1).LT.2011*10000).AND.
236         &    (ndaysrec.GE.365))  doReference=.TRUE.
237    
238            write(msgbuf,'(a,l)') ' sshv4:re-reference MDT=',doReference
239            call print_message( msgbuf, standardmessageunit,
240         &                          SQUEEZE_RIGHT , mythid)
241    
242         do bj = jtlo,jthi         do bj = jtlo,jthi
243          do bi = itlo,ithi          do bi = itlo,ithi
244           do j = jmin,jmax           do j = jmin,jmax
245            do i = imin,imax            do i = imin,imax
246                mean_slaobs(i,j,bi,bj)  = 0. _d 0                mean_slaobs_mdt(i,j,bi,bj)  = 0. _d 0
247                mean_slaobs_NUM(i,j,bi,bj)  = 0. _d 0                mean_slaobs_NUM(i,j,bi,bj)  = 0. _d 0
248            enddo            enddo
249           enddo           enddo
250          enddo          enddo
251         enddo         enddo
252    
253        do year=tpmean_y0,tpmean_y1        do year=mdt_y0,mdt_y1
254         do day=1,366         do day=1,366
255  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
256        call cost_sla_read_yd( topexfile,        call cost_sla_read_yd( topexfile, topexstartdate,
257       &                tpobs, tpmask,       &                tpobs, tpmask,
258       &                year, day, mythid )       &                year, day, mythid )
259  #endif  #endif
260  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
261        call cost_sla_read_yd( ersfile,        call cost_sla_read_yd( ersfile, ersstartdate,
262       &                ersobs, ersmask,       &                ersobs, ersmask,
263       &                year, day, mythid )       &                year, day, mythid )
264  #endif  #endif
265  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
266        call cost_sla_read_yd( gfofile,        call cost_sla_read_yd( gfofile, gfostartdate,
267       &                gfoobs, gfomask,       &                gfoobs, gfomask,
268       &                year, day, mythid )       &                year, day, mythid )
269  #endif  #endif
# Line 236  c      tpmean_y1=1999 Line 272  c      tpmean_y1=1999
272           do j = jmin,jmax           do j = jmin,jmax
273            do i = imin,imax            do i = imin,imax
274  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
275        if ( tpmask(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)*        if ( tpmask(i,j,bi,bj)*mdtmask(i,j,bi,bj)*
276       &    gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then       &    gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then
277            mean_slaobs(i,j,bi,bj)= mean_slaobs(i,j,bi,bj)+            mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+
278       &  tpobs(i,j,bi,bj)       &  tpobs(i,j,bi,bj)
279            mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0            mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
280        endif        endif
281  #endif  #endif
282  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
283        if ( ersmask(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)*        if ( ersmask(i,j,bi,bj)*mdtmask(i,j,bi,bj)*
284       &    gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then       &    gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then
285            mean_slaobs(i,j,bi,bj)= mean_slaobs(i,j,bi,bj)+            mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+
286       &  ersobs(i,j,bi,bj)       &  ersobs(i,j,bi,bj)
287            mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0            mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
288        endif        endif
289  #endif  #endif
290  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
291        if ( gfomask(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)*        if ( gfomask(i,j,bi,bj)*mdtmask(i,j,bi,bj)*
292       &    gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then       &    gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then
293            mean_slaobs(i,j,bi,bj)= mean_slaobs(i,j,bi,bj)+            mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+
294       &  gfoobs(i,j,bi,bj)       &  gfoobs(i,j,bi,bj)
295            mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0            mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
296        endif        endif
# Line 265  c      tpmean_y1=1999 Line 301  c      tpmean_y1=1999
301         enddo         enddo
302    
303         enddo !do day=1,366         enddo !do day=1,366
304        enddo !do year=tpmean_y0,tpmean_y1        enddo !do year=mdt_y0,mdt_y1
305    
306        do bj = jtlo,jthi        do bj = jtlo,jthi
307          do bi = itlo,ithi          do bi = itlo,ithi
# Line 273  c      tpmean_y1=1999 Line 309  c      tpmean_y1=1999
309              do i = imin,imax              do i = imin,imax
310                 if ( ( mean_slaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.                 if ( ( mean_slaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
311       &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.       &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
312    #ifndef ALLOW_HIGHLAT_ALTIMETRY
313       &              ( abs(YC(i,j,bi,bj)) .LE. 66. ).AND.       &              ( abs(YC(i,j,bi,bj)) .LE. 66. ).AND.
314       &              ( tpmeanmask(i,j,bi,bj) .NE. 0. ) ) then  #endif
315                    mean_slaobs(i,j,bi,bj) = mean_slaobs(i,j,bi,bj) /       &              ( mdtmask(i,j,bi,bj) .NE. 0. ) ) then
316                      mean_slaobs_mdt(i,j,bi,bj) =
317         &                 mean_slaobs_mdt(i,j,bi,bj) /
318       &                 mean_slaobs_NUM(i,j,bi,bj)       &                 mean_slaobs_NUM(i,j,bi,bj)
319                 else                 else
320                    mean_slaobs(i,j,bi,bj) = 0. _d 0                    mean_slaobs_mdt(i,j,bi,bj) = 0. _d 0
321                 endif                 endif
322              enddo              enddo
323            enddo            enddo
# Line 286  c      tpmean_y1=1999 Line 325  c      tpmean_y1=1999
325        enddo        enddo
326    
327    
328  c--   smooth mean_slaobs:  c--   smooth mean_slaobs_mdt:
329    
330        write(fname4test(1:80),'(1a)') 'sla2mdt_raw'        write(fname4test(1:80),'(1a)') 'sla2mdt_raw'
331        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
332       & 1,mean_slaobs,1,1,mythid)       & 1,mean_slaobs_mdt,1,1,mythid)
333    
334        call smooth_hetero2d(mean_slaobs,maskc,  #ifdef ALLOW_SMOOTH
335       &     gencost_scalefile(igen_mdt),300,mythid)        if ( useSMOOTH )
336         &  call smooth_hetero2d(mean_slaobs_mdt,maskc,
337         &     gencost_scalefile(igen_mdt),
338         &     gencost_smooth2Ddiffnbt(igen_mdt),mythid)
339    #endif
340    
341        write(fname4test(1:80),'(1a)') 'sla2mdt_smooth'        write(fname4test(1:80),'(1a)') 'sla2mdt_smooth'
342        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
343       & 1,mean_slaobs,1,1,mythid)       & 1,mean_slaobs_mdt,1,1,mythid)
344    
345  c--   re-reference tpmean:  c--   re-reference mdt:
346        do bj = jtlo,jthi        do bj = jtlo,jthi
347          do bi = itlo,ithi          do bi = itlo,ithi
348            do j = jmin,jmax            do j = jmin,jmax
349              do i = imin,imax              do i = imin,imax
350                 if ( ( tpmeanmask(i,j,bi,bj) .NE. 0. ).AND.                 if ( ( mdtmask(i,j,bi,bj) .NE. 0. ).AND.
351       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) ) then       &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
352                    tpmean(i,j,bi,bj) = tpmean(i,j,bi,bj)       &              ( doReference ) ) then
353       &                 -mean_slaobs(i,j,bi,bj)                    mdt(i,j,bi,bj) = mdt(i,j,bi,bj)
354         &                 -mean_slaobs_mdt(i,j,bi,bj)
355                 endif                 endif
356              enddo              enddo
357            enddo            enddo
# Line 330  cgf ==================================== Line 374  cgf ====================================
374              do i = imin,imax              do i = imin,imax
375    
376                psmean(i,j,bi,bj)    = 0. _d 0                psmean(i,j,bi,bj)    = 0. _d 0
377                  mean_psMslaobs(i,j,bi,bj)  = 0. _d 0
378                mean_psMtpobs(i,j,bi,bj)  = 0. _d 0                mean_psMtpobs(i,j,bi,bj)  = 0. _d 0
379                mean_psMersobs(i,j,bi,bj) = 0. _d 0                mean_psMersobs(i,j,bi,bj) = 0. _d 0
380                mean_psMgfoobs(i,j,bi,bj) = 0. _d 0                mean_psMgfoobs(i,j,bi,bj) = 0. _d 0
381                mean_psMssh_all(i,j,bi,bj) = 0. _d 0                mean_psMssh_all(i,j,bi,bj) = 0. _d 0
382                mean_slaobs2(i,j,bi,bj)  = 0. _d 0                mean_slaobs_model(i,j,bi,bj)  = 0. _d 0
383    
384                mean_psMtpobs_NUM(i,j,bi,bj)  = 0. _d 0                mean_psMtpobs_NUM(i,j,bi,bj)  = 0. _d 0
385                mean_psMersobs_NUM(i,j,bi,bj) = 0. _d 0                mean_psMersobs_NUM(i,j,bi,bj) = 0. _d 0
386                mean_psMgfoobs_NUM(i,j,bi,bj) = 0. _d 0                mean_psMgfoobs_NUM(i,j,bi,bj) = 0. _d 0
387                mean_psMssh_all_NUM(i,j,bi,bj) = 0. _d 0                mean_psMssh_all_NUM(i,j,bi,bj) = 0. _d 0
388                  mean_psMslaobs_MSK(i,j,bi,bj)  = 0. _d 0
               mean_psMtpobs_MSK(i,j,bi,bj)  = 0. _d 0  
               mean_psMersobs_MSK(i,j,bi,bj) = 0. _d 0  
               mean_psMgfoobs_MSK(i,j,bi,bj) = 0. _d 0  
389    
390              enddo              enddo
391            enddo            enddo
# Line 359  cgf ==================================== Line 401  cgf ====================================
401       &                       ladinit, optimcycle, mythid,       &                       ladinit, optimcycle, mythid,
402       &                       xx_psbar_mean_dummy )       &                       xx_psbar_mean_dummy )
403    
404    #ifndef ALLOW_PSBAR_MEAN
405            CALL REMOVE_MEAN_RL( 1, psbar, maskInC, maskInC, rA, drF,
406         &        'psbar', myTime, myThid )
407    #endif
408    
409  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
410        call cost_sla_read( topexfile, topexstartdate, topexperiod,        call cost_sla_read( topexfile, topexstartdate, topexperiod,
411       &                topexintercept, topexslope,       &                topexintercept, topexslope,
# Line 387  cgf ==================================== Line 434  cgf ====================================
434  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
435                  if ( tpmask(i,j,bi,bj)*                  if ( tpmask(i,j,bi,bj)*
436       &             gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then       &             gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then
437                     mean_slaobs2(i,j,bi,bj)=                     mean_slaobs_model(i,j,bi,bj)=
438       &                 mean_slaobs2(i,j,bi,bj)+tpobs(i,j,bi,bj)       &                 mean_slaobs_model(i,j,bi,bj)+tpobs(i,j,bi,bj)
439                     mean_psMtpobs(i,j,bi,bj) =                     mean_psMtpobs(i,j,bi,bj) =
440       &                 mean_psMtpobs(i,j,bi,bj) +       &                 mean_psMtpobs(i,j,bi,bj) +
441       &                 psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)       &                 psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)
# Line 399  cgf ==================================== Line 446  cgf ====================================
446  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
447                  if ( ersmask(i,j,bi,bj)*                  if ( ersmask(i,j,bi,bj)*
448       &             gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then       &             gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then
449                     mean_slaobs2(i,j,bi,bj)=                     mean_slaobs_model(i,j,bi,bj)=
450       &                 mean_slaobs2(i,j,bi,bj)+ersobs(i,j,bi,bj)       &                 mean_slaobs_model(i,j,bi,bj)+ersobs(i,j,bi,bj)
451                     mean_psMersobs(i,j,bi,bj) =                     mean_psMersobs(i,j,bi,bj) =
452       &                 mean_psMersobs(i,j,bi,bj) +       &                 mean_psMersobs(i,j,bi,bj) +
453       &                 psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)       &                 psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)
# Line 411  cgf ==================================== Line 458  cgf ====================================
458  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
459                  if ( gfomask(i,j,bi,bj)*                  if ( gfomask(i,j,bi,bj)*
460       &             gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then       &             gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then
461                     mean_slaobs2(i,j,bi,bj)=                     mean_slaobs_model(i,j,bi,bj)=
462       &                 mean_slaobs2(i,j,bi,bj)+gfoobs(i,j,bi,bj)       &                 mean_slaobs_model(i,j,bi,bj)+gfoobs(i,j,bi,bj)
463                     mean_psMgfoobs(i,j,bi,bj) =                     mean_psMgfoobs(i,j,bi,bj) =
464       &                 mean_psMgfoobs(i,j,bi,bj) +       &                 mean_psMgfoobs(i,j,bi,bj) +
465       &                 psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)       &                 psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)
# Line 433  c--   END loop over records for the firs Line 480  c--   END loop over records for the firs
480              do j = jmin,jmax              do j = jmin,jmax
481                do i = imin,imax                do i = imin,imax
482  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
483                 if ( mean_psMtpobs_NUM(i,j,bi,bj) .NE. 0. ) then                 if ( ( mean_psMtpobs_NUM(i,j,bi,bj) .NE. 0. )
484    #ifndef ALLOW_HIGHLAT_ALTIMETRY
485         &              .AND.( abs(YC(i,j,bi,bj)) .LE. 66. )
486    #endif
487         &              ) then
488                    mean_psMssh_all(i,j,bi,bj) =                    mean_psMssh_all(i,j,bi,bj) =
489       &                 mean_psMssh_all(i,j,bi,bj) +       &                 mean_psMssh_all(i,j,bi,bj) +
490       &                 mean_psMtpobs(i,j,bi,bj)       &                 mean_psMtpobs(i,j,bi,bj)
# Line 443  c--   END loop over records for the firs Line 494  c--   END loop over records for the firs
494                    mean_psMtpobs(i,j,bi,bj) =                    mean_psMtpobs(i,j,bi,bj) =
495       &                 mean_psMtpobs(i,j,bi,bj) /       &                 mean_psMtpobs(i,j,bi,bj) /
496       &                 mean_psMtpobs_NUM(i,j,bi,bj)       &                 mean_psMtpobs_NUM(i,j,bi,bj)
                   mean_psMtpobs_MSK(i,j,bi,bj) = 1. _d 0  
497                 endif                 endif
498  #endif  #endif
499  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
500                 if ( mean_psMersobs_NUM(i,j,bi,bj) .NE. 0. ) then                 if ( ( mean_psMersobs_NUM(i,j,bi,bj) .NE. 0. )
501    #ifndef ALLOW_HIGHLAT_ALTIMETRY    
502         &              .AND.( abs(YC(i,j,bi,bj)) .LE. 66. )
503    #endif
504         &              ) then
505                    mean_psMssh_all(i,j,bi,bj) =                    mean_psMssh_all(i,j,bi,bj) =
506       &                 mean_psMssh_all(i,j,bi,bj) +       &                 mean_psMssh_all(i,j,bi,bj) +
507       &                 mean_psMersobs(i,j,bi,bj)       &                 mean_psMersobs(i,j,bi,bj)
# Line 457  c--   END loop over records for the firs Line 511  c--   END loop over records for the firs
511                    mean_psMersobs(i,j,bi,bj) =                    mean_psMersobs(i,j,bi,bj) =
512       &                 mean_psMersobs(i,j,bi,bj) /       &                 mean_psMersobs(i,j,bi,bj) /
513       &                 mean_psMersobs_NUM(i,j,bi,bj)       &                 mean_psMersobs_NUM(i,j,bi,bj)
                   mean_psMersobs_MSK(i,j,bi,bj) = 1. _d 0  
514                 endif                 endif
515  #endif  #endif
516  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
517                 if ( mean_psMgfoobs_NUM(i,j,bi,bj) .NE. 0. ) then                 if ( ( mean_psMgfoobs_NUM(i,j,bi,bj) .NE. 0. )
518    #ifndef ALLOW_HIGHLAT_ALTIMETRY    
519         &              .AND.( abs(YC(i,j,bi,bj)) .LE. 66. )
520    #endif
521         &              ) then
522                    mean_psMssh_all(i,j,bi,bj) =                    mean_psMssh_all(i,j,bi,bj) =
523       &                 mean_psMssh_all(i,j,bi,bj) +       &                 mean_psMssh_all(i,j,bi,bj) +
524       &                 mean_psMgfoobs(i,j,bi,bj)       &                 mean_psMgfoobs(i,j,bi,bj)
# Line 471  c--   END loop over records for the firs Line 528  c--   END loop over records for the firs
528                    mean_psMgfoobs(i,j,bi,bj) =                    mean_psMgfoobs(i,j,bi,bj) =
529       &                 mean_psMgfoobs(i,j,bi,bj) /       &                 mean_psMgfoobs(i,j,bi,bj) /
530       &                 mean_psMgfoobs_NUM(i,j,bi,bj)       &                 mean_psMgfoobs_NUM(i,j,bi,bj)
                   mean_psMgfoobs_MSK(i,j,bi,bj) = 1. _d 0  
531                 endif                 endif
532  #endif  #endif
533                 if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .NE. 0. ).AND.  
534                   if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .GT. 0 ).AND.
535         &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
536    #ifndef ALLOW_HIGHLAT_ALTIMETRY    
537         &              ( abs(YC(i,j,bi,bj)) .LE. 66. )
538    #endif
539         &              ) then
540                      mean_psMslaobs(i,j,bi,bj) =
541         &                 mean_psMssh_all(i,j,bi,bj) /
542         &                 mean_psMssh_all_NUM(i,j,bi,bj)
543                      mean_psMslaobs_MSK(i,j,bi,bj) = 1. _d 0
544                   else
545                      mean_psMslaobs(i,j,bi,bj) = 0. _d 0
546                      mean_psMslaobs_MSK(i,j,bi,bj) = 0. _d 0
547                   endif
548    
549                   if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .GT. num0 ).AND.
550       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
551    #ifndef ALLOW_HIGHLAT_ALTIMETRY    
552       &              ( abs(YC(i,j,bi,bj)) .LE. 66. ).AND.       &              ( abs(YC(i,j,bi,bj)) .LE. 66. ).AND.
553       &              ( tpmeanmask(i,j,bi,bj) .NE. 0. ) ) then  #endif
554                    mean_slaobs2(i,j,bi,bj) =       &              ( mdtmask(i,j,bi,bj) .NE. 0. ).AND.
555       &                 mean_slaobs2(i,j,bi,bj) /       &              ( doReference ) ) then
556                      mean_slaobs_model(i,j,bi,bj) =
557         &                 mean_slaobs_model(i,j,bi,bj) /
558       &                 mean_psMssh_all_NUM(i,j,bi,bj)       &                 mean_psMssh_all_NUM(i,j,bi,bj)
559                    mean_psMssh_all(i,j,bi,bj) =                    mean_psMssh_all(i,j,bi,bj) =
560       &                 mean_psMssh_all(i,j,bi,bj) /       &                 mean_psMssh_all(i,j,bi,bj) /
561       &                 mean_psMssh_all_NUM(i,j,bi,bj)-tpmean(i,j,bi,bj)       &                 mean_psMssh_all_NUM(i,j,bi,bj)-mdt(i,j,bi,bj)
562                    mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0                    mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
563                    offset=offset+mean_psMssh_all(i,j,bi,bj)*                    offset=offset+RA(i,j,bi,bj)*mean_psMssh_all(i,j,bi,bj)
564       &                 mean_psMssh_all_NUM(i,j,bi,bj)                    offset_sum=offset_sum+RA(i,j,bi,bj)
565                    offset_sum=offset_sum+mean_psMssh_all_NUM(i,j,bi,bj)                 elseif ( ( mdtmask(i,j,bi,bj) .NE. 0. ) .AND.
566         &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
567    #ifndef ALLOW_HIGHLAT_ALTIMETRY
568         &              ( abs(YC(i,j,bi,bj)) .LE. 66. ).AND.
569    #endif
570         &              ( .NOT.doReference ) ) then
571                      mean_slaobs_model(i,j,bi,bj) = 0.d0
572                      mean_psMssh_all(i,j,bi,bj) = 0. _d 0
573                      mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
574                      offset=offset+RA(i,j,bi,bj)
575         &                  *( psmean(i,j,bi,bj) -mdt(i,j,bi,bj) )
576                      offset_sum=offset_sum+RA(i,j,bi,bj)
577                 else                 else
578                    mean_slaobs2(i,j,bi,bj) = 0.d0                    mean_slaobs_model(i,j,bi,bj) = 0.d0
579                    mean_psMssh_all(i,j,bi,bj) = 0. _d 0                    mean_psMssh_all(i,j,bi,bj) = 0. _d 0
580                    mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0                    mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0
581                 endif                 endif
# Line 498  c--   END loop over records for the firs Line 584  c--   END loop over records for the firs
584            enddo            enddo
585          enddo          enddo
586    
587    
588  c--   Do a global summation.  c--   Do a global summation.
589        _GLOBAL_SUM_RL( offset     , mythid )        _GLOBAL_SUM_RL( offset     , mythid )
590        _GLOBAL_SUM_RL( offset_sum , mythid )        _GLOBAL_SUM_RL( offset_sum , mythid )
591    
592           write(msgbuf,'(a,2d22.15)') ' sshv4:offset=',offset,offset_sum
593           call print_message( msgbuf, standardmessageunit,
594         &                          SQUEEZE_RIGHT , mythid)
595    
596        if (offset_sum .eq. 0.0) then        if (offset_sum .eq. 0.0) then
597          _BEGIN_MASTER( mythid )          _BEGIN_MASTER( mythid )
598          write(msgbuf,'(a)') ' cost_ssh: offset_sum = zero!'          write(msgbuf,'(a)') ' sshv4: offset_sum = zero!'
599          call print_message( msgbuf, standardmessageunit,          call print_message( msgbuf, standardmessageunit,
600       &                          SQUEEZE_RIGHT , mythid)       &                          SQUEEZE_RIGHT , mythid)
601          _END_MASTER( mythid )          _END_MASTER( mythid )
# Line 512  c        stop   '  ... stopped in cost_s Line 603  c        stop   '  ... stopped in cost_s
603        else        else
604          _BEGIN_MASTER( mythid )          _BEGIN_MASTER( mythid )
605          write(msgbuf,'(a,d22.15)')          write(msgbuf,'(a,d22.15)')
606       &          ' cost_ssh: offset_sum = ',offset_sum       &          ' sshv4: offset_sum = ',offset_sum
607          call print_message( msgbuf, standardmessageunit,          call print_message( msgbuf, standardmessageunit,
608       &                          SQUEEZE_RIGHT , mythid)       &                          SQUEEZE_RIGHT , mythid)
609          _END_MASTER( mythid )          _END_MASTER( mythid )
# Line 527  c--   subtract offset from mean_psMssh_a Line 618  c--   subtract offset from mean_psMssh_a
618            do j = jmin,jmax            do j = jmin,jmax
619              do i = imin,imax              do i = imin,imax
620    
621                 if ( (mean_psMssh_all_MSK(i,j,bi,bj) .NE. 0.) .AND.                 if ( (mean_psMssh_all_NUM(i,j,bi,bj) .GT. num0) .AND.
622       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.       &              ( mdtmask(i,j,bi,bj) .NE. 0. ) .AND.
623       &              ( abs(YC(i,j,bi,bj)) .LE. 66. ) ) then       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
624    #ifndef ALLOW_HIGHLAT_ALTIMETRY    
625         &              ( abs(YC(i,j,bi,bj)) .LE. 66. ) .AND.
626    #endif
627         &              ( doReference ) ) then
628  c use the re-referencing approach  c use the re-referencing approach
629                    mean_psMssh_all(i,j,bi,bj) =                     mean_psMssh_all(i,j,bi,bj) =
630       &                 mean_psMssh_all(i,j,bi,bj) - offset       &                 mean_psMssh_all(i,j,bi,bj) - offset
631                 elseif ( ( tpmeanmask(i,j,bi,bj) .NE. 0. ) .AND.                     mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
632    
633                   elseif ( ( mdtmask(i,j,bi,bj) .NE. 0. ) .AND.
634       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
635       &              ( abs(YC(i,j,bi,bj)) .GT. 66. ) ) then  #ifndef ALLOW_HIGHLAT_ALTIMETRY
636         &              ( abs(YC(i,j,bi,bj)) .LE. 66. ) .AND.
637    #endif
638         &              ( .NOT.doReference ) ) then    
639  c use the simpler approach  c use the simpler approach
640                    mean_psMssh_all(i,j,bi,bj) =                     mean_psMssh_all(i,j,bi,bj) =
641       &             psmean(i,j,bi,bj) -tpmean(i,j,bi,bj) - offset       &             psmean(i,j,bi,bj) -mdt(i,j,bi,bj) - offset
642                       mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
643    
644                 else                 else
645                    mean_psMssh_all(i,j,bi,bj) = 0. _d 0                     mean_psMssh_all(i,j,bi,bj) = 0. _d 0
646                       mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0
647                 endif                 endif
648    
649                 if ( maskc(i,j,1,bi,bj) .NE. 0. )                 if ( maskc(i,j,1,bi,bj) .NE. 0. )
650       &            psmean(i,j,bi,bj)=psmean(i,j,bi,bj)-offset       &            psmean(i,j,bi,bj)=psmean(i,j,bi,bj)-offset
651    
652              enddo              enddo
653            enddo            enddo
654          enddo          enddo
# Line 555  c--    smooth mean_psMssh_all Line 659  c--    smooth mean_psMssh_all
659        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
660       &    1,mean_psMssh_all,1,1,mythid)       &    1,mean_psMssh_all,1,1,mythid)
661    
662        call smooth_hetero2d(mean_psMssh_all,maskc,  #ifdef ALLOW_SMOOTH
663       &     gencost_scalefile(igen_mdt),300,mythid)        if ( useSMOOTH )
664         &  call smooth_hetero2d(mean_psMssh_all,maskc,
665         &     gencost_scalefile(igen_mdt),
666         &     gencost_smooth2Ddiffnbt(igen_mdt),mythid)
667    #endif
668    
669        write(fname4test(1:80),'(1a)') 'mdtdiff_smooth'        write(fname4test(1:80),'(1a)') 'mdtdiff_smooth'
670        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
# Line 564  c--    smooth mean_psMssh_all Line 672  c--    smooth mean_psMssh_all
672    
673        write(fname4test(1:80),'(1a)') 'sla2model_raw'        write(fname4test(1:80),'(1a)') 'sla2model_raw'
674        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
675       &    1,mean_slaobs2,1,1,mythid)       &    1,mean_slaobs_model,1,1,mythid)
676    
677        call smooth_hetero2d(mean_slaobs2,maskc,  #ifdef ALLOW_SMOOTH
678       &     gencost_scalefile(igen_mdt),300,mythid)        if ( useSMOOTH )
679         &  call smooth_hetero2d(mean_slaobs_model,maskc,
680         &     gencost_scalefile(igen_mdt),
681         &     gencost_smooth2Ddiffnbt(igen_mdt),mythid)
682    #endif
683    
684        write(fname4test(1:80),'(1a)') 'sla2model_smooth'        write(fname4test(1:80),'(1a)') 'sla2model_smooth'
685        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
686       &    1,mean_slaobs2,1,1,mythid)       &    1,mean_slaobs_model,1,1,mythid)
687    
688  cgf at this point:  cgf at this point:
689  cgf     1) mean_psMssh_all is the sample mean <psbar-slaobs-tpmean-offset>,  cgf     1) mean_psMssh_all is the sample mean <psbar-slaobs-mdt-offset>,
690  cgf             to which a smoothing filter has been applied.  cgf             to which a smoothing filter has been applied.
691  cgf     2) mean_psMtpobs is the (unsmoothed) sample mean <psbar-tpobs>.  cgf     2) mean_psMtpobs is the (unsmoothed) sample mean <psbar-tpobs>.
692  cgf             And similarly for ers and gfo, each treated separately.  cgf             And similarly for ers and gfo, each treated separately.
693    
694  #ifdef ALLOW_PROFILES  #ifdef ALLOW_PROFILES
695          if ( usePROFILES ) then
696        do bj = jtlo,jthi        do bj = jtlo,jthi
697          do bi = itlo,ithi          do bi = itlo,ithi
698            do j = 1,sny            do j = 1,sny
699              do i = 1,snx              do i = 1,snx
700    cgf this has not been tested recently
701                prof_etan_mean(i,j,bi,bj)=psmean(i,j,bi,bj)                prof_etan_mean(i,j,bi,bj)=psmean(i,j,bi,bj)
702              enddo              enddo
703            enddo            enddo
704          enddo          enddo
705        enddo        enddo
706        _EXCH_XY_RL( prof_etan_mean, mythid )        _EXCH_XY_RL( prof_etan_mean, mythid )
707          endif
708  #endif  #endif
709    
710    
# Line 604  cgf ==================================== Line 719  cgf ====================================
719          do bi = itlo,ithi          do bi = itlo,ithi
720            do j = jmin,jmax            do j = jmin,jmax
721              do i = imin,imax              do i = imin,imax
722         junk = mean_psMssh_all(i,j,bi,bj)         if (mean_psMssh_all_MSK(i,j,bi,bj).NE.0. _d 0) then
723         junkweight = gencost_weight(i,j,bi,bj,igen_mdt)           junk = mean_psMssh_all(i,j,bi,bj)
724       &      *tpmeanmask(i,j,bi,bj)           junkweight = gencost_weight(i,j,bi,bj,igen_mdt)
725         objf_gencost(igen_mdt,bi,bj) = objf_gencost(igen_mdt,bi,bj)       &              * mdtmask(i,j,bi,bj)
726           else
727             junk = 0. _d 0
728             junkweight = 0. _d 0
729           endif
730           objf_gencost(bi,bj,igen_mdt) = objf_gencost(bi,bj,igen_mdt)
731       &      + junk*junk*junkweight       &      + junk*junk*junkweight
732         if ( junkweight .ne. 0. ) num_gencost(igen_mdt,bi,bj) =         if ( junkweight .ne. 0. ) num_gencost(bi,bj,igen_mdt) =
733       &      num_gencost(igen_mdt,bi,bj) + 1. _d 0       &      num_gencost(bi,bj,igen_mdt) + 1. _d 0
734         diagnosfld(i,j,bi,bj) = junk*junk*junkweight         diagnosfld(i,j,bi,bj) = junk*junk*junkweight
735              enddo              enddo
736            enddo            enddo
737          enddo          enddo
738        enddo        enddo
739    
740        CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean', ' ', diagnosfld,        CALL WRITE_FLD_XY_RL( 'DiagnosMDT', ' ', diagnosfld,
741       &                           optimcycle, mythid )       &                           optimcycle, mythid )
742    
743  #endif /* ALLOW_SSH_MEAN_COST_CONTRIBUTION */  #endif /* ALLOW_SSH_MEAN_COST_CONTRIBUTION */
# Line 632  cgf ==================================== Line 752  cgf ====================================
752        ndaysave=35        ndaysave=35
753        ndaysaveRL=ndaysave        ndaysaveRL=ndaysave
754    
755        do irec = 1, ndaysrec-ndaysave+1, 7        do irec = 1, ndaysrec-ndaysave+1
756    
757         do bj = jtlo,jthi         do bj = jtlo,jthi
758          do bi = itlo,ithi          do bi = itlo,ithi
759           do j = jmin,jmax           do j = jmin,jmax
760            do i = imin,imax            do i = imin,imax
761                anom_psMslaobs(i,j,bi,bj)  = 0. _d 0                anom_psMslaobs(i,j,bi,bj)  = 0. _d 0
762                  anom_psMslaobs_MSK(i,j,bi,bj)  = 0. _d 0
763                  anom_psMslaobs_SUB(i,j,bi,bj)  = 0. _d 0
764                anom_slaobs(i,j,bi,bj)  = 0. _d 0                anom_slaobs(i,j,bi,bj)  = 0. _d 0
765                anom_psMtpobs(i,j,bi,bj)  = 0. _d 0                anom_slaobs_SUB(i,j,bi,bj)  = 0. _d 0
               anom_psMersobs(i,j,bi,bj) = 0. _d 0  
               anom_psMgfoobs(i,j,bi,bj) = 0. _d 0  
766                anom_psMslaobs_NUM(i,j,bi,bj)  = 0. _d 0                anom_psMslaobs_NUM(i,j,bi,bj)  = 0. _d 0
               anom_psMtpobs_NUM(i,j,bi,bj)  = 0. _d 0  
               anom_psMersobs_NUM(i,j,bi,bj) = 0. _d 0  
               anom_psMgfoobs_NUM(i,j,bi,bj) = 0. _d 0  
767            enddo            enddo
768           enddo           enddo
769          enddo          enddo
# Line 663  c -------------------------------------- Line 780  c --------------------------------------
780       &                       ladinit, optimcycle, mythid,       &                       ladinit, optimcycle, mythid,
781       &                       xx_psbar_mean_dummy )       &                       xx_psbar_mean_dummy )
782    
783    #ifndef ALLOW_PSBAR_MEAN
784            CALL REMOVE_MEAN_RL( 1, psbar, maskInC, maskInC, rA, drF,
785         &        'psbar', myTime, myThid )
786    #endif
787    
788  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
789        call cost_sla_read( topexfile, topexstartdate, topexperiod,        call cost_sla_read( topexfile, topexstartdate, topexperiod,
790       &                topexintercept, topexslope,       &                topexintercept, topexslope,
# Line 687  c -------------------------------------- Line 809  c --------------------------------------
809           do j = jmin,jmax           do j = jmin,jmax
810            do i = imin,imax            do i = imin,imax
811  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
812        if ( tpmask(i,j,bi,bj)*mean_psMtpobs_MSK(i,j,bi,bj)        if ( tpmask(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj)
813       &  .NE.0. ) then       &  .NE.0. ) then
            anom_psMtpobs(i,j,bi,bj)= anom_psMtpobs(i,j,bi,bj)+  
      &        psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)  
      &        -mean_psMtpobs(i,j,bi,bj)  
814             anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+             anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
815       &        psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)       &        psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)
816       &        -mean_psMtpobs(i,j,bi,bj)       &        -mean_psMslaobs(i,j,bi,bj)
817             anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+             anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
818       &        tpobs(i,j,bi,bj)       &        tpobs(i,j,bi,bj)
            anom_psMtpobs_NUM(i,j,bi,bj)=  
      &        anom_psMtpobs_NUM(i,j,bi,bj)+1. _d 0  
819             anom_psMslaobs_NUM(i,j,bi,bj)=             anom_psMslaobs_NUM(i,j,bi,bj)=
820       &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0       &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
821               if ( (2*jrec.EQ.ndaysave).OR.(2*jrec+1.EQ.ndaysave) )
822         &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0
823        endif        endif
824  #endif  #endif
825  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
826        if ( ersmask(i,j,bi,bj)*mean_psMersobs_MSK(i,j,bi,bj)        if ( ersmask(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj)
827       &  .NE.0. ) then       &  .NE.0. ) then
            anom_psMersobs(i,j,bi,bj)= anom_psMersobs(i,j,bi,bj)+  
      &        psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)  
      &        -mean_psMersobs(i,j,bi,bj)  
            anom_psMersobs_NUM(i,j,bi,bj)=  
      &        anom_psMersobs_NUM(i,j,bi,bj)+1. _d 0  
828             anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+             anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
829       &        psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)       &        psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)
830       &        -mean_psMersobs(i,j,bi,bj)       &        -mean_psMslaobs(i,j,bi,bj)
831             anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+             anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
832       &        ersobs(i,j,bi,bj)       &        ersobs(i,j,bi,bj)
833             anom_psMslaobs_NUM(i,j,bi,bj)=             anom_psMslaobs_NUM(i,j,bi,bj)=
834       &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0       &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
835               if ( (2*jrec.EQ.ndaysave).OR.(2*jrec+1.EQ.ndaysave) )
836         &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0
837        endif        endif
838  #endif  #endif
839  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
840        if ( gfomask(i,j,bi,bj)*mean_psMgfoobs_MSK(i,j,bi,bj)        if ( gfomask(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj)
841       &  .NE.0. ) then       &  .NE.0. ) then
            anom_psMgfoobs(i,j,bi,bj)= anom_psMgfoobs(i,j,bi,bj)+  
      &        psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)  
      &        -mean_psMgfoobs(i,j,bi,bj)  
            anom_psMgfoobs_NUM(i,j,bi,bj)=  
      &        anom_psMgfoobs_NUM(i,j,bi,bj)+1. _d 0  
842             anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+             anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
843       &        psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)       &        psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)
844       &        -mean_psMgfoobs(i,j,bi,bj)       &        -mean_psMslaobs(i,j,bi,bj)
845             anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+             anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
846       &        gfoobs(i,j,bi,bj)       &        gfoobs(i,j,bi,bj)
847             anom_psMslaobs_NUM(i,j,bi,bj)=             anom_psMslaobs_NUM(i,j,bi,bj)=
848       &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0       &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
849               if ( (2*jrec.EQ.ndaysave).OR.(2*jrec+1.EQ.ndaysave) )
850         &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0
851        endif        endif
852  #endif  #endif
853            enddo            enddo
# Line 748  c -------------------------------------- Line 861  c --------------------------------------
861            do bi = itlo,ithi            do bi = itlo,ithi
862              do j = jmin,jmax              do j = jmin,jmax
863                do i = imin,imax                do i = imin,imax
 #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  
                if ( anom_psMtpobs_NUM(i,j,bi,bj) .NE. 0. ) then  
                   anom_psMtpobs(i,j,bi,bj) =  
      &                 anom_psMtpobs(i,j,bi,bj) /  
      &                 anom_psMtpobs_NUM(i,j,bi,bj)  
                endif  
 #endif  
 #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION  
                if ( anom_psMersobs_NUM(i,j,bi,bj) .NE. 0. ) then  
                   anom_psMersobs(i,j,bi,bj) =  
      &                 anom_psMersobs(i,j,bi,bj) /  
      &                 anom_psMersobs_NUM(i,j,bi,bj)  
                endif  
 #endif  
 #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION  
                if ( anom_psMgfoobs_NUM(i,j,bi,bj) .NE. 0. ) then  
                   anom_psMgfoobs(i,j,bi,bj) =  
      &                 anom_psMgfoobs(i,j,bi,bj) /  
      &                 anom_psMgfoobs_NUM(i,j,bi,bj)  
                endif  
 #endif  
864                 if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.                 if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
865       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) ) then       &              ( maskc(i,j,1,bi,bj) .NE. 0. )
866    #ifndef ALLOW_HIGHLAT_ALTIMETRY    
867         &              .AND.( abs(YC(i,j,bi,bj)) .LE. 66. )
868    #endif
869         &            ) then
870                    anom_psMslaobs(i,j,bi,bj) =                    anom_psMslaobs(i,j,bi,bj) =
871       &                 anom_psMslaobs(i,j,bi,bj) /       &                 anom_psMslaobs(i,j,bi,bj) /
872       &                 anom_psMslaobs_NUM(i,j,bi,bj)       &                 anom_psMslaobs_NUM(i,j,bi,bj)
# Line 786  c -------------------------------------- Line 882  c --------------------------------------
882            enddo            enddo
883          enddo          enddo
884    
885    c PART 4.11: compute time variable offset
886    c ----------------------------------------
887    
888          slaoffset     = 0. _d 0
889          slaoffset_sum = 0. _d 0
890    
891            do bj = jtlo,jthi
892              do bi = itlo,ithi
893                do j = jmin,jmax
894                  do i = imin,imax
895                   if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
896         &              ( maskc(i,j,1,bi,bj) .NE. 0. )
897    #ifndef ALLOW_HIGHLAT_ALTIMETRY    
898         &              .AND.( abs(YC(i,j,bi,bj)) .LE. 66. )
899    #endif
900         &            ) then
901                      slaoffset=slaoffset+RA(i,j,bi,bj)*
902         &                      anom_psMslaobs(i,j,bi,bj)
903                      slaoffset_sum=slaoffset_sum+RA(i,j,bi,bj)
904                   endif
905                  enddo
906                enddo
907              enddo
908            enddo
909    
910          _GLOBAL_SUM_RL( slaoffset     , mythid )
911          _GLOBAL_SUM_RL( slaoffset_sum , mythid )
912    
913           write(msgbuf,'(a,2d22.15)') ' sshv4:slaoffset=',
914         &                         slaoffset,slaoffset_sum
915           call print_message( msgbuf, standardmessageunit,
916         &                          SQUEEZE_RIGHT , mythid)
917    
918          if (slaoffset_sum .eq. 0.0) then
919            _BEGIN_MASTER( mythid )
920            write(msgbuf,'(a)') ' sshv4: slaoffset_sum = zero!'
921            call print_message( msgbuf, standardmessageunit,
922         &                          SQUEEZE_RIGHT , mythid)
923            _END_MASTER( mythid )
924    c        stop   '  ... stopped in cost_ssh.'
925          else
926            _BEGIN_MASTER( mythid )
927            write(msgbuf,'(a,d22.15)')
928         &          ' sshv4: slaoffset_sum = ',slaoffset_sum
929            call print_message( msgbuf, standardmessageunit,
930         &                          SQUEEZE_RIGHT , mythid)
931            _END_MASTER( mythid )
932          endif
933    
934    c--   Compute (average) slaoffset
935          slaoffset = slaoffset / slaoffset_sum
936    
937            do bj = jtlo,jthi
938              do bi = itlo,ithi
939                do j = jmin,jmax
940                  do i = imin,imax
941                   if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
942         &              ( maskc(i,j,1,bi,bj) .NE. 0. )
943    #ifndef ALLOW_HIGHLAT_ALTIMETRY    
944         &              .AND.( abs(YC(i,j,bi,bj)) .LE. 66. )
945    #endif
946         &            ) then
947                      anom_psMslaobs(i,j,bi,bj) =
948         &                 anom_psMslaobs(i,j,bi,bj) - slaoffset
949                      anom_slaobs(i,j,bi,bj) =
950         &                 anom_slaobs(i,j,bi,bj) - slaoffset
951                   else
952                      anom_psMslaobs(i,j,bi,bj) = 0. _d 0
953                      anom_slaobs(i,j,bi,bj) = 0. _d 0
954                   endif
955                  enddo
956                enddo
957              enddo
958            enddo
959    
960  c PART 4.2: smooth anom_psMslaobs in space  c PART 4.2: smooth anom_psMslaobs in space
961  c ----------------------------------------  c ----------------------------------------
962    
# Line 799  c -------------------------------------- Line 970  c --------------------------------------
970       & 1,anom_slaobs,irec,1,mythid)       & 1,anom_slaobs,irec,1,mythid)
971  #endif  #endif
972    
973        call smooth_hetero2d(anom_psMslaobs,maskc,  #ifdef ALLOW_SMOOTH
974       &     gencost_scalefile(igen_lsc),300,mythid)        if ( useSMOOTH )
975         &  call smooth_hetero2d(anom_psMslaobs,maskc,
976         &     gencost_scalefile(igen_lsc),
977         &     gencost_smooth2Ddiffnbt(igen_lsc),mythid)
978    #endif
979    
980  #ifdef ALLOW_GENCOST_SSHV4_OUTPUT  #ifdef ALLOW_GENCOST_SSHV4_OUTPUT
981        call smooth_hetero2d(anom_slaobs,maskc,  #ifdef ALLOW_SMOOTH
982       &     gencost_scalefile(igen_lsc),300,mythid)        if ( useSMOOTH )
983         &  call smooth_hetero2d(anom_slaobs,maskc,
984         &     gencost_scalefile(igen_lsc),
985         &     gencost_smooth2Ddiffnbt(igen_lsc),mythid)
986    #endif
987    
988        write(fname4test(1:80),'(1a)') 'sladiff_smooth'        write(fname4test(1:80),'(1a)') 'sladiff_smooth'
989        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
# Line 825  c ------------------------------------ Line 1004  c ------------------------------------
1004  # if (defined (ALLOW_SSH_GFOANOM_COST_CONTRIBUTION) || \  # if (defined (ALLOW_SSH_GFOANOM_COST_CONTRIBUTION) || \
1005        defined (ALLOW_SSH_TPANOM_COST_CONTRIBUTION) || \        defined (ALLOW_SSH_TPANOM_COST_CONTRIBUTION) || \
1006        defined (ALLOW_SSH_ERSANOM_COST_CONTRIBUTION))        defined (ALLOW_SSH_ERSANOM_COST_CONTRIBUTION))
1007            junk = anom_psMslaobs(i,j,bi,bj)            if ( (gencost_weight(i,j,bi,bj,igen_lsc).GT.0.).AND.
1008            objf_gencost(igen_lsc,bi,bj) = objf_gencost(igen_lsc,bi,bj)       &         (anom_psMslaobs_MSK(i,j,bi,bj).GT.0.) ) then
1009                junk = anom_psMslaobs(i,j,bi,bj)
1010                anom_slaobs_SUB(i,j,bi,bj) = anom_slaobs(i,j,bi,bj)
1011                anom_psMslaobs_SUB(i,j,bi,bj) = anom_psMslaobs(i,j,bi,bj)
1012              else
1013                junk = 0. _d 0
1014                anom_slaobs_SUB(i,j,bi,bj)= 0. _d 0
1015                anom_psMslaobs_SUB(i,j,bi,bj)= 0. _d 0
1016              endif
1017              objf_gencost(bi,bj,igen_lsc) = objf_gencost(bi,bj,igen_lsc)
1018       &        + junk*junk*gencost_weight(i,j,bi,bj,igen_lsc)       &        + junk*junk*gencost_weight(i,j,bi,bj,igen_lsc)
      &        *maskc(i,j,1,bi,bj)/ndaysaveRL  
1019            if ( (gencost_weight(i,j,bi,bj,igen_lsc).GT.0.).AND.            if ( (gencost_weight(i,j,bi,bj,igen_lsc).GT.0.).AND.
1020       &         (anom_psMslaobs_NUM(i,j,bi,bj).GT.0.).AND.       &         (anom_psMslaobs_MSK(i,j,bi,bj).GT.0.) )
1021       &         (maskc(i,j,1,bi,bj) .ne. 0.) )       &         num_gencost(bi,bj,igen_lsc) =
1022       &         num_gencost(igen_lsc,bi,bj) =       &         num_gencost(bi,bj,igen_lsc) + 1. _d 0
      &         num_gencost(igen_lsc,bi,bj) + 1. _d 0 /ndaysaveRL  
1023  #endif  #endif
1024             enddo             enddo
1025           enddo           enddo
1026          enddo          enddo
1027         enddo         enddo
1028    
1029        enddo  #ifdef ALLOW_GENCOST_SSHV4_OUTPUT
1030          write(fname4test(1:80),'(1a)') 'sladiff_subsample'
1031          call mdswritefield(fname4test,32,.false.,'RL',
1032         & 1,anom_psMslaobs_SUB,irec,1,mythid)
1033    
1034          write(fname4test(1:80),'(1a)') 'slaobs_subsample'
1035          call mdswritefield(fname4test,32,.false.,'RL',
1036         & 1,anom_slaobs_SUB,irec,1,mythid)
1037    #endif
1038    
1039    
1040  cgf =======================================================  cgf =======================================================
1041  cgf PART 5: compute raw SLA cost term  cgf PART 5: compute raw SLA cost term
1042  cgf =======================================================  cgf =======================================================
1043    
1044    c time associated with this ndaysrec mean
1045            krec = irec + (ndaysave/2)
1046    
1047        do irec = 1, ndaysrec          call active_read_xy( fname, psbar, krec, doglobalread,
   
         call active_read_xy( fname, psbar, irec, doglobalread,  
1048       &                       ladinit, optimcycle, mythid,       &                       ladinit, optimcycle, mythid,
1049       &                       xx_psbar_mean_dummy )       &                       xx_psbar_mean_dummy )
1050    
1051    #ifndef ALLOW_PSBAR_MEAN
1052            CALL REMOVE_MEAN_RL( 1, psbar, maskInC, maskInC, rA, drF,
1053         &        'psbar', myTime, myThid )
1054    #endif
1055    
1056  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
1057          call cost_readtopex( irec, mythid )          call cost_readtopex( krec, mythid )
1058  #endif  #endif
1059  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
1060          call cost_readers( irec, mythid )          call cost_readers( krec, mythid )
1061  #endif  #endif
1062  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
1063          call cost_readgfo( irec, mythid )          call cost_readgfo( krec, mythid )
1064  #endif  #endif
1065    
1066         do bj = jtlo,jthi         do bj = jtlo,jthi
# Line 884  cgf ==================================== Line 1083  cgf ====================================
1083           do j = jmin,jmax           do j = jmin,jmax
1084            do i = imin,imax            do i = imin,imax
1085  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
1086        if ( tpmask(i,j,bi,bj)*mean_psMtpobs_MSK(i,j,bi,bj).NE.0. )        if ( tpmask(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj).NE.0. )
1087       & then       & then
1088           anom_psMtpobs(i,j,bi,bj)=           anom_psMtpobs(i,j,bi,bj)=
1089       &       psbar(i,j,bi,bj) - tpobs(i,j,bi,bj)       &       psbar(i,j,bi,bj) - tpobs(i,j,bi,bj)
1090       &       - mean_psMtpobs(i,j,bi,bj)       &       - mean_psMslaobs(i,j,bi,bj) - slaoffset
1091           anom_tpobs(i,j,bi,bj)=tpobs(i,j,bi,bj)       &       - anom_psMslaobs(i,j,bi,bj)
1092             anom_tpobs(i,j,bi,bj)=tpobs(i,j,bi,bj) - slaoffset
1093         &       - anom_slaobs(i,j,bi,bj)
1094        endif        endif
1095  #endif  #endif
1096  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
1097        if ( ersmask(i,j,bi,bj)*mean_psMersobs_MSK(i,j,bi,bj).NE.0. )        if ( ersmask(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj).NE.0. )
1098       & then       & then
1099           anom_psMersobs(i,j,bi,bj)=           anom_psMersobs(i,j,bi,bj)=
1100       &       psbar(i,j,bi,bj) - ersobs(i,j,bi,bj)       &       psbar(i,j,bi,bj) - ersobs(i,j,bi,bj)
1101       &       - mean_psMersobs(i,j,bi,bj)       &       - mean_psMslaobs(i,j,bi,bj) - slaoffset
1102           anom_ersobs(i,j,bi,bj)=ersobs(i,j,bi,bj)       &       - anom_psMslaobs(i,j,bi,bj)
1103             anom_ersobs(i,j,bi,bj)=ersobs(i,j,bi,bj) - slaoffset
1104         &       - anom_slaobs(i,j,bi,bj)
1105        endif        endif
1106  #endif  #endif
1107  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
1108        if ( gfomask(i,j,bi,bj)*mean_psMgfoobs_MSK(i,j,bi,bj).NE.0. )        if ( gfomask(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj).NE.0. )
1109       & then       & then
1110           anom_psMgfoobs(i,j,bi,bj)=           anom_psMgfoobs(i,j,bi,bj)=
1111       &       psbar(i,j,bi,bj) - gfoobs(i,j,bi,bj)       &       psbar(i,j,bi,bj) - gfoobs(i,j,bi,bj)
1112       &       - mean_psMgfoobs(i,j,bi,bj)       &       - mean_psMslaobs(i,j,bi,bj) - slaoffset
1113           anom_gfoobs(i,j,bi,bj)=gfoobs(i,j,bi,bj)       &       - anom_psMslaobs(i,j,bi,bj)
1114             anom_gfoobs(i,j,bi,bj)=gfoobs(i,j,bi,bj) - slaoffset
1115         &       - anom_slaobs(i,j,bi,bj)
1116        endif        endif
1117  #endif  #endif
1118            enddo            enddo
# Line 935  cgf ==================================== Line 1140  cgf ====================================
1140        write(fname4test(1:80),'(1a)') 'slaobs_gfo_raw'        write(fname4test(1:80),'(1a)') 'slaobs_gfo_raw'
1141        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
1142       & 1,anom_gfoobs,irec,1,mythid)       & 1,anom_gfoobs,irec,1,mythid)
1143  #endif  #endif
1144    
1145         do bj = jtlo,jthi         do bj = jtlo,jthi
1146          do bi = itlo,ithi          do bi = itlo,ithi
# Line 943  cgf ==================================== Line 1148  cgf ====================================
1148            do i = imin,imax            do i = imin,imax
1149  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
1150  c--             The array psobs contains SSH anomalies.  c--             The array psobs contains SSH anomalies.
1151                  junkweight = mean_psMtpobs_MSK(i,j,bi,bj)                  junkweight = mean_psMslaobs_MSK(i,j,bi,bj)
1152       &                      *gencost_weight(i,j,bi,bj,igen_tp)       &                      *gencost_weight(i,j,bi,bj,igen_tp)
1153       &                      *tpmask(i,j,bi,bj)       &                      *tpmask(i,j,bi,bj)
1154       &                      *cosphi(i,j,bi,bj)       &                      *cosphi(i,j,bi,bj)
1155                  junk = anom_psMtpobs(i,j,bi,bj)                  junk = anom_psMtpobs(i,j,bi,bj)
1156                  objf_gencost(igen_tp,bi,bj) =                  objf_gencost(bi,bj,igen_tp) =
1157       &              objf_gencost(igen_tp,bi,bj)+junk*junk*junkweight       &              objf_gencost(bi,bj,igen_tp)+junk*junk*junkweight
1158                  if ( junkweight .ne. 0. )                  if ( junkweight .ne. 0. )
1159       &              num_gencost(igen_tp,bi,bj) =       &              num_gencost(bi,bj,igen_tp) =
1160       &              num_gencost(igen_tp,bi,bj) + 1. _d 0       &              num_gencost(bi,bj,igen_tp) + 1. _d 0
1161  #endif  #endif
1162  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
1163  c--             The array ersobs contains SSH anomalies.  c--             The array ersobs contains SSH anomalies.
1164                  junkweight = mean_psMersobs_MSK(i,j,bi,bj)                  junkweight = mean_psMslaobs_MSK(i,j,bi,bj)
1165       &                      *gencost_weight(i,j,bi,bj,igen_ers)       &                      *gencost_weight(i,j,bi,bj,igen_ers)
1166       &                      *ersmask(i,j,bi,bj)       &                      *ersmask(i,j,bi,bj)
1167       &                      *cosphi(i,j,bi,bj)       &                      *cosphi(i,j,bi,bj)
1168                  junk = anom_psMersobs(i,j,bi,bj)                  junk = anom_psMersobs(i,j,bi,bj)
1169                  objf_gencost(igen_ers,bi,bj) =                  objf_gencost(bi,bj,igen_ers) =
1170       &               objf_gencost(igen_ers,bi,bj)+junk*junk*junkweight       &               objf_gencost(bi,bj,igen_ers)+junk*junk*junkweight
1171                  if ( junkweight .ne. 0. )                  if ( junkweight .ne. 0. )
1172       &              num_gencost(igen_ers,bi,bj) =       &              num_gencost(bi,bj,igen_ers) =
1173       &              num_gencost(igen_ers,bi,bj) + 1. _d 0       &              num_gencost(bi,bj,igen_ers) + 1. _d 0
1174  #endif  #endif
1175  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
1176  c--             The array gfoobs contains SSH anomalies.  c--             The array gfoobs contains SSH anomalies.
1177                  junkweight = mean_psMgfoobs_MSK(i,j,bi,bj)                  junkweight = mean_psMslaobs_MSK(i,j,bi,bj)
1178       &                      *gencost_weight(i,j,bi,bj,igen_gfo)       &                      *gencost_weight(i,j,bi,bj,igen_gfo)
1179       &                      *gfomask(i,j,bi,bj)       &                      *gfomask(i,j,bi,bj)
1180       &                      *cosphi(i,j,bi,bj)       &                      *cosphi(i,j,bi,bj)
1181                  junk = anom_psMgfoobs(i,j,bi,bj)                  junk = anom_psMgfoobs(i,j,bi,bj)
1182                  objf_gencost(igen_gfo,bi,bj) =                  objf_gencost(bi,bj,igen_gfo) =
1183       &              objf_gencost(igen_gfo,bi,bj)+junk*junk*junkweight       &              objf_gencost(bi,bj,igen_gfo)+junk*junk*junkweight
1184                  if ( junkweight .ne. 0. )                  if ( junkweight .ne. 0. )
1185       &              num_gencost(igen_gfo,bi,bj) =       &              num_gencost(bi,bj,igen_gfo) =
1186       &              num_gencost(igen_gfo,bi,bj) + 1. _d 0       &              num_gencost(bi,bj,igen_gfo) + 1. _d 0
1187  #endif  #endif
1188             enddo             enddo
1189           enddo           enddo
# Line 988  c--             The array gfoobs contain Line 1193  c--             The array gfoobs contain
1193        enddo        enddo
1194    
1195    
 #endif /* ifdef ALLOW_GENCOST_FREEFORM */  
1196  #endif /* ifdef ALLOW_GENCOST_CONTRIBUTION */  #endif /* ifdef ALLOW_GENCOST_CONTRIBUTION */
 #endif /* ifdef ALLOW_SMOOTH */  
1197  #endif /* ifdef ALLOW_SSH_COST_CONTRIBUTION */  #endif /* ifdef ALLOW_SSH_COST_CONTRIBUTION */
1198    
1199        end        end
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.21
  ViewVC Help
Powered by ViewVC 1.1.22