/[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.22 by gforget, Mon Jun 9 17:53:59 2014 UTC revision 1.34 by atn, Thu Nov 12 12:11:20 2015 UTC
# Line 5  C $Name$ Line 5  C $Name$
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_SIZE.h"  #endif
46  #include "ctrl.h"  #ifdef ALLOW_ECCO
47  #include "ctrl_dummy.h"  # include "ecco.h"
 #include "optim.h"  
 #include "DYNVARS.h"  
 #ifdef ALLOW_PROFILES  
 #include "profiles.h"  
48  #endif  #endif
 #include "cal.h"  
49  #ifdef ALLOW_SMOOTH  #ifdef ALLOW_SMOOTH
50  #include "SMOOTH.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  
57  #ifdef ALLOW_GENCOST_CONTRIBUTION  #ifdef ALLOW_GENCOST_CONTRIBUTION
58    
59  c     == functions ==  c     == functions ==
# Line 80  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
# Line 88  c     == local variables == Line 77  c     == local variables ==
77  c mapping to gencost  c mapping to gencost
78        integer igen_mdt, igen_lsc, igen_gmsl        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 spval, factor
83        _RL offset        _RL offset
84        _RL offset_sum        _RL offset_sum
85        _RL slaoffset        _RL slaoffset
86        _RL slaoffset_sum        _RL slaoffset_sum
87        _RL slaoffset_weight        _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          _RL mdtob(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
97          _RL tpob(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
98          _RL ersob(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
99          _RL gfoob(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
100          _RL mdtma(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
101          _RL tpma(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
102          _RL ersma(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
103          _RL gfoma(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
104          _RL etaday(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
105    
106          character*(MAX_LEN_FNAM) mdtfi, tpfi, ersfi, gfofi
107          integer tpsta(4), erssta(4), gfosta(4)
108          integer mdtsta(4), mdtend(4)
109          _RL tpper, ersper, gfoper
110    
111  c for PART 1: re-reference MDT (mdt) to the inferred SLA reference field  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 )        _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)        _RL mean_slaobs_mdt(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
# Line 138  c for PART 4/5: compute smooth/raw anoma Line 149  c for PART 4/5: compute smooth/raw anoma
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 mdt_y0,mdt_y1,year,day        integer mdt_y0,mdt_y1,year,day
152        integer num_var, num0        integer num0
153    
154        _RL junk,junkweight        _RL junk,junkweight
155    
# Line 150  c for PART 4/5: compute smooth/raw anoma Line 161  c for PART 4/5: compute smooth/raw anoma
161        character*(MAX_LEN_MBUF) msgbuf        character*(MAX_LEN_MBUF) msgbuf
162    
163        LOGICAL doReference        LOGICAL doReference
164          LOGICAL useEtaMean
165    
166  c     == external functions ==  c     == external functions ==
167    
# Line 174  c-- detect the relevant gencost indices Line 186  c-- detect the relevant gencost indices
186        igen_gfo=0        igen_gfo=0
187        igen_lsc=0        igen_lsc=0
188        igen_gmsl=0        igen_gmsl=0
189          igen_etaday=0
190        do k=1,NGENCOST        do k=1,NGENCOST
191          if (gencost_name(k).EQ.'sshv4-mdt') igen_mdt=k          if (gencost_name(k).EQ.'sshv4-mdt') igen_mdt=k
192          if (gencost_name(k).EQ.'sshv4-tp') igen_tp=k          if (gencost_name(k).EQ.'sshv4-tp') igen_tp=k
# Line 183  c-- detect the relevant gencost indices Line 196  c-- detect the relevant gencost indices
196          if (gencost_name(k).EQ.'sshv4-gmsl') igen_gmsl=k          if (gencost_name(k).EQ.'sshv4-gmsl') igen_gmsl=k
197        enddo        enddo
198    
199    
200          call ecco_zero(gencost_weight(:,:,1,1,igen_mdt),1,zeroRL,myThid)
201          call ecco_zero(gencost_weight(:,:,1,1,igen_lsc),1,zeroRL,myThid)
202          call ecco_zero(gencost_weight(:,:,1,1,igen_tp),1,zeroRL,myThid)
203          call ecco_zero(gencost_weight(:,:,1,1,igen_ers),1,zeroRL,myThid)
204          call ecco_zero(gencost_weight(:,:,1,1,igen_gfo),1,zeroRL,myThid)
205          if ( gencost_errfile(igen_mdt) .NE. ' ' )
206         &     call ecco_readwei(gencost_errfile(igen_mdt),
207         &     gencost_weight(:,:,1,1,igen_mdt),1,1,mythid)
208          if ( gencost_errfile(igen_lsc) .NE. ' ' )
209         &     call ecco_readwei(gencost_errfile(igen_lsc),
210         &     gencost_weight(:,:,1,1,igen_lsc),1,1,mythid)
211          if ( gencost_errfile(igen_tp) .NE. ' ' )
212         &     call ecco_readwei(gencost_errfile(igen_tp),
213         &     gencost_weight(:,:,1,1,igen_tp),1,1,mythid)
214          if ( gencost_errfile(igen_ers) .NE. ' ' )
215         &     call ecco_readwei(gencost_errfile(igen_ers),
216         &     gencost_weight(:,:,1,1,igen_ers),1,1,mythid)
217          if ( gencost_errfile(igen_gfo) .NE. ' ' )
218         &     call ecco_readwei(gencost_errfile(igen_gfo),
219         &     gencost_weight(:,:,1,1,igen_gfo),1,1,mythid)
220    
221    c switch for excluding global mean
222          useEtaMean=.TRUE.
223    
224          write(msgbuf,'(a,l)')
225         & ' sshv4: use global mean of eta useEtaMean=',useEtaMean
226          call print_message( msgbuf, standardmessageunit,
227         &                    SQUEEZE_RIGHT , mythid)
228    
229    c minimum number of observations to consider re-referenced MDT
230          num0=100
231    c determine whether or not to re-reference
232    c-- not only model period has to fall within the mdt period
233    c-- but that there has to be at least 365 days of model run
234    c-- so for short model run, this will always get set to FALSE
235          doReference=.FALSE.
236          if ((modelstartdate(1).GT.1992*10000).AND.
237         &    (modelstartdate(1).LT.2011*10000).AND.
238         &    (ndaysrec.GE.365))  doReference=.TRUE.
239    
240          write(msgbuf,'(a,l)') ' sshv4:re-reference MDT=',doReference
241          call print_message( msgbuf, standardmessageunit,
242         &                    SQUEEZE_RIGHT , mythid)
243    
244    c-- initialize local arrays
245          do bj = jtlo,jthi
246            do bi = itlo,ithi
247              do j = jmin,jmax
248                do i = imin,imax
249                  mdtma(i,j,bi,bj) = 0. _d 0
250                  mdtob(i,j,bi,bj) = 0. _d 0
251                  tpma(i,j,bi,bj) = 0. _d 0
252                  tpob(i,j,bi,bj) = 0. _d 0
253                  ersma(i,j,bi,bj) = 0. _d 0
254                  ersob(i,j,bi,bj) = 0. _d 0
255                  gfoma(i,j,bi,bj) = 0. _d 0
256                  gfoob(i,j,bi,bj) = 0. _d 0
257                enddo
258              enddo
259            enddo
260          enddo
261    
262    c mdtfi, mdtsta, mdtend
263          if (igen_mdt.GT.0) then
264           ilps=ilnblnk( gencost_datafile(igen_mdt) )
265           mdtfi=gencost_datafile(igen_mdt)(1:ilps)
266           call cal_CopyDate(gencost_startdate(1,igen_mdt),mdtsta,mythid)
267           call cal_CopyDate(gencost_enddate(1,igen_mdt), mdtend, mythid)
268          endif
269    c tpfi, tpsta, tpper
270          if (igen_tp.GT.0) then
271           ilps=ilnblnk( gencost_datafile(igen_tp) )
272           tpfi=gencost_datafile(igen_tp)(1:ilps)
273           call cal_CopyDate(gencost_startdate(1,igen_tp),tpsta,mythid)
274           tpper=gencost_period(igen_tp)
275           if (gencost_barfile(igen_tp)(1:5).EQ.'m_eta')
276         &    igen_etaday=igen_tp
277          endif
278    c ersfi, erssta, ersper
279          if (igen_ers.GT.0) then
280           ilps=ilnblnk( gencost_datafile(igen_ers) )
281           ersfi=gencost_datafile(igen_ers)(1:ilps)
282           call cal_CopyDate(gencost_startdate(1,igen_ers),erssta,mythid)
283           ersper=gencost_period(igen_ers)
284           if (gencost_barfile(igen_ers)(1:5).EQ.'m_eta')
285         &    igen_etaday=igen_ers
286          endif
287    c gfofi, gfosta, gfoper
288          if (igen_gfo.GT.0) then
289           ilps=ilnblnk( gencost_datafile(igen_gfo) )
290           gfofi=gencost_datafile(igen_gfo)(1:ilps)
291           call cal_CopyDate(gencost_startdate(1,igen_gfo),gfosta,mythid)
292           gfoper=gencost_period(igen_gfo)
293           if (gencost_barfile(igen_gfo)(1:5).EQ.'m_eta')
294         &    igen_etaday=igen_gfo
295          endif
296    
297    c mdt_y0, mdt_y1
298           mdt_y0=mdtsta(1)/10000
299           mdt_y1=mdtend(1)/10000
300    
301           write(msgbuf,'(a,i8,a,i8)')
302         &  'mdt[start,end]date(1): ', mdtsta(1), ',', mdtend(1)
303            call print_message( msgbuf, standardmessageunit,
304         &                      SQUEEZE_RIGHT , mythid)
305            write(msgbuf,'(a,i4,a,i4)')
306         &  '  TP MDT yrrange:          ', mdt_y0, ',', mdt_y1
307            call print_message( msgbuf, standardmessageunit,
308         &                      SQUEEZE_RIGHT , mythid)
309    
310  c--   First, read tiled data.  c--   First, read tiled data.
311        doglobalread = .false.        doglobalread = .false.
312        ladinit      = .false.        ladinit      = .false.
313    
314        write(fname(1:80),'(80a)') ' '        write(fname(1:80),'(80a)') ' '
315        ilps=ilnblnk( psbarfile )        ilps=ilnblnk( gencost_barfile(igen_etaday) )
316        write(fname(1:80),'(2a,i10.10)')        write(fname(1:80),'(2a,i10.10)')
317       &     psbarfile(1:ilps),'.',optimcycle       &     gencost_barfile(igen_etaday)(1:ilps),'.',eccoiter
   
318    
319  cgf =======================================================  cgf =======================================================
320  cgf PART 1:  cgf PART 1:
# Line 201  cgf        together) over the time inter Line 324  cgf        together) over the time inter
324  cgf        mean_slaobs_mdt from mdt.  cgf        mean_slaobs_mdt from mdt.
325  cgf        x At this point, mdt is the inferred SLA reference field.  cgf        x At this point, mdt is the inferred SLA reference field.
326  cgf        x From there on, mdt+sla will be directly comparable to  cgf        x From there on, mdt+sla will be directly comparable to
327  cgf        the model SSH (psbar).  cgf        the model SSH (etaday).
328  cgf =======================================================  cgf =======================================================
329    
330  c--   Read mean field and mask  c--   Read mean field and mask
       call cost_ReadTopexMean( mythid )  
   
 c--   Compute mean_slaobs_mdt: sample mean SLA over the time period of mdt.  
   
 c pavlis and ecco/rio  
 c      mdt_y0=1993  
 c      mdt_y1=2004  
 c maximenko  
 c      mdt_y0=1992  
 c      mdt_y1=2002  
 c rio  
 c      mdt_y0=1993  
 c      mdt_y1=1999  
 c generalized:  
        mdt_y0=mdtstartdate(1)/10000  
        mdt_y1=mdtenddate(1)/10000  
   
        write(msgbuf,'(a,i8,a,i8)')  
      &  'mdt[start,end]date(1): ',mdtstartdate(1),  
      &  ',',mdtenddate(1)  
         call print_message( msgbuf, standardmessageunit,  
      &                      SQUEEZE_RIGHT , mythid)  
331    
332          write(msgbuf,'(a,i4,a,i4)')        if (using_mdt)
333       &  '  TP MDT yrrange:          ',       &call mdsreadfield( mdtfi, cost_iprec, cost_yftype, 1,
334       &                        mdt_y0,',',       &                   mdtob, 1, mythid )
      &                        mdt_y1  
         call print_message( msgbuf, standardmessageunit,  
      &                      SQUEEZE_RIGHT , mythid)  
335    
336  c minimum number of observations to consider re-referenced MDT        factor =  0.01 _d 0
337        num0=100        spval  = -9990. _d 0
338    
339        doReference=.FALSE.        do bj = jtlo,jthi
340        if ((modelstartdate(1).GT.1992*10000).AND.          do bi = itlo,ithi
341       &    (modelstartdate(1).LT.2011*10000).AND.            do j = jmin,jmax
342       &    (ndaysrec.GE.365))  doReference=.TRUE.              do i = imin,imax
343                  if ( (maskC(i,j,1,bi,bj) .eq. 0.).OR.
344    #ifndef ALLOW_SHALLOW_ALTIMETRY
345         &             (R_low(i,j,bi,bj).GT.-200.).OR.
346    #endif
347         &             (mdtob(i,j,bi,bj) .lt. spval ).OR.
348         &             (mdtob(i,j,bi,bj) .eq. 0. _d 0) ) then
349                    mdtma(i,j,bi,bj) = 0. _d 0
350                    mdtob(i,j,bi,bj) = 0. _d 0
351                  else
352                    mdtma(i,j,bi,bj) = 1. _d 0
353                    mdtob(i,j,bi,bj) = mdtob(i,j,bi,bj)*factor
354                  endif
355                enddo
356              enddo
357            enddo
358          enddo
359    
360          write(msgbuf,'(a,l)') ' sshv4:re-reference MDT=',doReference  c--   Compute mean_slaobs_mdt: sample mean SLA over the time period of mdt.
         call print_message( msgbuf, standardmessageunit,  
      &                          SQUEEZE_RIGHT , mythid)  
361    
362         do bj = jtlo,jthi         do bj = jtlo,jthi
363          do bi = itlo,ithi          do bi = itlo,ithi
# Line 260  c minimum number of observations to cons Line 372  c minimum number of observations to cons
372    
373        do year=mdt_y0,mdt_y1        do year=mdt_y0,mdt_y1
374         do day=1,366         do day=1,366
375  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  
376        call cost_sla_read_yd( topexfile, topexstartdate,         do bj = jtlo,jthi
377       &                tpobs, tpmask,          do bi = itlo,ithi
378             do j = jmin,jmax
379              do i = imin,imax
380                  tpma(i,j,bi,bj) = 0. _d 0
381                  tpob(i,j,bi,bj) = 0. _d 0
382                  ersma(i,j,bi,bj) = 0. _d 0
383                  ersob(i,j,bi,bj) = 0. _d 0
384                  gfoma(i,j,bi,bj) = 0. _d 0
385                  gfoob(i,j,bi,bj) = 0. _d 0
386              enddo
387             enddo
388            enddo
389           enddo
390    
391           if (using_tpj)
392         & call cost_sla_read_yd( tpfi, tpsta,
393         &                tpob, tpma,
394       &                year, day, mythid )       &                year, day, mythid )
395  #endif         if (using_ers)
396  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION       & call cost_sla_read_yd( ersfi, erssta,
397        call cost_sla_read_yd( ersfile, ersstartdate,       &                ersob, ersma,
      &                ersobs, ersmask,  
398       &                year, day, mythid )       &                year, day, mythid )
399  #endif         if (using_gfo)
400  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION       & call cost_sla_read_yd( gfofi, gfosta,
401        call cost_sla_read_yd( gfofile, gfostartdate,       &                gfoob, gfoma,
      &                gfoobs, gfomask,  
402       &                year, day, mythid )       &                year, day, mythid )
403  #endif  
404         do bj = jtlo,jthi         do bj = jtlo,jthi
405          do bi = itlo,ithi          do bi = itlo,ithi
406           do j = jmin,jmax           do j = jmin,jmax
407            do i = imin,imax            do i = imin,imax
408  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION        if ( tpma(i,j,bi,bj)*mdtma(i,j,bi,bj)*
       if ( tpmask(i,j,bi,bj)*mdtmask(i,j,bi,bj)*  
409       &    gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then       &    gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then
410            mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+            mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+
411       &  tpobs(i,j,bi,bj)       &  tpob(i,j,bi,bj)
412            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
413        endif        endif
414  #endif        if ( ersma(i,j,bi,bj)*mdtma(i,j,bi,bj)*
 #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION  
       if ( ersmask(i,j,bi,bj)*mdtmask(i,j,bi,bj)*  
415       &    gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then       &    gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then
416            mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+            mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+
417       &  ersobs(i,j,bi,bj)       &  ersob(i,j,bi,bj)
418            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
419        endif        endif
420  #endif        if ( gfoma(i,j,bi,bj)*mdtma(i,j,bi,bj)*
 #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION  
       if ( gfomask(i,j,bi,bj)*mdtmask(i,j,bi,bj)*  
421       &    gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then       &    gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then
422            mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+            mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+
423       &  gfoobs(i,j,bi,bj)       &  gfoob(i,j,bi,bj)
424            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
425        endif        endif
 #endif  
426            enddo            enddo
427           enddo           enddo
428          enddo          enddo
# Line 317  c minimum number of observations to cons Line 437  c minimum number of observations to cons
437              do i = imin,imax              do i = imin,imax
438                 if ( ( mean_slaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.                 if ( ( mean_slaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
439       &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.       &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
440  #ifndef ALLOW_HIGHLAT_ALTIMETRY       &              ( mdtma(i,j,bi,bj) .NE. 0. ) ) then
      &              ( abs(YC(i,j,bi,bj)) .LE. 66. ).AND.  
 #endif  
      &              ( mdtmask(i,j,bi,bj) .NE. 0. ) ) then  
441                    mean_slaobs_mdt(i,j,bi,bj) =                    mean_slaobs_mdt(i,j,bi,bj) =
442       &                 mean_slaobs_mdt(i,j,bi,bj) /       &                 mean_slaobs_mdt(i,j,bi,bj) /
443       &                 mean_slaobs_NUM(i,j,bi,bj)       &                 mean_slaobs_NUM(i,j,bi,bj)
# Line 335  c minimum number of observations to cons Line 452  c minimum number of observations to cons
452    
453  c--   smooth mean_slaobs_mdt:  c--   smooth mean_slaobs_mdt:
454    
455          if (gencost_outputlevel(igen_mdt).GT.0) then
456        write(fname4test(1:80),'(1a)') 'sla2mdt_raw'        write(fname4test(1:80),'(1a)') 'sla2mdt_raw'
457        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
458       & 1,mean_slaobs_mdt,1,1,mythid)       & 1,mean_slaobs_mdt,1,1,mythid)
459          endif
460    
461  #ifdef ALLOW_SMOOTH  #ifdef ALLOW_SMOOTH
462        if ( useSMOOTH )        if ( useSMOOTH )
# Line 346  c--   smooth mean_slaobs_mdt: Line 465  c--   smooth mean_slaobs_mdt:
465       &     gencost_smooth2Ddiffnbt(igen_mdt),mythid)       &     gencost_smooth2Ddiffnbt(igen_mdt),mythid)
466  #endif  #endif
467    
468          if (gencost_outputlevel(igen_mdt).GT.0) then
469        write(fname4test(1:80),'(1a)') 'sla2mdt_smooth'        write(fname4test(1:80),'(1a)') 'sla2mdt_smooth'
470        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
471       & 1,mean_slaobs_mdt,1,1,mythid)       & 1,mean_slaobs_mdt,1,1,mythid)
472          endif
473    
474  c--   re-reference mdt:  c--   re-reference mdt:
475        do bj = jtlo,jthi        do bj = jtlo,jthi
476          do bi = itlo,ithi          do bi = itlo,ithi
477            do j = jmin,jmax            do j = jmin,jmax
478              do i = imin,imax              do i = imin,imax
479                 if ( ( mdtmask(i,j,bi,bj) .NE. 0. ).AND.                 if ( ( mdtma(i,j,bi,bj) .NE. 0. ).AND.
480       &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.       &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
481       &              ( doReference ) ) then       &              ( doReference ) ) then
482                    mdt(i,j,bi,bj) = mdt(i,j,bi,bj)                    mdtob(i,j,bi,bj) = mdtob(i,j,bi,bj)
483       &                 -mean_slaobs_mdt(i,j,bi,bj)       &                 -mean_slaobs_mdt(i,j,bi,bj)
484                 endif                 endif
485              enddo              enddo
# Line 368  c--   re-reference mdt: Line 489  c--   re-reference mdt:
489    
490    
491  cgf =======================================================  cgf =======================================================
492  cgf PART 2: compute sample means of psbar-slaobs over the  cgf PART 2: compute sample means of etaday-slaobs over the
493  cgf          period that is covered by the model (i.e. psbar).  cgf          period that is covered by the model (i.e. etaday).
494  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,
495  cgf             and offset will be used in PART 3 (MDT cost term).  cgf             and offset will be used in PART 3 (MDT cost term).
496  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 405  cgf ==================================== Line 526  cgf ====================================
526    
527        do irec = 1, ndaysrec        do irec = 1, ndaysrec
528    
529          call active_read_xy( fname, psbar, irec, doglobalread,  #ifdef ALLOW_AUTODIFF
530       &                       ladinit, optimcycle, mythid,          call active_read_xy( fname, etaday, irec, doglobalread,
531       &                       xx_psbar_mean_dummy )       &                       ladinit, eccoiter, mythid,
532         &                       gencost_dummy(igen_etaday) )
533  #ifndef ALLOW_PSBAR_MEAN  #else
534          CALL REMOVE_MEAN_RL( 1, psbar, maskInC, maskInC, rA, drF,          CALL READ_REC_XY_RL( fname, etaday, iRec, 1, myThid )
      &        'psbar', myTime, myThid )  
535  #endif  #endif
536    
537  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  
538        call cost_sla_read( topexfile, topexstartdate, topexperiod,          if (.NOT.useEtaMean)
539       &                topexintercept, topexslope,       &  CALL REMOVE_MEAN_RL( 1, etaday, maskInC, maskInC, rA, drF,
540       &                tpobs, tpmask,       &        'etaday', 0. _d 0, myThid )
541    
542           do bj = jtlo,jthi
543            do bi = itlo,ithi
544             do j = jmin,jmax
545              do i = imin,imax
546                  tpma(i,j,bi,bj) = 0. _d 0
547                  tpob(i,j,bi,bj) = 0. _d 0
548                  ersma(i,j,bi,bj) = 0. _d 0
549                  ersob(i,j,bi,bj) = 0. _d 0
550                  gfoma(i,j,bi,bj) = 0. _d 0
551                  gfoob(i,j,bi,bj) = 0. _d 0
552              enddo
553             enddo
554            enddo
555           enddo
556    
557            if (using_tpj)
558         &  call cost_sla_read( tpfi, tpsta, tpper,
559         &                zeroRL, zeroRL,
560         &                tpob, tpma,
561       &                irec, mythid )       &                irec, mythid )
562  #endif          if (using_ers)
563  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION       &  call cost_sla_read( ersfi, erssta, ersper,
564        call cost_sla_read( ersfile, ersstartdate, ersperiod,       &                zeroRL, zeroRL,
565       &                ersintercept, ersslope,       &                ersob, ersma,
      &                ersobs, ersmask,  
566       &                irec, mythid )       &                irec, mythid )
567  #endif          if (using_gfo)
568  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION       &  call cost_sla_read( gfofi, gfosta, gfoper,
569        call cost_sla_read( gfofile, gfostartdate, gfoperiod,       &                zeroRL, zeroRL,
570       &                gfointercept, gfoslope,       &                gfoob, gfoma,
      &                gfoobs, gfomask,  
571       &                irec, mythid )       &                irec, mythid )
 #endif  
572    
573          do bj = jtlo,jthi          do bj = jtlo,jthi
574            do bi = itlo,ithi            do bi = itlo,ithi
575              do j = jmin,jmax              do j = jmin,jmax
576                do i = imin,imax                do i = imin,imax
577                  psmean(i,j,bi,bj) = psmean(i,j,bi,bj) +                  psmean(i,j,bi,bj) = psmean(i,j,bi,bj) +
578       &                psbar(i,j,bi,bj) / float(ndaysrec)       &                etaday(i,j,bi,bj) / float(ndaysrec)
579  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION                  if ( tpma(i,j,bi,bj)*
                 if ( tpmask(i,j,bi,bj)*  
580       &             gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then       &             gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then
581                     mean_slaobs_model(i,j,bi,bj)=                     mean_slaobs_model(i,j,bi,bj)=
582       &                 mean_slaobs_model(i,j,bi,bj)+tpobs(i,j,bi,bj)       &              mean_slaobs_model(i,j,bi,bj)+tpob(i,j,bi,bj)
583                     mean_psMtpobs(i,j,bi,bj) =                     mean_psMtpobs(i,j,bi,bj) =
584       &                 mean_psMtpobs(i,j,bi,bj) +       &                 mean_psMtpobs(i,j,bi,bj) +
585       &                 psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)       &                 etaday(i,j,bi,bj)-tpob(i,j,bi,bj)
586                     mean_psMtpobs_NUM(i,j,bi,bj) =                     mean_psMtpobs_NUM(i,j,bi,bj) =
587       &                 mean_psMtpobs_NUM(i,j,bi,bj) + 1. _d 0       &                 mean_psMtpobs_NUM(i,j,bi,bj) + 1. _d 0
588                  endif                  endif
589  #endif                  if ( ersma(i,j,bi,bj)*
 #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION  
                 if ( ersmask(i,j,bi,bj)*  
590       &             gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then       &             gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then
591                     mean_slaobs_model(i,j,bi,bj)=                     mean_slaobs_model(i,j,bi,bj)=
592       &                 mean_slaobs_model(i,j,bi,bj)+ersobs(i,j,bi,bj)       &              mean_slaobs_model(i,j,bi,bj)+ersob(i,j,bi,bj)
593                     mean_psMersobs(i,j,bi,bj) =                     mean_psMersobs(i,j,bi,bj) =
594       &                 mean_psMersobs(i,j,bi,bj) +       &                 mean_psMersobs(i,j,bi,bj) +
595       &                 psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)       &                 etaday(i,j,bi,bj)-ersob(i,j,bi,bj)
596                     mean_psMersobs_NUM(i,j,bi,bj) =                     mean_psMersobs_NUM(i,j,bi,bj) =
597       &                 mean_psMersobs_NUM(i,j,bi,bj) + 1. _d 0       &                 mean_psMersobs_NUM(i,j,bi,bj) + 1. _d 0
598                  endif                  endif
599  #endif                  if ( gfoma(i,j,bi,bj)*
 #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION  
                 if ( gfomask(i,j,bi,bj)*  
600       &             gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then       &             gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then
601                     mean_slaobs_model(i,j,bi,bj)=                     mean_slaobs_model(i,j,bi,bj)=
602       &                 mean_slaobs_model(i,j,bi,bj)+gfoobs(i,j,bi,bj)       &              mean_slaobs_model(i,j,bi,bj)+gfoob(i,j,bi,bj)
603                     mean_psMgfoobs(i,j,bi,bj) =                     mean_psMgfoobs(i,j,bi,bj) =
604       &                 mean_psMgfoobs(i,j,bi,bj) +       &                 mean_psMgfoobs(i,j,bi,bj) +
605       &                 psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)       &                 etaday(i,j,bi,bj)-gfoob(i,j,bi,bj)
606                     mean_psMgfoobs_NUM(i,j,bi,bj) =                     mean_psMgfoobs_NUM(i,j,bi,bj) =
607       &                 mean_psMgfoobs_NUM(i,j,bi,bj) + 1. _d 0       &                 mean_psMgfoobs_NUM(i,j,bi,bj) + 1. _d 0
608                  endif                  endif
 #endif  
609                enddo                enddo
610              enddo              enddo
611            enddo            enddo
# Line 483  cgf ==================================== Line 614  cgf ====================================
614  c--   END loop over records for the first time.  c--   END loop over records for the first time.
615        enddo        enddo
616    
617          call ecco_zero(num0array,1,oneRL,mythid)
618    
619          do bj = jtlo,jthi          do bj = jtlo,jthi
620            do bi = itlo,ithi            do bi = itlo,ithi
621              do j = jmin,jmax              do j = jmin,jmax
622                do i = imin,imax                do i = imin,imax
 #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  
623                 if ( ( mean_psMtpobs_NUM(i,j,bi,bj) .NE. 0. )                 if ( ( mean_psMtpobs_NUM(i,j,bi,bj) .NE. 0. )
 #ifndef ALLOW_HIGHLAT_ALTIMETRY  
      &              .AND.( abs(YC(i,j,bi,bj)) .LE. 66. )  
 #endif  
624       &              ) then       &              ) then
625                    mean_psMssh_all(i,j,bi,bj) =                    mean_psMssh_all(i,j,bi,bj) =
626       &                 mean_psMssh_all(i,j,bi,bj) +       &                 mean_psMssh_all(i,j,bi,bj) +
# Line 503  c--   END loop over records for the firs Line 632  c--   END loop over records for the firs
632       &                 mean_psMtpobs(i,j,bi,bj) /       &                 mean_psMtpobs(i,j,bi,bj) /
633       &                 mean_psMtpobs_NUM(i,j,bi,bj)       &                 mean_psMtpobs_NUM(i,j,bi,bj)
634                 endif                 endif
 #endif  
 #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION  
635                 if ( ( mean_psMersobs_NUM(i,j,bi,bj) .NE. 0. )                 if ( ( mean_psMersobs_NUM(i,j,bi,bj) .NE. 0. )
 #ifndef ALLOW_HIGHLAT_ALTIMETRY      
      &              .AND.( abs(YC(i,j,bi,bj)) .LE. 66. )  
 #endif  
636       &              ) then       &              ) then
637                    mean_psMssh_all(i,j,bi,bj) =                    mean_psMssh_all(i,j,bi,bj) =
638       &                 mean_psMssh_all(i,j,bi,bj) +       &                 mean_psMssh_all(i,j,bi,bj) +
# Line 520  c--   END loop over records for the firs Line 644  c--   END loop over records for the firs
644       &                 mean_psMersobs(i,j,bi,bj) /       &                 mean_psMersobs(i,j,bi,bj) /
645       &                 mean_psMersobs_NUM(i,j,bi,bj)       &                 mean_psMersobs_NUM(i,j,bi,bj)
646                 endif                 endif
 #endif  
 #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION  
647                 if ( ( mean_psMgfoobs_NUM(i,j,bi,bj) .NE. 0. )                 if ( ( mean_psMgfoobs_NUM(i,j,bi,bj) .NE. 0. )
 #ifndef ALLOW_HIGHLAT_ALTIMETRY      
      &              .AND.( abs(YC(i,j,bi,bj)) .LE. 66. )  
 #endif  
648       &              ) then       &              ) then
649                    mean_psMssh_all(i,j,bi,bj) =                    mean_psMssh_all(i,j,bi,bj) =
650       &                 mean_psMssh_all(i,j,bi,bj) +       &                 mean_psMssh_all(i,j,bi,bj) +
# Line 537  c--   END loop over records for the firs Line 656  c--   END loop over records for the firs
656       &                 mean_psMgfoobs(i,j,bi,bj) /       &                 mean_psMgfoobs(i,j,bi,bj) /
657       &                 mean_psMgfoobs_NUM(i,j,bi,bj)       &                 mean_psMgfoobs_NUM(i,j,bi,bj)
658                 endif                 endif
 #endif  
659    
660                 if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .GT. 0 ).AND.                 if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .GT. 0 ).AND.
661       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.       &              ( maskc(i,j,1,bi,bj) .NE. 0. )
 #ifndef ALLOW_HIGHLAT_ALTIMETRY      
      &              ( abs(YC(i,j,bi,bj)) .LE. 66. )  
 #endif  
662       &              ) then       &              ) then
663                    mean_psMslaobs(i,j,bi,bj) =                    mean_psMslaobs(i,j,bi,bj) =
664       &                 mean_psMssh_all(i,j,bi,bj) /       &                 mean_psMssh_all(i,j,bi,bj) /
# Line 556  c--   END loop over records for the firs Line 671  c--   END loop over records for the firs
671    
672                 if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .GT. num0 ).AND.                 if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .GT. num0 ).AND.
673       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
674  #ifndef ALLOW_HIGHLAT_ALTIMETRY           &              ( mdtma(i,j,bi,bj) .NE. 0. ).AND.
      &              ( abs(YC(i,j,bi,bj)) .LE. 66. ).AND.  
 #endif  
      &              ( mdtmask(i,j,bi,bj) .NE. 0. ).AND.  
675       &              ( doReference ) ) then       &              ( doReference ) ) then
676                    mean_slaobs_model(i,j,bi,bj) =                    mean_slaobs_model(i,j,bi,bj) =
677       &                 mean_slaobs_model(i,j,bi,bj) /       &                 mean_slaobs_model(i,j,bi,bj) /
678       &                 mean_psMssh_all_NUM(i,j,bi,bj)       &                 mean_psMssh_all_NUM(i,j,bi,bj)
679                    mean_psMssh_all(i,j,bi,bj) =                    mean_psMssh_all(i,j,bi,bj) =
680       &                 mean_psMssh_all(i,j,bi,bj) /       &                 mean_psMssh_all(i,j,bi,bj) /
681       &                 mean_psMssh_all_NUM(i,j,bi,bj)-mdt(i,j,bi,bj)       &                 mean_psMssh_all_NUM(i,j,bi,bj)-mdtob(i,j,bi,bj)
682                    mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0                    mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
683                    offset=offset+RA(i,j,bi,bj)*mean_psMssh_all(i,j,bi,bj)                    offset=offset+RA(i,j,bi,bj)*mean_psMssh_all(i,j,bi,bj)
684                    offset_sum=offset_sum+RA(i,j,bi,bj)                    offset_sum=offset_sum+RA(i,j,bi,bj)
685                 elseif ( ( mdtmask(i,j,bi,bj) .NE. 0. ) .AND.                 elseif ( ( mdtma(i,j,bi,bj) .NE. 0. ) .AND.
686       &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.       &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
 #ifndef ALLOW_HIGHLAT_ALTIMETRY  
      &              ( abs(YC(i,j,bi,bj)) .LE. 66. ).AND.  
 #endif  
687       &              ( .NOT.doReference ) ) then       &              ( .NOT.doReference ) ) then
688                    mean_slaobs_model(i,j,bi,bj) = 0.d0                    mean_slaobs_model(i,j,bi,bj) = 0.d0
689                    mean_psMssh_all(i,j,bi,bj) = 0. _d 0                    mean_psMssh_all(i,j,bi,bj) = 0. _d 0
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+RA(i,j,bi,bj)                    offset=offset+RA(i,j,bi,bj)
692       &                  *( psmean(i,j,bi,bj) -mdt(i,j,bi,bj) )       &                  *( psmean(i,j,bi,bj) -mdtob(i,j,bi,bj) )
693                    offset_sum=offset_sum+RA(i,j,bi,bj)                    offset_sum=offset_sum+RA(i,j,bi,bj)
694                      num0array(i,j,bi,bj)=0. _d 0
695                 else                 else
696                    mean_slaobs_model(i,j,bi,bj) = 0.d0                    mean_slaobs_model(i,j,bi,bj) = 0.d0
697                    mean_psMssh_all(i,j,bi,bj) = 0. _d 0                    mean_psMssh_all(i,j,bi,bj) = 0. _d 0
698                    mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0                    mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0
699                      num0array(i,j,bi,bj)=0. _d 0
700                 endif                 endif
701                enddo                enddo
702              enddo              enddo
# Line 597  c--   Do a global summation. Line 708  c--   Do a global summation.
708        _GLOBAL_SUM_RL( offset     , mythid )        _GLOBAL_SUM_RL( offset     , mythid )
709        _GLOBAL_SUM_RL( offset_sum , mythid )        _GLOBAL_SUM_RL( offset_sum , mythid )
710    
711    catn--- add warning that written mean_slaobs_model is all zero
712    c       due to having less data points that hardcoded num0
713          num0total=0. _d 0
714            do bj = jtlo,jthi
715              do bi = itlo,ithi
716                do j = jmin,jmax
717                  do i = imin,imax
718                    num0total=num0total+num0array(i,j,bi,bj)
719                  enddo
720                enddo
721              enddo
722            enddo
723          if(num0total.lt.1. _d 0) then
724            write(msgbuf,'(A,i5,A)')
725         &  '** WARNING ** S/R COST_GENCOST_SSHV4: There are <',num0,' pts'
726          call print_message( msgbuf, errormessageunit,
727         &                    SQUEEZE_RIGHT , mythid)
728            write(msgbuf,'(A)')
729         &  '                  at all grid points for combined tp+ers+gfo.'
730          call print_message( msgbuf, errormessageunit,
731         &                    SQUEEZE_RIGHT , mythid)
732            write(msgbuf,'(A)')
733         &  'So, all model slaob minus model sla2model_raw are set to 0.'
734          call print_message( msgbuf, errormessageunit,
735         &                    SQUEEZE_RIGHT , mythid)
736          endif
737    catn-------
738    
739        if (offset_sum .eq. 0.0) then        if (offset_sum .eq. 0.0) then
740            if (gencost_outputlevel(igen_gmsl).GT.0) then
741          _BEGIN_MASTER( mythid )          _BEGIN_MASTER( mythid )
742          write(msgbuf,'(a)') ' sshv4: offset_sum = zero!'          write(msgbuf,'(a)') ' sshv4: offset_sum = zero!'
743          call print_message( msgbuf, standardmessageunit,          call print_message( msgbuf, standardmessageunit,
744       &                          SQUEEZE_RIGHT , mythid)       &                          SQUEEZE_RIGHT , mythid)
745          _END_MASTER( mythid )          _END_MASTER( mythid )
746            endif
747  c        stop   '  ... stopped in cost_ssh.'  c        stop   '  ... stopped in cost_ssh.'
748        else        else
749            if (gencost_outputlevel(igen_gmsl).GT.0) then
750          _BEGIN_MASTER( mythid )          _BEGIN_MASTER( mythid )
751          write(msgbuf,'(a,2d22.15)') ' sshv4:offset,sum=',          write(msgbuf,'(a,2d22.15)') ' sshv4:offset,sum=',
752       &      offset,offset_sum       &      offset,offset_sum
753          call print_message( msgbuf, standardmessageunit,          call print_message( msgbuf, standardmessageunit,
754       &                          SQUEEZE_RIGHT , mythid)       &                          SQUEEZE_RIGHT , mythid)
755          _END_MASTER( mythid )          _END_MASTER( mythid )
756            endif
757        endif        endif
758    
759  c--   Compute (average) offset  c--   Compute (average) offset
# Line 623  c--   subtract offset from mean_psMssh_a Line 766  c--   subtract offset from mean_psMssh_a
766              do i = imin,imax              do i = imin,imax
767    
768                 if ( (mean_psMssh_all_NUM(i,j,bi,bj) .GT. num0) .AND.                 if ( (mean_psMssh_all_NUM(i,j,bi,bj) .GT. num0) .AND.
769       &              ( mdtmask(i,j,bi,bj) .NE. 0. ) .AND.       &              ( mdtma(i,j,bi,bj) .NE. 0. ) .AND.
770       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
 #ifndef ALLOW_HIGHLAT_ALTIMETRY      
      &              ( abs(YC(i,j,bi,bj)) .LE. 66. ) .AND.  
 #endif  
771       &              ( doReference ) ) then       &              ( doReference ) ) then
772  c use the re-referencing approach  c use the re-referencing approach
773                     mean_psMssh_all(i,j,bi,bj) =                     mean_psMssh_all(i,j,bi,bj) =
774       &                 mean_psMssh_all(i,j,bi,bj) - offset       &                 mean_psMssh_all(i,j,bi,bj) - offset
775                     mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0                     mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
776    
777                 elseif ( ( mdtmask(i,j,bi,bj) .NE. 0. ) .AND.                 elseif ( ( mdtma(i,j,bi,bj) .NE. 0. ) .AND.
778       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.       &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
779  #ifndef ALLOW_HIGHLAT_ALTIMETRY       &              ( .NOT.doReference ) ) then
      &              ( abs(YC(i,j,bi,bj)) .LE. 66. ) .AND.  
 #endif  
      &              ( .NOT.doReference ) ) then      
780  c use the simpler approach  c use the simpler approach
781                     mean_psMssh_all(i,j,bi,bj) =                     mean_psMssh_all(i,j,bi,bj) =
782       &             psmean(i,j,bi,bj) -mdt(i,j,bi,bj) - offset       &             psmean(i,j,bi,bj) -mdtob(i,j,bi,bj) - offset
783                     mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0                     mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
784    
785                 else                 else
# Line 659  c use the simpler approach Line 796  c use the simpler approach
796         enddo         enddo
797    
798  c--    smooth mean_psMssh_all  c--    smooth mean_psMssh_all
799          if (gencost_outputlevel(igen_mdt).GT.0) then
800        write(fname4test(1:80),'(1a)') 'mdtdiff_raw'        write(fname4test(1:80),'(1a)') 'mdtdiff_raw'
801        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
802       &    1,mean_psMssh_all,1,1,mythid)       &    1,mean_psMssh_all,1,1,mythid)
803          endif
804    
805  #ifdef ALLOW_SMOOTH  #ifdef ALLOW_SMOOTH
806        if ( useSMOOTH )        if ( useSMOOTH )
# Line 670  c--    smooth mean_psMssh_all Line 809  c--    smooth mean_psMssh_all
809       &     gencost_smooth2Ddiffnbt(igen_mdt),mythid)       &     gencost_smooth2Ddiffnbt(igen_mdt),mythid)
810  #endif  #endif
811    
812          if (gencost_outputlevel(igen_mdt).GT.0) then
813        write(fname4test(1:80),'(1a)') 'mdtdiff_smooth'        write(fname4test(1:80),'(1a)') 'mdtdiff_smooth'
814        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
815       &    1,mean_psMssh_all,1,1,mythid)       &    1,mean_psMssh_all,1,1,mythid)
816          endif
817    
818          if (gencost_outputlevel(igen_mdt).GT.0) then
819        write(fname4test(1:80),'(1a)') 'sla2model_raw'        write(fname4test(1:80),'(1a)') 'sla2model_raw'
820        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
821       &    1,mean_slaobs_model,1,1,mythid)       &    1,mean_slaobs_model,1,1,mythid)
822          endif
823    
824  #ifdef ALLOW_SMOOTH  #ifdef ALLOW_SMOOTH
825        if ( useSMOOTH )        if ( useSMOOTH )
# Line 685  c--    smooth mean_psMssh_all Line 828  c--    smooth mean_psMssh_all
828       &     gencost_smooth2Ddiffnbt(igen_mdt),mythid)       &     gencost_smooth2Ddiffnbt(igen_mdt),mythid)
829  #endif  #endif
830    
831          if (gencost_outputlevel(igen_mdt).GT.0) then
832        write(fname4test(1:80),'(1a)') 'sla2model_smooth'        write(fname4test(1:80),'(1a)') 'sla2model_smooth'
833        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
834       &    1,mean_slaobs_model,1,1,mythid)       &    1,mean_slaobs_model,1,1,mythid)
835          endif
836    
837  cgf at this point:  cgf at this point:
838  cgf     1) mean_psMssh_all is the sample mean <psbar-slaobs-mdt-offset>,  cgf     1) mean_psMssh_all is the sample mean <etaday-slaobs-mdt-offset>,
839  cgf             to which a smoothing filter has been applied.  cgf             to which a smoothing filter has been applied.
840  cgf     2) mean_psMtpobs is the (unsmoothed) sample mean <psbar-tpobs>.  cgf     2) mean_psMtpobs is the (unsmoothed) sample mean <etaday-tpobs>.
841  cgf             And similarly for ers and gfo, each treated separately.  cgf             And similarly for ers and gfo, each treated separately.
842    
 #ifdef ALLOW_PROFILES  
       if ( usePROFILES ) then  
       do bj = jtlo,jthi  
         do bi = itlo,ithi  
           do j = 1,sny  
             do i = 1,snx  
 cgf this has not been tested recently  
               prof_etan_mean(i,j,bi,bj)=psmean(i,j,bi,bj)  
             enddo  
           enddo  
         enddo  
       enddo  
       _EXCH_XY_RL( prof_etan_mean, mythid )  
       endif  
 #endif  
   
   
