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