843  cgf =======================================================  cgf =======================================================
844  cgf PART 3: compute MDT cost term  cgf PART 3: compute MDT cost term
845  cgf =======================================================  cgf =======================================================
846    
847    
 #ifdef ALLOW_SSH_MEAN_COST_CONTRIBUTION  
   
848        do bj = jtlo,jthi        do bj = jtlo,jthi
849          do bi = itlo,ithi          do bi = itlo,ithi
850            do j = jmin,jmax            do j = jmin,jmax
# Line 726  cgf ==================================== Line 852  cgf ====================================
852         if (mean_psMssh_all_MSK(i,j,bi,bj).NE.0. _d 0) then         if (mean_psMssh_all_MSK(i,j,bi,bj).NE.0. _d 0) then
853           junk = mean_psMssh_all(i,j,bi,bj)           junk = mean_psMssh_all(i,j,bi,bj)
854           junkweight = gencost_weight(i,j,bi,bj,igen_mdt)           junkweight = gencost_weight(i,j,bi,bj,igen_mdt)
855       &              * mdtmask(i,j,bi,bj)       &              * mdtma(i,j,bi,bj)
856         else         else
857           junk = 0. _d 0           junk = 0. _d 0
858           junkweight = 0. _d 0           junkweight = 0. _d 0
# Line 741  cgf ==================================== Line 867  cgf ====================================
867          enddo          enddo
868        enddo        enddo
869    
870        CALL WRITE_FLD_XY_RL( 'DiagnosMDT', ' ', diagnosfld,        if (gencost_outputlevel(igen_mdt).GT.0) then
871       &                           optimcycle, mythid )        write(fname4test(1:80),'(1a)') 'misfits_mdt'
872          call mdswritefield(fname4test,32,.false.,'RL',
873  #endif /* ALLOW_SSH_MEAN_COST_CONTRIBUTION */       & 1,diagnosfld,1,1,mythid)
874          endif
   
875    
876  cgf =======================================================  cgf =======================================================
877  cgf PART 4: compute smooth SLA cost term  cgf PART 4: compute smooth SLA cost term
# Line 756  cgf ==================================== Line 881  cgf ====================================
881        ndaysave=35        ndaysave=35
882        ndaysaveRL=ndaysave        ndaysaveRL=ndaysave
883    
884    catn add a warning if have less nrecday than the hardcoded 35days
885          tempinteger=ndaysrec-ndaysave+1
886          if(tempinteger.lt.1) then
887            write(msgbuf,'(A,i5)')
888         &  '** WARNING ** S/R COST_GENCOST_SSHV4: There are < ',ndaysave
889          call print_message( msgbuf, errormessageunit,
890         &                    SQUEEZE_RIGHT , mythid)
891            write(msgbuf,'(A)')
892         &   '        days required to calculate running mean tp+ers+gfo.'
893          call print_message( msgbuf, errormessageunit,
894         &                    SQUEEZE_RIGHT , mythid)
895            write(msgbuf,'(A)')
896         &  'PART 4 in cost_gencost_sshv4 is skipped entirely, and there'
897          call print_message( msgbuf, errormessageunit,
898         &                    SQUEEZE_RIGHT , mythid)
899            write(msgbuf,'(A)')
900         &  'will be NO report of [tp,ers,gfo,lsc,gmsl] in costfunction.'
901          call print_message( msgbuf, errormessageunit,
902         &                    SQUEEZE_RIGHT , mythid)
903          endif
904    catn -----------
905    
906        do irec = 1, ndaysrec-ndaysave+1        do irec = 1, ndaysrec-ndaysave+1
907    
908         do bj = jtlo,jthi         do bj = jtlo,jthi
# Line 780  c -------------------------------------- Line 927  c --------------------------------------
927    
928          krec=irec+jrec-1          krec=irec+jrec-1
929    
930          call active_read_xy( fname, psbar, krec, doglobalread,  #ifdef ALLOW_AUTODIFF
931       &                       ladinit, optimcycle, mythid,          call active_read_xy( fname, etaday, krec, doglobalread,
932       &                       xx_psbar_mean_dummy )       &                       ladinit, eccoiter, mythid,
933         &                       gencost_dummy(igen_etaday) )
934  #ifndef ALLOW_PSBAR_MEAN  #else
935          CALL REMOVE_MEAN_RL( 1, psbar, maskInC, maskInC, rA, drF,          CALL READ_REC_XY_RL( fname, etaday, kRec, 1, myThid )
      &        'psbar', myTime, myThid )  
936  #endif  #endif
937    
938  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION          if (.NOT.useEtaMean)
939        call cost_sla_read( topexfile, topexstartdate, topexperiod,       &  CALL REMOVE_MEAN_RL( 1, etaday, maskInC, maskInC, rA, drF,
940       &                topexintercept, topexslope,       &        'etaday', 0. _d 0, myThid )
941       &                tpobs, tpmask,  
942           do bj = jtlo,jthi
943            do bi = itlo,ithi
944             do j = jmin,jmax
945              do i = imin,imax
946                  tpma(i,j,bi,bj) = 0. _d 0
947                  tpob(i,j,bi,bj) = 0. _d 0
948                  ersma(i,j,bi,bj) = 0. _d 0
949                  ersob(i,j,bi,bj) = 0. _d 0
950                  gfoma(i,j,bi,bj) = 0. _d 0
951                  gfoob(i,j,bi,bj) = 0. _d 0
952              enddo
953             enddo
954            enddo
955           enddo
956    
957            if (using_tpj)
958         &  call cost_sla_read( tpfi, tpsta, tpper,
959         &                zeroRL, zeroRL,
960         &                tpob, tpma,
961       &                krec, mythid )       &                krec, mythid )
962  #endif          if (using_ers)
963  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION       &  call cost_sla_read( ersfi, erssta, ersper,
964        call cost_sla_read( ersfile, ersstartdate, ersperiod,       &                zeroRL, zeroRL,
965       &                ersintercept, ersslope,       &                ersob, ersma,
      &                ersobs, ersmask,  
966       &                krec, mythid )       &                krec, mythid )
967  #endif          if (using_gfo)
968  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION       &  call cost_sla_read( gfofi, gfosta, gfoper,
969        call cost_sla_read( gfofile, gfostartdate, gfoperiod,       &                zeroRL, zeroRL,
970       &                gfointercept, gfoslope,       &                gfoob, gfoma,
      &                gfoobs, gfomask,  
971       &                krec, mythid )       &                krec, mythid )
 #endif  
972    
973         do bj = jtlo,jthi         do bj = jtlo,jthi
974          do bi = itlo,ithi          do bi = itlo,ithi
975           do j = jmin,jmax           do j = jmin,jmax
976            do i = imin,imax            do i = imin,imax
977  #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_psMslaobs_MSK(i,j,bi,bj)  
978       &  .NE.0. ) then       &  .NE.0. ) then
979             anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+             anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
980       &        psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)       &        etaday(i,j,bi,bj)-tpob(i,j,bi,bj)
981       &        -mean_psMslaobs(i,j,bi,bj)       &        -mean_psMslaobs(i,j,bi,bj)
982             anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+             anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
983       &        tpobs(i,j,bi,bj)       &        tpob(i,j,bi,bj)
984             anom_psMslaobs_NUM(i,j,bi,bj)=             anom_psMslaobs_NUM(i,j,bi,bj)=
985       &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0       &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
986             if ( (2*jrec.EQ.ndaysave).OR.(2*jrec+1.EQ.ndaysave) )             if ( (2*jrec.EQ.ndaysave).OR.(2*jrec+1.EQ.ndaysave) )
987       &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0       &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0
988        endif        endif
989  #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_psMslaobs_MSK(i,j,bi,bj)  
990       &  .NE.0. ) then       &  .NE.0. ) then
991             anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+             anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
992       &        psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)       &        etaday(i,j,bi,bj)-ersob(i,j,bi,bj)
993       &        -mean_psMslaobs(i,j,bi,bj)       &        -mean_psMslaobs(i,j,bi,bj)
994             anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+             anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
995       &        ersobs(i,j,bi,bj)       &        ersob(i,j,bi,bj)
996             anom_psMslaobs_NUM(i,j,bi,bj)=             anom_psMslaobs_NUM(i,j,bi,bj)=
997       &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0       &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
998             if ( (2*jrec.EQ.ndaysave).OR.(2*jrec+1.EQ.ndaysave) )             if ( (2*jrec.EQ.ndaysave).OR.(2*jrec+1.EQ.ndaysave) )
999       &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0       &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0
1000        endif        endif
1001  #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_psMslaobs_MSK(i,j,bi,bj)  
1002       &  .NE.0. ) then       &  .NE.0. ) then
1003             anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+             anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
1004       &        psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)       &        etaday(i,j,bi,bj)-gfoob(i,j,bi,bj)
1005       &        -mean_psMslaobs(i,j,bi,bj)       &        -mean_psMslaobs(i,j,bi,bj)
1006             anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+             anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
1007       &        gfoobs(i,j,bi,bj)       &        gfoob(i,j,bi,bj)
1008             anom_psMslaobs_NUM(i,j,bi,bj)=             anom_psMslaobs_NUM(i,j,bi,bj)=
1009       &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0       &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
1010             if ( (2*jrec.EQ.ndaysave).OR.(2*jrec+1.EQ.ndaysave) )             if ( (2*jrec.EQ.ndaysave).OR.(2*jrec+1.EQ.ndaysave) )
1011       &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0       &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0
1012        endif        endif
 #endif  
1013            enddo            enddo
1014           enddo           enddo
1015          enddo          enddo
# Line 867  c -------------------------------------- Line 1023  c --------------------------------------
1023                do i = imin,imax                do i = imin,imax
1024                 if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.                 if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
1025       &              ( maskc(i,j,1,bi,bj) .NE. 0. )       &              ( maskc(i,j,1,bi,bj) .NE. 0. )
 #ifndef ALLOW_HIGHLAT_ALTIMETRY      
      &              .AND.( abs(YC(i,j,bi,bj)) .LE. 66. )  
 #endif  
1026       &            ) then       &            ) then
1027                    anom_psMslaobs(i,j,bi,bj) =                    anom_psMslaobs(i,j,bi,bj) =
1028       &                 anom_psMslaobs(i,j,bi,bj) /       &                 anom_psMslaobs(i,j,bi,bj) /
# Line 899  c -------------------------------------- Line 1052  c --------------------------------------
1052                do i = imin,imax                do i = imin,imax
1053                 if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.                 if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
1054       &              ( maskc(i,j,1,bi,bj) .NE. 0. )       &              ( maskc(i,j,1,bi,bj) .NE. 0. )
 #ifndef ALLOW_HIGHLAT_ALTIMETRY      
      &              .AND.( abs(YC(i,j,bi,bj)) .LE. 66. )  
 #endif  
1055       &            ) then       &            ) then
1056                    slaoffset=slaoffset+RA(i,j,bi,bj)*                    slaoffset=slaoffset+RA(i,j,bi,bj)*
1057       &                      anom_psMslaobs(i,j,bi,bj)       &                      anom_psMslaobs(i,j,bi,bj)
# Line 921  c Since slaoffset is itself area weighte Line 1071  c Since slaoffset is itself area weighte
1071        _GLOBAL_SUM_RL( slaoffset_sum , mythid )        _GLOBAL_SUM_RL( slaoffset_sum , mythid )
1072    
1073        if (slaoffset_sum .eq. 0.0) then        if (slaoffset_sum .eq. 0.0) then
1074            if (gencost_outputlevel(igen_gmsl).GT.0) then
1075          _BEGIN_MASTER( mythid )          _BEGIN_MASTER( mythid )
1076          write(msgbuf,'(a)') ' sshv4: slaoffset_sum = zero!'          write(msgbuf,'(a)') ' sshv4: slaoffset_sum = zero!'
1077          call print_message( msgbuf, standardmessageunit,          call print_message( msgbuf, standardmessageunit,
1078       &                          SQUEEZE_RIGHT , mythid)       &                          SQUEEZE_RIGHT , mythid)
1079          _END_MASTER( mythid )          _END_MASTER( mythid )
1080            endif
1081  c        stop   '  ... stopped in cost_ssh.'  c        stop   '  ... stopped in cost_ssh.'
1082        else        else
1083            if (gencost_outputlevel(igen_gmsl).GT.0) then
1084          _BEGIN_MASTER( mythid )          _BEGIN_MASTER( mythid )
1085          write(msgbuf,'(a,3d22.15)') ' sshv4:slaoffset,sum,weight=',          write(msgbuf,'(a,3d22.15)') ' sshv4:slaoffset,sum,weight=',
1086       &         slaoffset,slaoffset_sum,slaoffset_weight       &         slaoffset,slaoffset_sum,slaoffset_weight
1087          call print_message( msgbuf, standardmessageunit,          call print_message( msgbuf, standardmessageunit,
1088       &                          SQUEEZE_RIGHT , mythid)       &                          SQUEEZE_RIGHT , mythid)
1089          _END_MASTER( mythid )          _END_MASTER( mythid )
1090            endif
1091        endif        endif
1092    
1093  c--   Compute (average) slaoffset  c--   Compute (average) slaoffset
# Line 946  c--   Subtract slaoffset from anom_psMsl Line 1100  c--   Subtract slaoffset from anom_psMsl
1100                do i = imin,imax                do i = imin,imax
1101                 if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.                 if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
1102       &              ( maskc(i,j,1,bi,bj) .NE. 0. )       &              ( maskc(i,j,1,bi,bj) .NE. 0. )
 #ifndef ALLOW_HIGHLAT_ALTIMETRY      
      &              .AND.( abs(YC(i,j,bi,bj)) .LE. 66. )  
 #endif  
1103       &            ) then       &            ) then
1104                    anom_psMslaobs(i,j,bi,bj) =                    anom_psMslaobs(i,j,bi,bj) =
1105       &                 anom_psMslaobs(i,j,bi,bj) - slaoffset       &                 anom_psMslaobs(i,j,bi,bj) - slaoffset
# Line 966  c--   Subtract slaoffset from anom_psMsl Line 1117  c--   Subtract slaoffset from anom_psMsl
1117  c PART 4.2: smooth anom_psMslaobs in space  c PART 4.2: smooth anom_psMslaobs in space
1118  c ----------------------------------------  c ----------------------------------------
1119    
1120  #ifdef ALLOW_GENCOST_SSHV4_OUTPUT        if (gencost_outputlevel(igen_lsc).GT.0) then
1121        write(fname4test(1:80),'(1a)') 'sladiff_raw'        write(fname4test(1:80),'(1a)') 'sladiff_raw'
1122        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
1123       & 1,anom_psMslaobs,irec,1,mythid)       & 1,anom_psMslaobs,irec,1,mythid)
# Line 974  c -------------------------------------- Line 1125  c --------------------------------------
1125        write(fname4test(1:80),'(1a)') 'slaobs_raw'        write(fname4test(1:80),'(1a)') 'slaobs_raw'
1126        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
1127       & 1,anom_slaobs,irec,1,mythid)       & 1,anom_slaobs,irec,1,mythid)
1128  #endif        endif
1129    
1130  #ifdef ALLOW_SMOOTH  #ifdef ALLOW_SMOOTH
1131        if ( useSMOOTH )        if ( useSMOOTH )
# Line 983  c -------------------------------------- Line 1134  c --------------------------------------
1134       &     gencost_smooth2Ddiffnbt(igen_lsc),mythid)       &     gencost_smooth2Ddiffnbt(igen_lsc),mythid)
1135  #endif  #endif
1136    
1137  #ifdef ALLOW_GENCOST_SSHV4_OUTPUT        if (gencost_outputlevel(igen_lsc).GT.0) then
1138  #ifdef ALLOW_SMOOTH  #ifdef ALLOW_SMOOTH
1139        if ( useSMOOTH )        if ( useSMOOTH )
1140       &  call smooth_hetero2d(anom_slaobs,maskc,       &  call smooth_hetero2d(anom_slaobs,maskc,
1141       &     gencost_scalefile(igen_lsc),       &     gencost_scalefile(igen_lsc),
1142       &     gencost_smooth2Ddiffnbt(igen_lsc),mythid)       &     gencost_smooth2Ddiffnbt(igen_lsc),mythid)
1143  #endif  #endif
1144    c
1145        write(fname4test(1:80),'(1a)') 'sladiff_smooth'        write(fname4test(1:80),'(1a)') 'sladiff_smooth'
1146        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
1147       & 1,anom_psMslaobs,irec,1,mythid)       & 1,anom_psMslaobs,irec,1,mythid)
1148    c
1149        write(fname4test(1:80),'(1a)') 'slaobs_smooth'        write(fname4test(1:80),'(1a)') 'slaobs_smooth'
1150        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
1151       & 1,anom_slaobs,irec,1,mythid)       & 1,anom_slaobs,irec,1,mythid)
1152  #endif        endif
1153    
1154  c PART 4.3: compute cost function term  c PART 4.3: compute cost function term
1155  c ------------------------------------  c ------------------------------------
# Line 1007  c ------------------------------------ Line 1158  c ------------------------------------
1158          do bi = itlo,ithi          do bi = itlo,ithi
1159           do j = jmin,jmax           do j = jmin,jmax
1160            do i = imin,imax            do i = imin,imax
 # if (defined (ALLOW_SSH_GFOANOM_COST_CONTRIBUTION) || \  
       defined (ALLOW_SSH_TPANOM_COST_CONTRIBUTION) || \  
       defined (ALLOW_SSH_ERSANOM_COST_CONTRIBUTION))  
1161            if ( (gencost_weight(i,j,bi,bj,igen_lsc).GT.0.).AND.            if ( (gencost_weight(i,j,bi,bj,igen_lsc).GT.0.).AND.
1162       &         (anom_psMslaobs_MSK(i,j,bi,bj).GT.0.) ) then       &         (anom_psMslaobs_MSK(i,j,bi,bj).GT.0.) ) then
1163              junk = anom_psMslaobs(i,j,bi,bj)              junk = anom_psMslaobs(i,j,bi,bj)
# Line 1026  c ------------------------------------ Line 1174  c ------------------------------------
1174       &         (anom_psMslaobs_MSK(i,j,bi,bj).GT.0.) )       &         (anom_psMslaobs_MSK(i,j,bi,bj).GT.0.) )
1175       &         num_gencost(bi,bj,igen_lsc) =       &         num_gencost(bi,bj,igen_lsc) =
1176       &         num_gencost(bi,bj,igen_lsc) + 1. _d 0       &         num_gencost(bi,bj,igen_lsc) + 1. _d 0
 #endif  
1177             enddo             enddo
1178           enddo           enddo
1179          enddo          enddo
1180         enddo         enddo
1181    
1182  #ifdef ALLOW_GENCOST_SSHV4_OUTPUT        if (gencost_outputlevel(igen_lsc).GT.0) then
1183        write(fname4test(1:80),'(1a)') 'sladiff_subsample'        write(fname4test(1:80),'(1a)') 'sladiff_subsample'
1184        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
1185       & 1,anom_psMslaobs_SUB,irec,1,mythid)       & 1,anom_psMslaobs_SUB,irec,1,mythid)
1186    c
1187        write(fname4test(1:80),'(1a)') 'slaobs_subsample'        write(fname4test(1:80),'(1a)') 'slaobs_subsample'
1188        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
1189       & 1,anom_slaobs_SUB,irec,1,mythid)       & 1,anom_slaobs_SUB,irec,1,mythid)
1190  #endif        endif
1191    
1192  c PART 4.4: compute cost function term for global mean sea level  c PART 4.4: compute cost function term for global mean sea level
1193  c --------------------------------------------------------------  c --------------------------------------------------------------
1194    
1195        IF ( MASTER_CPU_THREAD(myThid).AND.(igen_gmsl.NE.0) ) THEN        IF ( ( MASTER_CPU_THREAD(myThid).AND.(igen_gmsl.NE.0) ).AND.
1196         &     ( slaoffset_sum.GT.0.25 _d 0 * globalArea ) ) THEN
1197            junk=slaoffset            junk=slaoffset
1198            junkweight=slaoffset_weight            junkweight=1. _d 0
1199            objf_gencost(1,1,igen_gmsl) = objf_gencost(1,1,igen_gmsl)            objf_gencost(1,1,igen_gmsl) = objf_gencost(1,1,igen_gmsl)
1200       &        + junk*junk*junkweight       &        + junk*junk*junkweight
1201            num_gencost(1,1,igen_gmsl) = num_gencost(1,1,igen_gmsl)            num_gencost(1,1,igen_gmsl) = num_gencost(1,1,igen_gmsl)
# Line 1062  cgf ==================================== Line 1210  cgf ====================================
1210  c time associated with this ndaysrec mean  c time associated with this ndaysrec mean
1211          krec = irec + (ndaysave/2)          krec = irec + (ndaysave/2)
1212    
1213          call active_read_xy( fname, psbar, krec, doglobalread,  #ifdef ALLOW_AUTODIFF
1214       &                       ladinit, optimcycle, mythid,          call active_read_xy( fname, etaday, krec, doglobalread,
1215       &                       xx_psbar_mean_dummy )       &                       ladinit, eccoiter, mythid,
1216         &                       gencost_dummy(igen_etaday) )
1217  #ifndef ALLOW_PSBAR_MEAN  #else
1218          CALL REMOVE_MEAN_RL( 1, psbar, maskInC, maskInC, rA, drF,          CALL READ_REC_XY_RL( fname, etaday, kRec, 1, myThid )
      &        'psbar', myTime, myThid )  
1219  #endif  #endif
1220    
1221  #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION          if (.NOT.useEtaMean)
1222          call cost_readtopex( krec, mythid )       &  CALL REMOVE_MEAN_RL( 1, etaday, maskInC, maskInC, rA, drF,
1223  #endif       &        'etaday', 0. _d 0, myThid )
1224  #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION  
1225          call cost_readers( krec, mythid )         do bj = jtlo,jthi
1226  #endif          do bi = itlo,ithi
1227  #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION           do j = jmin,jmax
1228          call cost_readgfo( krec, mythid )            do i = imin,imax
1229  #endif                tpma(i,j,bi,bj) = 0. _d 0
1230                  tpob(i,j,bi,bj) = 0. _d 0
1231                  ersma(i,j,bi,bj) = 0. _d 0
1232                  ersob(i,j,bi,bj) = 0. _d 0
1233                  gfoma(i,j,bi,bj) = 0. _d 0
1234                  gfoob(i,j,bi,bj) = 0. _d 0
1235              enddo
1236             enddo
1237            enddo
1238           enddo
1239    
1240            if (using_tpj)
1241         &  call cost_sla_read( tpfi, tpsta, tpper,
1242         &                zeroRL, zeroRL,
1243         &                tpob, tpma,
1244         &                krec, mythid )
1245            if (using_ers)
1246         &  call cost_sla_read( ersfi, erssta, ersper,
1247         &                zeroRL, zeroRL,
1248         &                ersob, ersma,
1249         &                krec, mythid )
1250            if (using_gfo)
1251         &  call cost_sla_read( gfofi, gfosta, gfoper,
1252         &                zeroRL, zeroRL,
1253         &                gfoob, gfoma,
1254         &                krec, mythid )
1255    
1256         do bj = jtlo,jthi         do bj = jtlo,jthi
1257          do bi = itlo,ithi          do bi = itlo,ithi
# Line 1100  c time associated with this ndaysrec mea Line 1272  c time associated with this ndaysrec mea
1272          do bi = itlo,ithi          do bi = itlo,ithi
1273           do j = jmin,jmax           do j = jmin,jmax
1274            do i = imin,imax            do i = imin,imax
1275  #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_psMslaobs_MSK(i,j,bi,bj).NE.0. )  
1276       & then       & then
1277           anom_psMtpobs(i,j,bi,bj)=           anom_psMtpobs(i,j,bi,bj)=
1278       &       psbar(i,j,bi,bj) - tpobs(i,j,bi,bj)       &       etaday(i,j,bi,bj) - tpob(i,j,bi,bj)
1279       &       - mean_psMslaobs(i,j,bi,bj) - slaoffset       &       - mean_psMslaobs(i,j,bi,bj) - slaoffset
1280       &       - anom_psMslaobs(i,j,bi,bj)       &       - anom_psMslaobs(i,j,bi,bj)
1281           anom_tpobs(i,j,bi,bj)=tpobs(i,j,bi,bj) - slaoffset           anom_tpobs(i,j,bi,bj)=tpob(i,j,bi,bj) - slaoffset
1282       &       - anom_slaobs(i,j,bi,bj)       &       - anom_slaobs(i,j,bi,bj)
1283        endif        endif
1284  #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_psMslaobs_MSK(i,j,bi,bj).NE.0. )  
1285       & then       & then
1286           anom_psMersobs(i,j,bi,bj)=           anom_psMersobs(i,j,bi,bj)=
1287       &       psbar(i,j,bi,bj) - ersobs(i,j,bi,bj)       &       etaday(i,j,bi,bj) - ersob(i,j,bi,bj)
1288       &       - mean_psMslaobs(i,j,bi,bj) - slaoffset       &       - mean_psMslaobs(i,j,bi,bj) - slaoffset
1289       &       - anom_psMslaobs(i,j,bi,bj)       &       - anom_psMslaobs(i,j,bi,bj)
1290           anom_ersobs(i,j,bi,bj)=ersobs(i,j,bi,bj) - slaoffset           anom_ersobs(i,j,bi,bj)=ersob(i,j,bi,bj) - slaoffset
1291       &       - anom_slaobs(i,j,bi,bj)       &       - anom_slaobs(i,j,bi,bj)
1292        endif        endif
1293  #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_psMslaobs_MSK(i,j,bi,bj).NE.0. )  
1294       & then       & then
1295           anom_psMgfoobs(i,j,bi,bj)=           anom_psMgfoobs(i,j,bi,bj)=
1296       &       psbar(i,j,bi,bj) - gfoobs(i,j,bi,bj)       &       etaday(i,j,bi,bj) - gfoob(i,j,bi,bj)
1297       &       - mean_psMslaobs(i,j,bi,bj) - slaoffset       &       - mean_psMslaobs(i,j,bi,bj) - slaoffset
1298       &       - anom_psMslaobs(i,j,bi,bj)       &       - anom_psMslaobs(i,j,bi,bj)
1299           anom_gfoobs(i,j,bi,bj)=gfoobs(i,j,bi,bj) - slaoffset           anom_gfoobs(i,j,bi,bj)=gfoob(i,j,bi,bj) - slaoffset
1300       &       - anom_slaobs(i,j,bi,bj)       &       - anom_slaobs(i,j,bi,bj)
1301        endif        endif
 #endif  
1302            enddo            enddo
1303           enddo           enddo
1304          enddo          enddo
1305         enddo         enddo
1306    
1307  #ifdef ALLOW_GENCOST_SSHV4_OUTPUT        if (gencost_outputlevel(igen_tp).GT.0) then
1308        write(fname4test(1:80),'(1a)') 'sladiff_tp_raw'        write(fname4test(1:80),'(1a)') 'sladiff_tp_raw'
1309        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
1310       & 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)  
   
1311        write(fname4test(1:80),'(1a)') 'slaobs_tp_raw'        write(fname4test(1:80),'(1a)') 'slaobs_tp_raw'
1312        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
1313       & 1,anom_tpobs,irec,1,mythid)       & 1,anom_tpobs,irec,1,mythid)
1314          endif
1315    
1316          if (gencost_outputlevel(igen_ers).GT.0) then
1317          write(fname4test(1:80),'(1a)') 'sladiff_ers_raw'
1318          call mdswritefield(fname4test,32,.false.,'RL',
1319         & 1,anom_psMersobs,irec,1,mythid)
1320        write(fname4test(1:80),'(1a)') 'slaobs_ers_raw'        write(fname4test(1:80),'(1a)') 'slaobs_ers_raw'
1321        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
1322       & 1,anom_ersobs,irec,1,mythid)       & 1,anom_ersobs,irec,1,mythid)
1323          endif
1324    
1325          if (gencost_outputlevel(igen_gfo).GT.0) then
1326          write(fname4test(1:80),'(1a)') 'sladiff_gfo_raw'
1327          call mdswritefield(fname4test,32,.false.,'RL',
1328         & 1,anom_psMgfoobs,irec,1,mythid)
1329        write(fname4test(1:80),'(1a)') 'slaobs_gfo_raw'        write(fname4test(1:80),'(1a)') 'slaobs_gfo_raw'
1330        call mdswritefield(fname4test,32,.false.,'RL',        call mdswritefield(fname4test,32,.false.,'RL',
1331       & 1,anom_gfoobs,irec,1,mythid)       & 1,anom_gfoobs,irec,1,mythid)
1332  #endif        endif
1333    
1334         do bj = jtlo,jthi         do bj = jtlo,jthi
1335          do bi = itlo,ithi          do bi = itlo,ithi
1336           do j = jmin,jmax           do j = jmin,jmax
1337            do i = imin,imax            do i = imin,imax
 #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION  
1338  c--             The array psobs contains SSH anomalies.  c--             The array psobs contains SSH anomalies.
1339                  junkweight = mean_psMslaobs_MSK(i,j,bi,bj)                  junkweight = mean_psMslaobs_MSK(i,j,bi,bj)
1340       &                      *gencost_weight(i,j,bi,bj,igen_tp)       &                      *gencost_weight(i,j,bi,bj,igen_tp)
1341       &                      *tpmask(i,j,bi,bj)       &                      *tpma(i,j,bi,bj)
      &                      *cosphi(i,j,bi,bj)  
1342                  junk = anom_psMtpobs(i,j,bi,bj)                  junk = anom_psMtpobs(i,j,bi,bj)
1343                  objf_gencost(bi,bj,igen_tp) =                  objf_gencost(bi,bj,igen_tp) =
1344       &              objf_gencost(bi,bj,igen_tp)+junk*junk*junkweight       &              objf_gencost(bi,bj,igen_tp)+junk*junk*junkweight
1345                  if ( junkweight .ne. 0. )                  if ( junkweight .ne. 0. )
1346       &              num_gencost(bi,bj,igen_tp) =       &              num_gencost(bi,bj,igen_tp) =
1347       &              num_gencost(bi,bj,igen_tp) + 1. _d 0       &              num_gencost(bi,bj,igen_tp) + 1. _d 0
 #endif  
 #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION  
1348  c--             The array ersobs contains SSH anomalies.  c--             The array ersobs contains SSH anomalies.
1349                  junkweight = mean_psMslaobs_MSK(i,j,bi,bj)                  junkweight = mean_psMslaobs_MSK(i,j,bi,bj)
1350       &                      *gencost_weight(i,j,bi,bj,igen_ers)       &                      *gencost_weight(i,j,bi,bj,igen_ers)
1351       &                      *ersmask(i,j,bi,bj)       &                      *ersma(i,j,bi,bj)
      &                      *cosphi(i,j,bi,bj)  
1352                  junk = anom_psMersobs(i,j,bi,bj)                  junk = anom_psMersobs(i,j,bi,bj)
1353                  objf_gencost(bi,bj,igen_ers) =                  objf_gencost(bi,bj,igen_ers) =
1354       &               objf_gencost(bi,bj,igen_ers)+junk*junk*junkweight       &               objf_gencost(bi,bj,igen_ers)+junk*junk*junkweight
1355                  if ( junkweight .ne. 0. )                  if ( junkweight .ne. 0. )
1356       &              num_gencost(bi,bj,igen_ers) =       &              num_gencost(bi,bj,igen_ers) =
1357       &              num_gencost(bi,bj,igen_ers) + 1. _d 0       &              num_gencost(bi,bj,igen_ers) + 1. _d 0
 #endif  
 #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION  
1358  c--             The array gfoobs contains SSH anomalies.  c--             The array gfoobs contains SSH anomalies.
1359                  junkweight = mean_psMslaobs_MSK(i,j,bi,bj)                  junkweight = mean_psMslaobs_MSK(i,j,bi,bj)
1360       &                      *gencost_weight(i,j,bi,bj,igen_gfo)       &                      *gencost_weight(i,j,bi,bj,igen_gfo)
1361       &                      *gfomask(i,j,bi,bj)       &                      *gfoma(i,j,bi,bj)
      &                      *cosphi(i,j,bi,bj)  
1362                  junk = anom_psMgfoobs(i,j,bi,bj)                  junk = anom_psMgfoobs(i,j,bi,bj)
1363                  objf_gencost(bi,bj,igen_gfo) =                  objf_gencost(bi,bj,igen_gfo) =
1364       &              objf_gencost(bi,bj,igen_gfo)+junk*junk*junkweight       &              objf_gencost(bi,bj,igen_gfo)+junk*junk*junkweight
1365                  if ( junkweight .ne. 0. )                  if ( junkweight .ne. 0. )
1366       &              num_gencost(bi,bj,igen_gfo) =       &              num_gencost(bi,bj,igen_gfo) =
1367       &              num_gencost(bi,bj,igen_gfo) + 1. _d 0       &              num_gencost(bi,bj,igen_gfo) + 1. _d 0
 #endif  
1368             enddo             enddo
1369           enddo           enddo
1370          enddo          enddo
# Line 1210  c--             The array gfoobs contain Line 1372  c--             The array gfoobs contain
1372    
1373        enddo        enddo
1374    
   
1375  #endif /* ifdef ALLOW_GENCOST_CONTRIBUTION */  #endif /* ifdef ALLOW_GENCOST_CONTRIBUTION */
 #endif /* ifdef ALLOW_SSH_COST_CONTRIBUTION */  
1376    
1377        end        end
Legend:
Removed from v.1.22  
changed lines
  Added in v.1.34
  ViewVC Help
Powered by ViewVC 1.1.22