Back to home page

MITgcm
 
 

    

File indexing completed on 2023-07-16 05:10:05 UTC

view on githubraw file Latest commit 025a9bb1 on 2023-07-15 15:12:22 UTC
8f7d13d0c9 Jean*0001 #include "ECCO_OPTIONS.h"
f1e8c1d7ee Gael*0002 
                0003       subroutine cost_gencost_sshv4(
0a1233e81a Ou W*0004      I                        myThid )
f1e8c1d7ee Gael*0005 
                0006 c     ==================================================================
                0007 c     SUBROUTINE cost_gencost_sshv4
                0008 c     ==================================================================
                0009 c
                0010 c     o Evaluate cost function contributions of sea surface height.
                0011 c
                0012 c        started: Gael Forget, Oct-2009
                0013 c
11c3150c71 Mart*0014 c        working assumption for the time mean dynamic topography (MDT)
                0015 c        constraint: the various SLA data sets (tp, ers, gfo) have been
                0016 c        consistenty cross-referenced, as done in the RADS data sets. We
                0017 c        do not need to know the reference dynamic topography (or
                0018 c        SSH/Geoid). Simply it is assumed that the biases between
                0019 c        instruments have been taken care of. This is only a matter for
                0020 c        the MDT constraint, not for SLA constraints (see below).
f1e8c1d7ee Gael*0021 c
                0022 cgf 1) there are a few hardcoded numbers that will eventually be put in common
                0023 cgf     blocks/namelists
                0024 cgf 2) there are a several refinements that should be considered, such as
                0025 cgf     modulating weights with respect to numbers of samples
                0026 c
                0027 c     ==================================================================
                0028 c     SUBROUTINE cost_gencost_sshv4
                0029 c     ==================================================================
                0030 
                0031       implicit none
                0032 
                0033 c     == global variables ==
                0034 
                0035 #include "EEPARAMS.h"
                0036 #include "SIZE.h"
                0037 #include "PARAMS.h"
                0038 #include "GRID.h"
f09238ab8f Gael*0039 #ifdef ALLOW_CAL
                0040 # include "cal.h"
                0041 #endif
                0042 #ifdef ALLOW_ECCO
13d362b8c1 Ou W*0043 # include "ECCO_SIZE.h"
                0044 # include "ECCO.h"
f09238ab8f Gael*0045 #endif
ebac42f582 An T*0046 #ifdef ALLOW_SMOOTH
68ab1d598d Gael*0047 # include "SMOOTH.h"
ebac42f582 An T*0048 #endif
f1e8c1d7ee Gael*0049 
                0050 c     == routine arguments ==
f8e779c983 antn*0051       integer myThid
f1e8c1d7ee Gael*0052 
                0053 #ifdef ALLOW_GENCOST_CONTRIBUTION
bc9bce397d Gael*0054 
                0055 c     == functions ==
                0056       LOGICAL  MASTER_CPU_THREAD
                0057       EXTERNAL MASTER_CPU_THREAD
0a1233e81a Ou W*0058       INTEGER  ILNBLNK
                0059       EXTERNAL ILNBLNK
bc9bce397d Gael*0060 
f1e8c1d7ee Gael*0061 c     == local variables ==
                0062 
                0063       integer bi,bj
                0064       integer i,j,k
                0065       integer itlo,ithi
                0066       integer jtlo,jthi
                0067       integer jmin,jmax
                0068       integer imin,imax
                0069       integer irec,jrec,krec
                0070       integer ilps
025a9bb173 antn*0071       integer nrec_year, nrec_total, nrec_mdt_ref
                0072 c introduce use_mon to switch from using daily to monthly data if true
                0073       logical use_mon
f1e8c1d7ee Gael*0074 
                0075       logical doglobalread
                0076       logical ladinit
                0077 
                0078 c mapping to gencost
bc9bce397d Gael*0079       integer igen_mdt, igen_lsc, igen_gmsl
8f7d13d0c9 Jean*0080       integer igen_tp, igen_ers, igen_gfo
025a9bb173 antn*0081       integer igen_etagcm
f1e8c1d7ee Gael*0082 
f09238ab8f Gael*0083       _RL spval, factor
21dbf3b3cd Gael*0084       _RL offset
f1e8c1d7ee Gael*0085       _RL offset_sum
21dbf3b3cd Gael*0086       _RL slaoffset
                0087       _RL slaoffset_sum
bc9bce397d Gael*0088       _RL slaoffset_weight
21dbf3b3cd Gael*0089 
6b881d9117 Gael*0090 c local arrays
f8e779c983 antn*0091       _RL num0array ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, nSx, nSy )
6c0b33a90f An T*0092       _RL num0total
                0093       integer tempinteger
                0094 
f8e779c983 antn*0095       _RL diagnosfld ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, nSx, nSy )
f09238ab8f Gael*0096 
f8e779c983 antn*0097       _RL mdtob(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0098       _RL tpob(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0099       _RL ersob(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0100       _RL gfoob(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0101       _RL mdtma(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0102       _RL tpma(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0103       _RL ersma(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0104       _RL gfoma(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
025a9bb173 antn*0105       _RL etagcm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f09238ab8f Gael*0106 
                0107       character*(MAX_LEN_FNAM) mdtfi, tpfi, ersfi, gfofi
                0108       integer tpsta(4), erssta(4), gfosta(4)
                0109       integer mdtsta(4), mdtend(4)
                0110       _RL tpper, ersper, gfoper
f1e8c1d7ee Gael*0111 
d39bcf1629 An T*0112 c for PART 1: re-reference MDT (mdt) to the inferred SLA reference field
f8e779c983 antn*0113       _RL psmean    ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, nSx, nSy )
                0114       _RL mean_slaobs_mdt(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0115       _RL mean_slaobs_NUM(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f1e8c1d7ee Gael*0116 
                0117 c for PART 2: compute time mean differences over the model period
f8e779c983 antn*0118       _RL mean_slaobs_model(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0119       _RL mean_psMssh_all(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0120       _RL mean_psMssh_all_NUM(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0121       _RL mean_psMssh_all_MSK(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f1e8c1d7ee Gael*0122 
f8e779c983 antn*0123       _RL mean_psMslaobs(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0124       _RL mean_psMslaobs_MSK(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
52e9fa3441 Gael*0125 
f8e779c983 antn*0126       _RL mean_psMtpobs(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0127       _RL mean_psMtpobs_NUM(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f1e8c1d7ee Gael*0128 
f8e779c983 antn*0129       _RL mean_psMersobs(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0130       _RL mean_psMersobs_NUM(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f1e8c1d7ee Gael*0131 
f8e779c983 antn*0132       _RL mean_psMgfoobs(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0133       _RL mean_psMgfoobs_NUM(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f1e8c1d7ee Gael*0134 
                0135 c for PART 4/5: compute smooth/raw anomalies
f8e779c983 antn*0136       _RL anom_psMslaobs(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0137       _RL anom_psMslaobs_NUM (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0138       _RL anom_psMslaobs_MSK (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0139       _RL anom_psMslaobs_SUB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0140       _RL anom_slaobs (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0141       _RL anom_slaobs_SUB  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f1e8c1d7ee Gael*0142 
f8e779c983 antn*0143       _RL anom_psMtpobs(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0144       _RL anom_tpobs(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f1e8c1d7ee Gael*0145 
f8e779c983 antn*0146       _RL anom_psMersobs(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0147       _RL anom_ersobs(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f1e8c1d7ee Gael*0148 
f8e779c983 antn*0149       _RL anom_psMgfoobs(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0150       _RL anom_gfoobs(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f1e8c1d7ee Gael*0151 
025a9bb173 antn*0152       integer mdt_y0,mdt_y1,year
f09238ab8f Gael*0153       integer num0
f1e8c1d7ee Gael*0154 
                0155       _RL junk,junkweight
                0156 
025a9bb173 antn*0157       integer nrecsave
f1e8c1d7ee Gael*0158 
3248fb8b5a Gael*0159       integer k2, k2_mdt, k2_lsc
b9b1c4d04c Timo*0160       logical dosumsq
3248fb8b5a Gael*0161 
de57a2ec4b Mart*0162       character*(MAX_LEN_FNAM) fname
                0163       character*(MAX_LEN_FNAM) fname4test
f1e8c1d7ee Gael*0164       character*(MAX_LEN_MBUF) msgbuf
                0165 
2a0641f1b7 Gael*0166       LOGICAL doReference
6b881d9117 Gael*0167       LOGICAL useEtaMean
2a0641f1b7 Gael*0168 
f1e8c1d7ee Gael*0169 c     == end of interface ==
                0170 
0a1233e81a Ou W*0171 C-    Initialise (again) gencost_weight
                0172       DO k = 1, NGENCOST
                0173        DO bj = myByLo(myThid), myByHi(myThid)
                0174         DO bi = myBxLo(myThid), myBxHi(myThid)
                0175          DO j=1-OLy,sNy+OLy
                0176           DO i=1-OLx,sNx+OLx
                0177            gencost_weight(i,j,bi,bj,k) = 0. _d 0
                0178           ENDDO
                0179          ENDDO
                0180         ENDDO
                0181        ENDDO
                0182       ENDDO
                0183 
f8e779c983 antn*0184       jtlo = myByLo(myThid)
                0185       jthi = myByHi(myThid)
                0186       itlo = myBxLo(myThid)
                0187       ithi = myBxHi(myThid)
f1e8c1d7ee Gael*0188       jmin = 1
f8e779c983 antn*0189       jmax = sNy
f1e8c1d7ee Gael*0190       imin = 1
f8e779c983 antn*0191       imax = sNx
f1e8c1d7ee Gael*0192 
                0193 c-- detect the relevant gencost indices
                0194       igen_mdt=0
                0195       igen_tp =0
                0196       igen_ers=0
                0197       igen_gfo=0
                0198       igen_lsc=0
bc9bce397d Gael*0199       igen_gmsl=0
025a9bb173 antn*0200       igen_etagcm=0
f1e8c1d7ee Gael*0201       do k=1,NGENCOST
                0202         if (gencost_name(k).EQ.'sshv4-mdt') igen_mdt=k
                0203         if (gencost_name(k).EQ.'sshv4-tp') igen_tp=k
                0204         if (gencost_name(k).EQ.'sshv4-ers') igen_ers=k
                0205         if (gencost_name(k).EQ.'sshv4-gfo') igen_gfo=k
                0206         if (gencost_name(k).EQ.'sshv4-lsc') igen_lsc=k
bc9bce397d Gael*0207         if (gencost_name(k).EQ.'sshv4-gmsl') igen_gmsl=k
f1e8c1d7ee Gael*0208       enddo
                0209 
3248fb8b5a Gael*0210       k2_mdt=0
                0211       k2_lsc=0
11c3150c71 Mart*0212       if (igen_mdt.GT.0) then
                0213        do k2 = 1, NGENPPROC
3248fb8b5a Gael*0214         if (gencost_posproc(k2,igen_mdt).EQ.'smooth') k2_mdt=k2
11c3150c71 Mart*0215        enddo
                0216       endif
                0217       if (igen_lsc.GT.0) then
                0218        do k2 = 1, NGENPPROC
3248fb8b5a Gael*0219         if (gencost_posproc(k2,igen_lsc).EQ.'smooth') k2_lsc=k2
11c3150c71 Mart*0220        enddo
                0221       endif
025a9bb173 antn*0222       write(msgbuf,'(a,2(1x,i5))')
                0223      &  'sshv4: k2_[mdt,lsc]: ', k2_mdt, k2_lsc
                0224       call print_message( msgbuf, standardmessageunit,
                0225      &                    SQUEEZE_RIGHT , myThid)
cbd85e4123 Gael*0226 
b9b1c4d04c Timo*0227       dosumsq = .TRUE.
11c3150c71 Mart*0228       if (igen_mdt.GT.0) then
0a1233e81a Ou W*0229        call ecco_zero( gencost_weight(1-OLx,1-OLy,1,1,igen_mdt),
                0230      &                 1, zeroRL, myThid )
11c3150c71 Mart*0231        if ( gencost_errfile(igen_mdt) .NE. ' ' )
                0232      &   call ecco_readwei( gencost_errfile(igen_mdt),
                0233      &     gencost_weight(1-OLx,1-OLy,1,1,igen_mdt),
                0234      &     1, 1, 1, dosumsq, myThid )
                0235       endif
                0236       if (igen_lsc.GT.0) then
0a1233e81a Ou W*0237        call ecco_zero( gencost_weight(1-OLx,1-OLy,1,1,igen_lsc),
                0238      &                 1, zeroRL, myThid )
11c3150c71 Mart*0239        if ( gencost_errfile(igen_lsc) .NE. ' ' )
                0240      &   call ecco_readwei( gencost_errfile(igen_lsc),
                0241      &     gencost_weight(1-OLx,1-OLy,1,1,igen_lsc),
                0242      &     1, 1, 1, dosumsq, myThid )
                0243       endif
                0244       if (igen_tp.GT.0) then
0a1233e81a Ou W*0245        call ecco_zero( gencost_weight(1-OLx,1-OLy,1,1,igen_tp),
                0246      &                 1, zeroRL, myThid )
11c3150c71 Mart*0247        if ( gencost_errfile(igen_tp) .NE. ' ' )
                0248      &   call ecco_readwei( gencost_errfile(igen_tp),
                0249      &     gencost_weight(1-OLx,1-OLy,1,1,igen_tp),
                0250      &     1, 1, 1, dosumsq, myThid )
                0251       endif
                0252       if (igen_ers.GT.0) then
0a1233e81a Ou W*0253        call ecco_zero( gencost_weight(1-OLx,1-OLy,1,1,igen_ers),
                0254      &                 1, zeroRL, myThid )
11c3150c71 Mart*0255        if ( gencost_errfile(igen_ers) .NE. ' ' )
                0256      &   call ecco_readwei( gencost_errfile(igen_ers),
                0257      &     gencost_weight(1-OLx,1-OLy,1,1,igen_ers),
                0258      &     1, 1, 1, dosumsq, myThid )
                0259       endif
                0260       if (igen_gfo.GT.0) then
0a1233e81a Ou W*0261        call ecco_zero( gencost_weight(1-OLx,1-OLy,1,1,igen_gfo),
                0262      &                 1, zeroRL, myThid )
11c3150c71 Mart*0263        if ( gencost_errfile(igen_gfo) .NE. ' ' )
                0264      &   call ecco_readwei( gencost_errfile(igen_gfo),
                0265      &     gencost_weight(1-OLx,1-OLy,1,1,igen_gfo),
                0266      &     1, 1, 1, dosumsq, myThid )
                0267       endif
6b881d9117 Gael*0268 
                0269 c switch for excluding global mean
                0270       useEtaMean=.TRUE.
                0271 
f8e779c983 antn*0272       write(msgbuf,'(a,l)')
6b881d9117 Gael*0273      & ' sshv4: use global mean of eta useEtaMean=',useEtaMean
                0274       call print_message( msgbuf, standardmessageunit,
f8e779c983 antn*0275      &                    SQUEEZE_RIGHT , myThid)
6b881d9117 Gael*0276 
025a9bb173 antn*0277 C Check to determine use_mon and nrec_year. Check all igen_[tp,ers,gfo]
                0278 C in case some of the data files are not defined.
                0279       use_mon = .FALSE.
                0280       IF ( igen_tp.GT.0 ) use_mon = use_mon .OR.
                0281      &     gencost_avgperiod(igen_tp) .EQ.'month' .OR.
                0282      &     gencost_avgperiod(igen_tp) .EQ.'MONTH'
                0283       IF ( igen_ers.GT.0 ) use_mon = use_mon .OR.
                0284      &     gencost_avgperiod(igen_ers) .EQ.'month' .OR.
                0285      &     gencost_avgperiod(igen_ers) .EQ.'MONTH'
                0286       IF ( igen_gfo.GT.0 ) use_mon = use_mon .OR.
                0287      &      gencost_avgperiod(igen_gfo).EQ.'month' .OR.
                0288      &      gencost_avgperiod(igen_gfo).EQ.'MONTH'
                0289 
6b881d9117 Gael*0290 c minimum number of observations to consider re-referenced MDT
025a9bb173 antn*0291       nrec_year=366
                0292       nrec_total=ndaysrec
                0293       IF (use_mon) THEN
                0294         num0=3
                0295         nrec_mdt_ref=12
                0296         nrec_year=12
                0297         nrec_total=nmonsrec
                0298       else
                0299         num0=100
                0300         nrec_mdt_ref=365
                0301       endif
                0302 
                0303       write(msgbuf,'(a,L5,2(1x,i10),2(1x,i5))')
                0304      & ' sshv4: use_mon, nrec_year, nrec_total, num0, nrec_mdt_ref: ',
                0305      &          use_mon, nrec_year, nrec_total, num0, nrec_mdt_ref
                0306       call print_message( msgbuf, standardmessageunit,
                0307      &                    SQUEEZE_RIGHT , myThid)
6b881d9117 Gael*0308 c determine whether or not to re-reference
6c0b33a90f An T*0309 c-- not only model period has to fall within the mdt period
025a9bb173 antn*0310 c-- but that there has to be at least 365 days (or 12 months) of model
                0311 c-- run so for short model run, this will always get set to FALSE
                0312 C Note also the hard-coded years 1992,2011 here that should be
                0313 c removed at some point in the future if the reference MDT start
                0314 c date no longer falls between 1992 and 2011
6b881d9117 Gael*0315       doReference=.FALSE.
025a9bb173 antn*0316       IF ((modelstartdate(1).GT.1992*10000).AND.
6b881d9117 Gael*0317      &    (modelstartdate(1).LT.2011*10000).AND.
025a9bb173 antn*0318      &    (nrec_total.GE.nrec_mdt_ref))  doReference=.TRUE.
6b881d9117 Gael*0319 
                0320       write(msgbuf,'(a,l)') ' sshv4:re-reference MDT=',doReference
                0321       call print_message( msgbuf, standardmessageunit,
f8e779c983 antn*0322      &                    SQUEEZE_RIGHT , myThid)
6b881d9117 Gael*0323 
f09238ab8f Gael*0324 c-- initialize local arrays
                0325       do bj = jtlo,jthi
                0326         do bi = itlo,ithi
                0327           do j = jmin,jmax
                0328             do i = imin,imax
                0329               mdtma(i,j,bi,bj) = 0. _d 0
                0330               mdtob(i,j,bi,bj) = 0. _d 0
                0331               tpma(i,j,bi,bj) = 0. _d 0
                0332               tpob(i,j,bi,bj) = 0. _d 0
                0333               ersma(i,j,bi,bj) = 0. _d 0
                0334               ersob(i,j,bi,bj) = 0. _d 0
                0335               gfoma(i,j,bi,bj) = 0. _d 0
                0336               gfoob(i,j,bi,bj) = 0. _d 0
                0337             enddo
                0338           enddo
                0339         enddo
                0340       enddo
                0341 
                0342 c mdtfi, mdtsta, mdtend
                0343       if (igen_mdt.GT.0) then
                0344        ilps=ilnblnk( gencost_datafile(igen_mdt) )
                0345        mdtfi=gencost_datafile(igen_mdt)(1:ilps)
f8e779c983 antn*0346        call cal_CopyDate(gencost_startdate(1,igen_mdt),mdtsta,myThid)
                0347        call cal_CopyDate(gencost_enddate(1,igen_mdt), mdtend, myThid)
11c3150c71 Mart*0348        if (gencost_barfile(igen_mdt)(1:5).EQ.'m_eta')
025a9bb173 antn*0349      &    igen_etagcm=igen_mdt
f09238ab8f Gael*0350       endif
                0351 c tpfi, tpsta, tpper
                0352       if (igen_tp.GT.0) then
                0353        ilps=ilnblnk( gencost_datafile(igen_tp) )
                0354        tpfi=gencost_datafile(igen_tp)(1:ilps)
f8e779c983 antn*0355        call cal_CopyDate(gencost_startdate(1,igen_tp),tpsta,myThid)
f09238ab8f Gael*0356        tpper=gencost_period(igen_tp)
f8e779c983 antn*0357        if (gencost_barfile(igen_tp)(1:5).EQ.'m_eta')
025a9bb173 antn*0358      &    igen_etagcm=igen_tp
f09238ab8f Gael*0359       endif
                0360 c ersfi, erssta, ersper
                0361       if (igen_ers.GT.0) then
                0362        ilps=ilnblnk( gencost_datafile(igen_ers) )
                0363        ersfi=gencost_datafile(igen_ers)(1:ilps)
f8e779c983 antn*0364        call cal_CopyDate(gencost_startdate(1,igen_ers),erssta,myThid)
f09238ab8f Gael*0365        ersper=gencost_period(igen_ers)
f8e779c983 antn*0366        if (gencost_barfile(igen_ers)(1:5).EQ.'m_eta')
025a9bb173 antn*0367      &    igen_etagcm=igen_ers
f09238ab8f Gael*0368       endif
                0369 c gfofi, gfosta, gfoper
                0370       if (igen_gfo.GT.0) then
                0371        ilps=ilnblnk( gencost_datafile(igen_gfo) )
                0372        gfofi=gencost_datafile(igen_gfo)(1:ilps)
f8e779c983 antn*0373        call cal_CopyDate(gencost_startdate(1,igen_gfo),gfosta,myThid)
f09238ab8f Gael*0374        gfoper=gencost_period(igen_gfo)
f8e779c983 antn*0375        if (gencost_barfile(igen_gfo)(1:5).EQ.'m_eta')
025a9bb173 antn*0376      &    igen_etagcm=igen_gfo
f09238ab8f Gael*0377       endif
                0378 
6b881d9117 Gael*0379 c mdt_y0, mdt_y1
025a9bb173 antn*0380 C in the case the model run is within the valid mdt years,
                0381 C we use the max/min check below to get the bound instead of
                0382 C basing on mdt years alone.
                0383 c       mdt_y0=mdtsta(1)/10000
                0384 c       mdt_y1=mdtend(1)/10000
                0385        mdt_y0=max(mdtsta(1)/10000,modelstartdate(1)/10000)
                0386        mdt_y1=min(mdtend(1)/10000,modelenddate(1)/10000)
6b881d9117 Gael*0387 
                0388        write(msgbuf,'(a,i8,a,i8)')
                0389      &  'mdt[start,end]date(1): ', mdtsta(1), ',', mdtend(1)
                0390         call print_message( msgbuf, standardmessageunit,
f8e779c983 antn*0391      &                      SQUEEZE_RIGHT , myThid)
6b881d9117 Gael*0392         write(msgbuf,'(a,i4,a,i4)')
                0393      &  '  TP MDT yrrange:          ', mdt_y0, ',', mdt_y1
                0394         call print_message( msgbuf, standardmessageunit,
f8e779c983 antn*0395      &                      SQUEEZE_RIGHT , myThid)
6b881d9117 Gael*0396 
f1e8c1d7ee Gael*0397 c--   First, read tiled data.
                0398       doglobalread = .false.
                0399       ladinit      = .false.
                0400 
025a9bb173 antn*0401       ilps=ilnblnk( gencost_barfile(igen_etagcm) )
de57a2ec4b Mart*0402       write(fname,'(2a,i10.10)')
025a9bb173 antn*0403      &     gencost_barfile(igen_etagcm)(1:ilps),'.',eccoiter
d699b51f70 Gael*0404 
f1e8c1d7ee Gael*0405 cgf =======================================================
                0406 cgf PART 1:
d39bcf1629 An T*0407 cgf        x Get the MDT (mdt) ... compute the sample mean
52e9fa3441 Gael*0408 cgf        (mean_slaobs_mdt) of the SLA data (i.e. RADS for tp, ers, and gfo
f1e8c1d7ee Gael*0409 cgf        together) over the time interval of the MDT ... subtract
52e9fa3441 Gael*0410 cgf        mean_slaobs_mdt from mdt.
d39bcf1629 An T*0411 cgf        x At this point, mdt is the inferred SLA reference field.
                0412 cgf        x From there on, mdt+sla will be directly comparable to
025a9bb173 antn*0413 cgf        the model SSH (etagcm).
f1e8c1d7ee Gael*0414 cgf =======================================================
                0415 
                0416 c--   Read mean field and mask
f09238ab8f Gael*0417 
025a9bb173 antn*0418       write(msgbuf,'(a,l)') ' sshv4:using_mdt=',using_mdt
                0419       call print_message( msgbuf, standardmessageunit,
                0420      &                    SQUEEZE_RIGHT , myThid)
f09238ab8f Gael*0421       if (using_mdt)
f8e779c983 antn*0422      &  CALL READ_REC_3D_RL( mdtfi, cost_iprec, 1,
                0423      &                       mdtob, 1, 1, myThid )
f09238ab8f Gael*0424 
                0425       factor =  0.01 _d 0
                0426       spval  = -9990. _d 0
                0427 
                0428       do bj = jtlo,jthi
                0429         do bi = itlo,ithi
                0430           do j = jmin,jmax
                0431             do i = imin,imax
0a1233e81a Ou W*0432               if ( (maskInC(i,j,bi,bj) .EQ. 0.).OR.
f09238ab8f Gael*0433 #ifndef ALLOW_SHALLOW_ALTIMETRY
                0434      &             (R_low(i,j,bi,bj).GT.-200.).OR.
                0435 #endif
0a1233e81a Ou W*0436      &             (mdtob(i,j,bi,bj) .LT. spval ).OR.
                0437      &             (mdtob(i,j,bi,bj) .EQ. 0. _d 0) ) then
f09238ab8f Gael*0438                 mdtma(i,j,bi,bj) = 0. _d 0
                0439                 mdtob(i,j,bi,bj) = 0. _d 0
                0440               else
                0441                 mdtma(i,j,bi,bj) = 1. _d 0
                0442                 mdtob(i,j,bi,bj) = mdtob(i,j,bi,bj)*factor
                0443               endif
                0444             enddo
                0445           enddo
                0446         enddo
                0447       enddo
f1e8c1d7ee Gael*0448 
52e9fa3441 Gael*0449 c--   Compute mean_slaobs_mdt: sample mean SLA over the time period of mdt.
f1e8c1d7ee Gael*0450 
                0451        do bj = jtlo,jthi
                0452         do bi = itlo,ithi
                0453          do j = jmin,jmax
                0454           do i = imin,imax
52e9fa3441 Gael*0455               mean_slaobs_mdt(i,j,bi,bj)  = 0. _d 0
f1e8c1d7ee Gael*0456               mean_slaobs_NUM(i,j,bi,bj)  = 0. _d 0
                0457           enddo
                0458          enddo
                0459         enddo
                0460        enddo
                0461 
d39bcf1629 An T*0462       do year=mdt_y0,mdt_y1
025a9bb173 antn*0463        do irec=1,nrec_year
f09238ab8f Gael*0464 
14be08a18d Gael*0465        do bj = jtlo,jthi
                0466         do bi = itlo,ithi
                0467          do j = jmin,jmax
                0468           do i = imin,imax
                0469               tpma(i,j,bi,bj) = 0. _d 0
                0470               tpob(i,j,bi,bj) = 0. _d 0
                0471               ersma(i,j,bi,bj) = 0. _d 0
                0472               ersob(i,j,bi,bj) = 0. _d 0
                0473               gfoma(i,j,bi,bj) = 0. _d 0
                0474               gfoob(i,j,bi,bj) = 0. _d 0
                0475           enddo
                0476          enddo
                0477         enddo
                0478        enddo
                0479 
025a9bb173 antn*0480 #ifdef ECCO_VERBOSE
                0481        write(msgbuf,'(a,3(1x,l))') ' sshv4:using_[tpj,ers,gfo]: ',
                0482      &   using_tpj, using_ers, using_gfo
                0483        call print_message( msgbuf, standardmessageunit,
                0484      &                     SQUEEZE_RIGHT , myThid)
                0485 #endif
f09238ab8f Gael*0486        if (using_tpj)
025a9bb173 antn*0487      & call cost_sla_read_yd( tpfi, tpsta, use_mon,
f09238ab8f Gael*0488      &                tpob, tpma,
025a9bb173 antn*0489      &                year, irec, myThid )
f09238ab8f Gael*0490        if (using_ers)
025a9bb173 antn*0491      & call cost_sla_read_yd( ersfi, erssta, use_mon,
f09238ab8f Gael*0492      &                ersob, ersma,
025a9bb173 antn*0493      &                year, irec, myThid )
f09238ab8f Gael*0494        if (using_gfo)
025a9bb173 antn*0495      & call cost_sla_read_yd( gfofi, gfosta, use_mon,
f09238ab8f Gael*0496      &                gfoob, gfoma,
025a9bb173 antn*0497      &                year, irec, myThid )
f09238ab8f Gael*0498 
0a1233e81a Ou W*0499        if (igen_tp .GT. 0) then
11c3150c71 Mart*0500         do bj = jtlo,jthi
                0501          do bi = itlo,ithi
                0502           do j = jmin,jmax
                0503            do i = imin,imax
f09238ab8f Gael*0504       if ( tpma(i,j,bi,bj)*mdtma(i,j,bi,bj)*
f1e8c1d7ee Gael*0505      &    gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then
52e9fa3441 Gael*0506           mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+
f09238ab8f Gael*0507      &  tpob(i,j,bi,bj)
f1e8c1d7ee Gael*0508           mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
                0509       endif
11c3150c71 Mart*0510            enddo
                0511           enddo
                0512          enddo
                0513         enddo
                0514        endif
0a1233e81a Ou W*0515        if (igen_ers .GT. 0) then
11c3150c71 Mart*0516         do bj = jtlo,jthi
                0517          do bi = itlo,ithi
                0518           do j = jmin,jmax
                0519            do i = imin,imax
f09238ab8f Gael*0520       if ( ersma(i,j,bi,bj)*mdtma(i,j,bi,bj)*
f1e8c1d7ee Gael*0521      &    gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then
52e9fa3441 Gael*0522           mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+
f09238ab8f Gael*0523      &  ersob(i,j,bi,bj)
f1e8c1d7ee Gael*0524           mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
                0525       endif
11c3150c71 Mart*0526            enddo
                0527           enddo
                0528          enddo
                0529         enddo
                0530        endif
0a1233e81a Ou W*0531        if (igen_gfo .GT. 0) then
11c3150c71 Mart*0532         do bj = jtlo,jthi
                0533          do bi = itlo,ithi
                0534           do j = jmin,jmax
                0535            do i = imin,imax
f09238ab8f Gael*0536       if ( gfoma(i,j,bi,bj)*mdtma(i,j,bi,bj)*
f1e8c1d7ee Gael*0537      &    gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then
52e9fa3441 Gael*0538           mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+
f09238ab8f Gael*0539      &  gfoob(i,j,bi,bj)
f1e8c1d7ee Gael*0540           mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
                0541       endif
11c3150c71 Mart*0542            enddo
f1e8c1d7ee Gael*0543           enddo
                0544          enddo
                0545         enddo
11c3150c71 Mart*0546        endif
f1e8c1d7ee Gael*0547 
025a9bb173 antn*0548        enddo !do irec=1,nrec_year
d39bcf1629 An T*0549       enddo !do year=mdt_y0,mdt_y1
f1e8c1d7ee Gael*0550 
                0551       do bj = jtlo,jthi
                0552         do bi = itlo,ithi
                0553           do j = jmin,jmax
                0554             do i = imin,imax
                0555                if ( ( mean_slaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
                0556      &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
f09238ab8f Gael*0557      &              ( mdtma(i,j,bi,bj) .NE. 0. ) ) then
52e9fa3441 Gael*0558                   mean_slaobs_mdt(i,j,bi,bj) =
                0559      &                 mean_slaobs_mdt(i,j,bi,bj) /
f1e8c1d7ee Gael*0560      &                 mean_slaobs_NUM(i,j,bi,bj)
                0561                else
52e9fa3441 Gael*0562                   mean_slaobs_mdt(i,j,bi,bj) = 0. _d 0
f1e8c1d7ee Gael*0563                endif
                0564             enddo
                0565           enddo
                0566         enddo
                0567       enddo
                0568 
52e9fa3441 Gael*0569 c--   smooth mean_slaobs_mdt:
f1e8c1d7ee Gael*0570 
0a1233e81a Ou W*0571       if (igen_mdt.GT.0) then
11c3150c71 Mart*0572        if (gencost_outputlevel(igen_mdt).GT.0) then
de57a2ec4b Mart*0573         write(fname4test,'(1a)') 'sla2mdt_raw'
f8e779c983 antn*0574         CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                0575      &                        mean_slaobs_mdt, 1, 1, myThid )
11c3150c71 Mart*0576        endif
f1e8c1d7ee Gael*0577 
adf79a5d72 Gael*0578 #ifdef ALLOW_SMOOTH
0a1233e81a Ou W*0579        if ( useSMOOTH.AND.(k2_mdt.GT.0) ) then
                0580         call smooth_hetero2d(mean_slaobs_mdt,maskInC,
3248fb8b5a Gael*0581      &     gencost_posproc_c(k2_mdt,igen_mdt),
f8e779c983 antn*0582      &     gencost_posproc_i(k2_mdt,igen_mdt),myThid)
f1e8c1d7ee Gael*0583 
0a1233e81a Ou W*0584         if (gencost_outputlevel(igen_mdt).GT.0) then
de57a2ec4b Mart*0585          write(fname4test,'(1a)') 'sla2mdt_smooth'
0a1233e81a Ou W*0586          CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                0587      &                         mean_slaobs_mdt, 1, 1, myThid )
                0588         endif
11c3150c71 Mart*0589        endif
0a1233e81a Ou W*0590 #endif
                0591 c     endif igen_mdt .gt. 0
6b881d9117 Gael*0592       endif
f1e8c1d7ee Gael*0593 
d39bcf1629 An T*0594 c--   re-reference mdt:
f1e8c1d7ee Gael*0595       do bj = jtlo,jthi
                0596         do bi = itlo,ithi
                0597           do j = jmin,jmax
                0598             do i = imin,imax
f09238ab8f Gael*0599                if ( ( mdtma(i,j,bi,bj) .NE. 0. ).AND.
2a0641f1b7 Gael*0600      &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
                0601      &              ( doReference ) ) then
f09238ab8f Gael*0602                   mdtob(i,j,bi,bj) = mdtob(i,j,bi,bj)
52e9fa3441 Gael*0603      &                 -mean_slaobs_mdt(i,j,bi,bj)
f1e8c1d7ee Gael*0604                endif
                0605             enddo
                0606           enddo
                0607         enddo
                0608       enddo
                0609 
                0610 cgf =======================================================
025a9bb173 antn*0611 cgf PART 2: compute sample means of etagcm-slaobs over the
                0612 cgf         period that is covered by the model (i.e. etagcm).
11c3150c71 Mart*0613 cgf       x for all SLA data sets together: mean_psMssh_all, mean_psMssh_all_MSK,
                0614 cgf         and offset will be used in PART 3 (MDT cost term).
                0615 cgf       x for each SLA data individually. mean_psMtpobs, mean_psMtpobs_MS, etc.
                0616 cgf         will be used in PART 4&5 (SLA cost terms).
f1e8c1d7ee Gael*0617 cgf =======================================================
                0618 
                0619       do bj = jtlo,jthi
                0620         do bi = itlo,ithi
                0621           do j = jmin,jmax
                0622             do i = imin,imax
                0623 
                0624               psmean(i,j,bi,bj)    = 0. _d 0
52e9fa3441 Gael*0625               mean_psMslaobs(i,j,bi,bj)  = 0. _d 0
f1e8c1d7ee Gael*0626               mean_psMtpobs(i,j,bi,bj)  = 0. _d 0
                0627               mean_psMersobs(i,j,bi,bj) = 0. _d 0
                0628               mean_psMgfoobs(i,j,bi,bj) = 0. _d 0
                0629               mean_psMssh_all(i,j,bi,bj) = 0. _d 0
52e9fa3441 Gael*0630               mean_slaobs_model(i,j,bi,bj)  = 0. _d 0
f1e8c1d7ee Gael*0631 
                0632               mean_psMtpobs_NUM(i,j,bi,bj)  = 0. _d 0
                0633               mean_psMersobs_NUM(i,j,bi,bj) = 0. _d 0
                0634               mean_psMgfoobs_NUM(i,j,bi,bj) = 0. _d 0
                0635               mean_psMssh_all_NUM(i,j,bi,bj) = 0. _d 0
52e9fa3441 Gael*0636               mean_psMslaobs_MSK(i,j,bi,bj)  = 0. _d 0
f1e8c1d7ee Gael*0637 
                0638             enddo
                0639           enddo
                0640         enddo
                0641       enddo
                0642       offset     = 0. _d 0
                0643       offset_sum = 0. _d 0
                0644 
025a9bb173 antn*0645 #ifdef ECCO_VERBOSE
                0646       write(msgbuf,'(a,3(1x,i6))')
                0647      & 'sshv4: before entering irec loop, n[mon,day]srec,nrec_total: ',
                0648      & nmonsrec,ndaysrec,nrec_total
                0649       call print_message( msgbuf, standardmessageunit,
                0650      &                    SQUEEZE_RIGHT , myThid)
                0651 #endif
                0652 
                0653       do irec = 1, nrec_total
f1e8c1d7ee Gael*0654 
101f75e5cd Gael*0655 #ifdef ALLOW_AUTODIFF
025a9bb173 antn*0656         call active_read_xy( fname, etagcm, irec, doglobalread,
f8e779c983 antn*0657      &                       ladinit, eccoiter, myThid,
025a9bb173 antn*0658      &                       gencost_dummy(igen_etagcm) )
101f75e5cd Gael*0659 #else
025a9bb173 antn*0660         CALL READ_REC_XY_RL( fname, etagcm, iRec, 1, myThid )
101f75e5cd Gael*0661 #endif
                0662 
f8e779c983 antn*0663         if (.NOT.useEtaMean)
025a9bb173 antn*0664      &  CALL REMOVE_MEAN_RL( 1, etagcm, maskInC, maskInC, rA, drF,
                0665      &                       'etagcm', zeroRL, 1, myThid )
0699f1d2f2 Gael*0666 
14be08a18d Gael*0667        do bj = jtlo,jthi
                0668         do bi = itlo,ithi
                0669          do j = jmin,jmax
                0670           do i = imin,imax
                0671               tpma(i,j,bi,bj) = 0. _d 0
                0672               tpob(i,j,bi,bj) = 0. _d 0
                0673               ersma(i,j,bi,bj) = 0. _d 0
                0674               ersob(i,j,bi,bj) = 0. _d 0
                0675               gfoma(i,j,bi,bj) = 0. _d 0
                0676               gfoob(i,j,bi,bj) = 0. _d 0
                0677           enddo
                0678          enddo
                0679         enddo
                0680        enddo
                0681 
f09238ab8f Gael*0682         if (using_tpj)
025a9bb173 antn*0683      &  call cost_sla_read( tpfi, tpsta, tpper, use_mon,
f09238ab8f Gael*0684      &                zeroRL, zeroRL,
                0685      &                tpob, tpma,
f8e779c983 antn*0686      &                irec, myThid )
f09238ab8f Gael*0687         if (using_ers)
025a9bb173 antn*0688      &  call cost_sla_read( ersfi, erssta, ersper, use_mon,
f09238ab8f Gael*0689      &                zeroRL, zeroRL,
                0690      &                ersob, ersma,
f8e779c983 antn*0691      &                irec, myThid )
f09238ab8f Gael*0692         if (using_gfo)
025a9bb173 antn*0693      &  call cost_sla_read( gfofi, gfosta, gfoper, use_mon,
f09238ab8f Gael*0694      &                zeroRL, zeroRL,
                0695      &                gfoob, gfoma,
f8e779c983 antn*0696      &                irec, myThid )
f1e8c1d7ee Gael*0697 
                0698         do bj = jtlo,jthi
                0699           do bi = itlo,ithi
                0700             do j = jmin,jmax
                0701               do i = imin,imax
                0702                 psmean(i,j,bi,bj) = psmean(i,j,bi,bj) +
025a9bb173 antn*0703      &                etagcm(i,j,bi,bj) / float(nrec_total)
0a1233e81a Ou W*0704                 if ( igen_tp .GT. 0 ) then
11c3150c71 Mart*0705                   if ( tpma(i,j,bi,bj)*
f1e8c1d7ee Gael*0706      &             gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then
52e9fa3441 Gael*0707                    mean_slaobs_model(i,j,bi,bj)=
f09238ab8f Gael*0708      &              mean_slaobs_model(i,j,bi,bj)+tpob(i,j,bi,bj)
f1e8c1d7ee Gael*0709                    mean_psMtpobs(i,j,bi,bj) =
                0710      &                 mean_psMtpobs(i,j,bi,bj) +
025a9bb173 antn*0711      &                 etagcm(i,j,bi,bj)-tpob(i,j,bi,bj)
f1e8c1d7ee Gael*0712                    mean_psMtpobs_NUM(i,j,bi,bj) =
                0713      &                 mean_psMtpobs_NUM(i,j,bi,bj) + 1. _d 0
11c3150c71 Mart*0714                   endif
f1e8c1d7ee Gael*0715                 endif
0a1233e81a Ou W*0716                 if ( igen_ers .GT. 0 ) then
11c3150c71 Mart*0717                   if ( ersma(i,j,bi,bj)*
f1e8c1d7ee Gael*0718      &             gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then
52e9fa3441 Gael*0719                    mean_slaobs_model(i,j,bi,bj)=
f09238ab8f Gael*0720      &              mean_slaobs_model(i,j,bi,bj)+ersob(i,j,bi,bj)
f1e8c1d7ee Gael*0721                    mean_psMersobs(i,j,bi,bj) =
                0722      &                 mean_psMersobs(i,j,bi,bj) +
025a9bb173 antn*0723      &                 etagcm(i,j,bi,bj)-ersob(i,j,bi,bj)
f1e8c1d7ee Gael*0724                    mean_psMersobs_NUM(i,j,bi,bj) =
                0725      &                 mean_psMersobs_NUM(i,j,bi,bj) + 1. _d 0
11c3150c71 Mart*0726                   endif
f1e8c1d7ee Gael*0727                 endif
0a1233e81a Ou W*0728                 if ( igen_gfo .GT. 0 ) then
11c3150c71 Mart*0729                   if ( gfoma(i,j,bi,bj)*
f1e8c1d7ee Gael*0730      &             gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then
52e9fa3441 Gael*0731                    mean_slaobs_model(i,j,bi,bj)=
f09238ab8f Gael*0732      &              mean_slaobs_model(i,j,bi,bj)+gfoob(i,j,bi,bj)
f1e8c1d7ee Gael*0733                    mean_psMgfoobs(i,j,bi,bj) =
                0734      &                 mean_psMgfoobs(i,j,bi,bj) +
025a9bb173 antn*0735      &                 etagcm(i,j,bi,bj)-gfoob(i,j,bi,bj)
f1e8c1d7ee Gael*0736                    mean_psMgfoobs_NUM(i,j,bi,bj) =
                0737      &                 mean_psMgfoobs_NUM(i,j,bi,bj) + 1. _d 0
11c3150c71 Mart*0738                   endif
f1e8c1d7ee Gael*0739                 endif
                0740               enddo
                0741             enddo
                0742           enddo
                0743         enddo
                0744 
025a9bb173 antn*0745 c--   END loop over irec records for the first time.
f1e8c1d7ee Gael*0746       enddo
                0747 
f8e779c983 antn*0748       call ecco_zero(num0array,1,oneRL,myThid)
6c0b33a90f An T*0749 
f1e8c1d7ee Gael*0750         do bj = jtlo,jthi
                0751           do bi = itlo,ithi
                0752             do j = jmin,jmax
                0753               do i = imin,imax
70fde697ea Gael*0754                if ( ( mean_psMtpobs_NUM(i,j,bi,bj) .NE. 0. )
                0755      &              ) then
f1e8c1d7ee Gael*0756                   mean_psMssh_all(i,j,bi,bj) =
                0757      &                 mean_psMssh_all(i,j,bi,bj) +
                0758      &                 mean_psMtpobs(i,j,bi,bj)
                0759                   mean_psMssh_all_NUM(i,j,bi,bj) =
                0760      &                 mean_psMssh_all_NUM(i,j,bi,bj) +
                0761      &                 mean_psMtpobs_NUM(i,j,bi,bj)
                0762                   mean_psMtpobs(i,j,bi,bj) =
                0763      &                 mean_psMtpobs(i,j,bi,bj) /
                0764      &                 mean_psMtpobs_NUM(i,j,bi,bj)
                0765                endif
70fde697ea Gael*0766                if ( ( mean_psMersobs_NUM(i,j,bi,bj) .NE. 0. )
                0767      &              ) then
f1e8c1d7ee Gael*0768                   mean_psMssh_all(i,j,bi,bj) =
                0769      &                 mean_psMssh_all(i,j,bi,bj) +
                0770      &                 mean_psMersobs(i,j,bi,bj)
                0771                   mean_psMssh_all_NUM(i,j,bi,bj) =
                0772      &                 mean_psMssh_all_NUM(i,j,bi,bj) +
                0773      &                 mean_psMersobs_NUM(i,j,bi,bj)
                0774                   mean_psMersobs(i,j,bi,bj) =
                0775      &                 mean_psMersobs(i,j,bi,bj) /
                0776      &                 mean_psMersobs_NUM(i,j,bi,bj)
                0777                endif
70fde697ea Gael*0778                if ( ( mean_psMgfoobs_NUM(i,j,bi,bj) .NE. 0. )
                0779      &              ) then
f1e8c1d7ee Gael*0780                   mean_psMssh_all(i,j,bi,bj) =
                0781      &                 mean_psMssh_all(i,j,bi,bj) +
                0782      &                 mean_psMgfoobs(i,j,bi,bj)
                0783                   mean_psMssh_all_NUM(i,j,bi,bj) =
                0784      &                 mean_psMssh_all_NUM(i,j,bi,bj) +
                0785      &                 mean_psMgfoobs_NUM(i,j,bi,bj)
                0786                   mean_psMgfoobs(i,j,bi,bj) =
                0787      &                 mean_psMgfoobs(i,j,bi,bj) /
                0788      &                 mean_psMgfoobs_NUM(i,j,bi,bj)
                0789                endif
52e9fa3441 Gael*0790 
                0791                if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .GT. 0 ).AND.
76d5f72256 Gael*0792      &              ( maskc(i,j,1,bi,bj) .NE. 0. )
52e9fa3441 Gael*0793      &              ) then
                0794                   mean_psMslaobs(i,j,bi,bj) =
                0795      &                 mean_psMssh_all(i,j,bi,bj) /
                0796      &                 mean_psMssh_all_NUM(i,j,bi,bj)
                0797                   mean_psMslaobs_MSK(i,j,bi,bj) = 1. _d 0
                0798                else
                0799                   mean_psMslaobs(i,j,bi,bj) = 0. _d 0
                0800                   mean_psMslaobs_MSK(i,j,bi,bj) = 0. _d 0
                0801                endif
                0802 
204eb34df4 Gael*0803                if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .GT. num0 ).AND.
f1e8c1d7ee Gael*0804      &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
f09238ab8f Gael*0805      &              ( mdtma(i,j,bi,bj) .NE. 0. ).AND.
2a0641f1b7 Gael*0806      &              ( doReference ) ) then
52e9fa3441 Gael*0807                   mean_slaobs_model(i,j,bi,bj) =
                0808      &                 mean_slaobs_model(i,j,bi,bj) /
f1e8c1d7ee Gael*0809      &                 mean_psMssh_all_NUM(i,j,bi,bj)
                0810                   mean_psMssh_all(i,j,bi,bj) =
                0811      &                 mean_psMssh_all(i,j,bi,bj) /
f09238ab8f Gael*0812      &                 mean_psMssh_all_NUM(i,j,bi,bj)-mdtob(i,j,bi,bj)
f1e8c1d7ee Gael*0813                   mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
204eb34df4 Gael*0814                   offset=offset+RA(i,j,bi,bj)*mean_psMssh_all(i,j,bi,bj)
                0815                   offset_sum=offset_sum+RA(i,j,bi,bj)
f09238ab8f Gael*0816                elseif ( ( mdtma(i,j,bi,bj) .NE. 0. ) .AND.
2a0641f1b7 Gael*0817      &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
                0818      &              ( .NOT.doReference ) ) then
52e9fa3441 Gael*0819                   mean_slaobs_model(i,j,bi,bj) = 0.d0
2a0641f1b7 Gael*0820                   mean_psMssh_all(i,j,bi,bj) = 0. _d 0
                0821                   mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
                0822                   offset=offset+RA(i,j,bi,bj)
f09238ab8f Gael*0823      &                  *( psmean(i,j,bi,bj) -mdtob(i,j,bi,bj) )
2a0641f1b7 Gael*0824                   offset_sum=offset_sum+RA(i,j,bi,bj)
6c0b33a90f An T*0825                   num0array(i,j,bi,bj)=0. _d 0
f1e8c1d7ee Gael*0826                else
52e9fa3441 Gael*0827                   mean_slaobs_model(i,j,bi,bj) = 0.d0
f1e8c1d7ee Gael*0828                   mean_psMssh_all(i,j,bi,bj) = 0. _d 0
                0829                   mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0
6c0b33a90f An T*0830                   num0array(i,j,bi,bj)=0. _d 0
f1e8c1d7ee Gael*0831                endif
                0832               enddo
                0833             enddo
                0834           enddo
                0835         enddo
                0836 
                0837 c--   Do a global summation.
f8e779c983 antn*0838       _GLOBAL_SUM_RL( offset     , myThid )
                0839       _GLOBAL_SUM_RL( offset_sum , myThid )
f1e8c1d7ee Gael*0840 
025a9bb173 antn*0841 C --  add warning that written mean_slaobs_model is all zero
                0842 c     due to having less data points that hardcoded num0
0a1233e81a Ou W*0843         num0total=0. _d 0
6c0b33a90f An T*0844         do bj = jtlo,jthi
                0845           do bi = itlo,ithi
                0846             do j = jmin,jmax
                0847               do i = imin,imax
                0848                 num0total=num0total+num0array(i,j,bi,bj)
                0849               enddo
                0850             enddo
                0851           enddo
                0852         enddo
0a1233e81a Ou W*0853       if(num0total.LT.1. _d 0) then
f8e779c983 antn*0854         write(msgbuf,'(A,i5,A)')
6c0b33a90f An T*0855      &  '** WARNING ** S/R COST_GENCOST_SSHV4: There are <',num0,' pts'
                0856       call print_message( msgbuf, errormessageunit,
f8e779c983 antn*0857      &                    SQUEEZE_RIGHT , myThid)
                0858         write(msgbuf,'(A)')
6c0b33a90f An T*0859      &  '                  at all grid points for combined tp+ers+gfo.'
                0860       call print_message( msgbuf, errormessageunit,
f8e779c983 antn*0861      &                    SQUEEZE_RIGHT , myThid)
                0862         write(msgbuf,'(A)')
6c0b33a90f An T*0863      &  'So, all model slaob minus model sla2model_raw are set to 0.'
                0864       call print_message( msgbuf, errormessageunit,
f8e779c983 antn*0865      &                    SQUEEZE_RIGHT , myThid)
6c0b33a90f An T*0866       endif
                0867 
0a1233e81a Ou W*0868       if (igen_gmsl.GT.0) then
                0869        if (offset_sum .EQ. 0.0) then
cf137f827f Gael*0870         if (gencost_outputlevel(igen_gmsl).GT.0) then
f8e779c983 antn*0871         _BEGIN_MASTER( myThid )
e5310f7c13 Gael*0872         write(msgbuf,'(a)') ' sshv4: offset_sum = zero!'
f1e8c1d7ee Gael*0873         call print_message( msgbuf, standardmessageunit,
f8e779c983 antn*0874      &                          SQUEEZE_RIGHT , myThid)
                0875         _END_MASTER( myThid )
cf137f827f Gael*0876         endif
0a1233e81a Ou W*0877 c       stop   '  ... stopped in cost_ssh.'
11c3150c71 Mart*0878        else
cf137f827f Gael*0879         if (gencost_outputlevel(igen_gmsl).GT.0) then
f8e779c983 antn*0880         _BEGIN_MASTER( myThid )
bc9bce397d Gael*0881         write(msgbuf,'(a,2d22.15)') ' sshv4:offset,sum=',
                0882      &      offset,offset_sum
f1e8c1d7ee Gael*0883         call print_message( msgbuf, standardmessageunit,
f8e779c983 antn*0884      &                          SQUEEZE_RIGHT , myThid)
                0885         _END_MASTER( myThid )
cf137f827f Gael*0886         endif
11c3150c71 Mart*0887        endif
0a1233e81a Ou W*0888 c     endif igen_gmsl .gt. 0
f1e8c1d7ee Gael*0889       endif
                0890 
                0891 c--   Compute (average) offset
0a1233e81a Ou W*0892       if ( offset_sum.NE.0. _d 0 ) then
                0893         offset = offset / offset_sum
                0894       else
                0895         offset = 0. _d 0
                0896       endif
f1e8c1d7ee Gael*0897 
                0898 c--   subtract offset from mean_psMssh_all and psmean
                0899       do bj = jtlo,jthi
                0900         do bi = itlo,ithi
                0901           do j = jmin,jmax
                0902             do i = imin,imax
                0903 
204eb34df4 Gael*0904                if ( (mean_psMssh_all_NUM(i,j,bi,bj) .GT. num0) .AND.
f09238ab8f Gael*0905      &              ( mdtma(i,j,bi,bj) .NE. 0. ) .AND.
f8e779c983 antn*0906      &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
204eb34df4 Gael*0907      &              ( doReference ) ) then
f1e8c1d7ee Gael*0908 c use the re-referencing approach
204eb34df4 Gael*0909                    mean_psMssh_all(i,j,bi,bj) =
f1e8c1d7ee Gael*0910      &                 mean_psMssh_all(i,j,bi,bj) - offset
204eb34df4 Gael*0911                    mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
                0912 
f09238ab8f Gael*0913                elseif ( ( mdtma(i,j,bi,bj) .NE. 0. ) .AND.
204eb34df4 Gael*0914      &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
76d5f72256 Gael*0915      &              ( .NOT.doReference ) ) then
2a0641f1b7 Gael*0916 c use the simpler approach
                0917                    mean_psMssh_all(i,j,bi,bj) =
f09238ab8f Gael*0918      &             psmean(i,j,bi,bj) -mdtob(i,j,bi,bj) - offset
204eb34df4 Gael*0919                    mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
                0920 
2a0641f1b7 Gael*0921                else
                0922                    mean_psMssh_all(i,j,bi,bj) = 0. _d 0
204eb34df4 Gael*0923                    mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0
f883c17701 Gael*0924                endif
2a0641f1b7 Gael*0925 
f1e8c1d7ee Gael*0926                if ( maskc(i,j,1,bi,bj) .NE. 0. )
                0927      &            psmean(i,j,bi,bj)=psmean(i,j,bi,bj)-offset
f883c17701 Gael*0928 
f1e8c1d7ee Gael*0929             enddo
                0930           enddo
                0931         enddo
11c3150c71 Mart*0932       enddo
f1e8c1d7ee Gael*0933 
                0934 c--    smooth mean_psMssh_all
0a1233e81a Ou W*0935       if (igen_mdt.GT.0) then
11c3150c71 Mart*0936        if (gencost_outputlevel(igen_mdt).GT.0) then
de57a2ec4b Mart*0937         write(fname4test,'(1a)') 'mdtdiff_raw'
f8e779c983 antn*0938         CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                0939      &                        mean_psMssh_all, 1, 1, myThid )
de57a2ec4b Mart*0940         write(fname4test,'(1a)') 'sla2model_raw'
f8e779c983 antn*0941         CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                0942      &                        mean_slaobs_model, 1, 1, myThid )
11c3150c71 Mart*0943        endif
f1e8c1d7ee Gael*0944 
adf79a5d72 Gael*0945 #ifdef ALLOW_SMOOTH
0a1233e81a Ou W*0946        if ( useSMOOTH.AND.(k2_mdt.GT.0) ) then
                0947         call smooth_hetero2d(mean_psMssh_all,maskInC,
                0948      &     gencost_posproc_c(k2_mdt,igen_mdt),
                0949      &     gencost_posproc_i(k2_mdt,igen_mdt),myThid)
                0950         call smooth_hetero2d(mean_slaobs_model,maskInC,
3248fb8b5a Gael*0951      &     gencost_posproc_c(k2_mdt,igen_mdt),
f8e779c983 antn*0952      &     gencost_posproc_i(k2_mdt,igen_mdt),myThid)
f1e8c1d7ee Gael*0953 
0a1233e81a Ou W*0954         if (gencost_outputlevel(igen_mdt).GT.0) then
de57a2ec4b Mart*0955          write(fname4test,'(1a)') 'mdtdiff_smooth'
0a1233e81a Ou W*0956          CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                0957      &                         mean_psMssh_all, 1, 1, myThid )
de57a2ec4b Mart*0958          write(fname4test,'(1a)') 'sla2model_smooth'
0a1233e81a Ou W*0959          CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                0960      &                         mean_slaobs_model, 1, 1, myThid )
                0961         endif
11c3150c71 Mart*0962        endif
0a1233e81a Ou W*0963 #endif
f1e8c1d7ee Gael*0964 
                0965 cgf at this point:
025a9bb173 antn*0966 cgf     1) mean_psMssh_all is the sample mean <etagcm-slaobs-mdt-offset>,
f1e8c1d7ee Gael*0967 cgf             to which a smoothing filter has been applied.
025a9bb173 antn*0968 cgf     2) mean_psMtpobs is the (unsmoothed) sample mean <etagcm-tpobs>.
f1e8c1d7ee Gael*0969 cgf             And similarly for ers and gfo, each treated separately.
                0970 
                0971 cgf =======================================================
                0972 cgf PART 3: compute MDT cost term
                0973 cgf =======================================================
                0974 
11c3150c71 Mart*0975        do bj = jtlo,jthi
f1e8c1d7ee Gael*0976         do bi = itlo,ithi
11c3150c71 Mart*0977          do j = jmin,jmax
                0978           do i = imin,imax
204eb34df4 Gael*0979        if (mean_psMssh_all_MSK(i,j,bi,bj).NE.0. _d 0) then
                0980          junk = mean_psMssh_all(i,j,bi,bj)
                0981          junkweight = gencost_weight(i,j,bi,bj,igen_mdt)
f09238ab8f Gael*0982      &              * mdtma(i,j,bi,bj)
204eb34df4 Gael*0983        else
                0984          junk = 0. _d 0
                0985          junkweight = 0. _d 0
                0986        endif
61bb52ec0d Gael*0987        objf_gencost(bi,bj,igen_mdt) = objf_gencost(bi,bj,igen_mdt)
f1e8c1d7ee Gael*0988      &      + junk*junk*junkweight
0a1233e81a Ou W*0989        if ( junkweight .NE. 0. ) num_gencost(bi,bj,igen_mdt) =
61bb52ec0d Gael*0990      &      num_gencost(bi,bj,igen_mdt) + 1. _d 0
f1e8c1d7ee Gael*0991        diagnosfld(i,j,bi,bj) = junk*junk*junkweight
                0992           enddo
11c3150c71 Mart*0993          enddo
f1e8c1d7ee Gael*0994         enddo
11c3150c71 Mart*0995        enddo
f1e8c1d7ee Gael*0996 
11c3150c71 Mart*0997        if (gencost_outputlevel(igen_mdt).GT.0) then
de57a2ec4b Mart*0998         write(fname4test,'(1a)') 'misfits_mdt'
f8e779c983 antn*0999         CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                1000      &                        diagnosfld, 1, 1, myThid )
11c3150c71 Mart*1001        endif
0a1233e81a Ou W*1002 c     endif igen_mdt .gt. 0
6b881d9117 Gael*1003       endif
f1e8c1d7ee Gael*1004 
                1005 cgf =======================================================
                1006 cgf PART 4: compute smooth SLA cost term
                1007 cgf =======================================================
                1008 
025a9bb173 antn*1009       if (use_mon) then
                1010         nrecsave=1
                1011       else
                1012         nrecsave=35
                1013       endif
f1e8c1d7ee Gael*1014 
025a9bb173 antn*1015 C add a warning if nrec_total is less than the hardcoded nrecsave
                1016       tempinteger=nrec_total-nrecsave+1
                1017 #ifdef ECCO_VERBOSE
                1018       write(msgbuf,'(A,i5,1x,i5,1x,i5)')
                1019      & 'sshv4: nrec[save,_total],tempinteger',
                1020      &     nrecsave,nrec_total,tempinteger
                1021 #endif
0a1233e81a Ou W*1022       if(tempinteger.LT.1) then
f8e779c983 antn*1023         write(msgbuf,'(A,i5)')
025a9bb173 antn*1024      &  '** WARNING ** S/R COST_GENCOST_SSHV4: There are < ',nrecsave
                1025         call print_message( msgbuf, errormessageunit,
                1026      &                      SQUEEZE_RIGHT , myThid)
                1027         if (use_mon) then
                1028          write(msgbuf,'(A)')
                1029      &   '        months required to calculate running mean tp+ers+gfo.'
                1030          call print_message( msgbuf, errormessageunit,
                1031      &                       SQUEEZE_RIGHT , myThid)
                1032         else
                1033          write(msgbuf,'(A)')
6c0b33a90f An T*1034      &   '        days required to calculate running mean tp+ers+gfo.'
025a9bb173 antn*1035          call print_message( msgbuf, errormessageunit,
                1036      &                       SQUEEZE_RIGHT , myThid)
                1037         endif
f8e779c983 antn*1038         write(msgbuf,'(A)')
6c0b33a90f An T*1039      &  'PART 4 in cost_gencost_sshv4 is skipped entirely, and there'
                1040       call print_message( msgbuf, errormessageunit,
f8e779c983 antn*1041      &                    SQUEEZE_RIGHT , myThid)
                1042         write(msgbuf,'(A)')
6c0b33a90f An T*1043      &  'will be NO report of [tp,ers,gfo,lsc,gmsl] in costfunction.'
                1044       call print_message( msgbuf, errormessageunit,
f8e779c983 antn*1045      &                    SQUEEZE_RIGHT , myThid)
6c0b33a90f An T*1046       endif
                1047 
025a9bb173 antn*1048       do irec = 1, nrec_total-nrecsave+1
f1e8c1d7ee Gael*1049 
                1050        do bj = jtlo,jthi
                1051         do bi = itlo,ithi
                1052          do j = jmin,jmax
                1053           do i = imin,imax
                1054               anom_psMslaobs(i,j,bi,bj)  = 0. _d 0
d699b51f70 Gael*1055               anom_psMslaobs_MSK(i,j,bi,bj)  = 0. _d 0
                1056               anom_psMslaobs_SUB(i,j,bi,bj)  = 0. _d 0
f1e8c1d7ee Gael*1057               anom_slaobs(i,j,bi,bj)  = 0. _d 0
d699b51f70 Gael*1058               anom_slaobs_SUB(i,j,bi,bj)  = 0. _d 0
f1e8c1d7ee Gael*1059               anom_psMslaobs_NUM(i,j,bi,bj)  = 0. _d 0
                1060           enddo
                1061          enddo
                1062         enddo
                1063        enddo
                1064 
025a9bb173 antn*1065 c PART 4.1: compute running sample average over nrecsave
f1e8c1d7ee Gael*1066 c ------------------------------------------------------
                1067 
025a9bb173 antn*1068       do jrec=1,nrecsave
f1e8c1d7ee Gael*1069 
                1070         krec=irec+jrec-1
                1071 
101f75e5cd Gael*1072 #ifdef ALLOW_AUTODIFF
025a9bb173 antn*1073         call active_read_xy( fname, etagcm, krec, doglobalread,
f8e779c983 antn*1074      &                       ladinit, eccoiter, myThid,
025a9bb173 antn*1075      &                       gencost_dummy(igen_etagcm) )
101f75e5cd Gael*1076 #else
025a9bb173 antn*1077         CALL READ_REC_XY_RL( fname, etagcm, kRec, 1, myThid )
101f75e5cd Gael*1078 #endif
f1e8c1d7ee Gael*1079 
6b881d9117 Gael*1080         if (.NOT.useEtaMean)
025a9bb173 antn*1081      &  CALL REMOVE_MEAN_RL( 1, etagcm, maskInC, maskInC, rA, drF,
                1082      &                       'etagcm', zeroRL, 1, myThid )
0699f1d2f2 Gael*1083 
14be08a18d Gael*1084        do bj = jtlo,jthi
                1085         do bi = itlo,ithi
                1086          do j = jmin,jmax
                1087           do i = imin,imax
                1088               tpma(i,j,bi,bj) = 0. _d 0
                1089               tpob(i,j,bi,bj) = 0. _d 0
                1090               ersma(i,j,bi,bj) = 0. _d 0
                1091               ersob(i,j,bi,bj) = 0. _d 0
                1092               gfoma(i,j,bi,bj) = 0. _d 0
                1093               gfoob(i,j,bi,bj) = 0. _d 0
                1094           enddo
                1095          enddo
                1096         enddo
                1097        enddo
                1098 
f09238ab8f Gael*1099         if (using_tpj)
025a9bb173 antn*1100      &  call cost_sla_read( tpfi, tpsta, tpper, use_mon,
f09238ab8f Gael*1101      &                zeroRL, zeroRL,
                1102      &                tpob, tpma,
f8e779c983 antn*1103      &                krec, myThid )
f09238ab8f Gael*1104         if (using_ers)
025a9bb173 antn*1105      &  call cost_sla_read( ersfi, erssta, ersper, use_mon,
f09238ab8f Gael*1106      &                zeroRL, zeroRL,
                1107      &                ersob, ersma,
f8e779c983 antn*1108      &                krec, myThid )
f09238ab8f Gael*1109         if (using_gfo)
025a9bb173 antn*1110      &  call cost_sla_read( gfofi, gfosta, gfoper, use_mon,
f09238ab8f Gael*1111      &                zeroRL, zeroRL,
                1112      &                gfoob, gfoma,
f8e779c983 antn*1113      &                krec, myThid )
f1e8c1d7ee Gael*1114 
                1115        do bj = jtlo,jthi
                1116         do bi = itlo,ithi
                1117          do j = jmin,jmax
                1118           do i = imin,imax
f09238ab8f Gael*1119       if ( tpma(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj)
f1e8c1d7ee Gael*1120      &  .NE.0. ) then
                1121            anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
025a9bb173 antn*1122      &        etagcm(i,j,bi,bj)-tpob(i,j,bi,bj)
52e9fa3441 Gael*1123      &        -mean_psMslaobs(i,j,bi,bj)
f1e8c1d7ee Gael*1124            anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
f09238ab8f Gael*1125      &        tpob(i,j,bi,bj)
f1e8c1d7ee Gael*1126            anom_psMslaobs_NUM(i,j,bi,bj)=
                1127      &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
025a9bb173 antn*1128            if ( (2*jrec.EQ.nrecsave).OR.(2*jrec-1.EQ.nrecsave) )
d699b51f70 Gael*1129      &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0
f1e8c1d7ee Gael*1130       endif
f09238ab8f Gael*1131       if ( ersma(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj)
f1e8c1d7ee Gael*1132      &  .NE.0. ) then
                1133            anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
025a9bb173 antn*1134      &        etagcm(i,j,bi,bj)-ersob(i,j,bi,bj)
52e9fa3441 Gael*1135      &        -mean_psMslaobs(i,j,bi,bj)
f1e8c1d7ee Gael*1136            anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
f09238ab8f Gael*1137      &        ersob(i,j,bi,bj)
f1e8c1d7ee Gael*1138            anom_psMslaobs_NUM(i,j,bi,bj)=
                1139      &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
025a9bb173 antn*1140            if ( (2*jrec.EQ.nrecsave).OR.(2*jrec-1.EQ.nrecsave) )
d699b51f70 Gael*1141      &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0
f1e8c1d7ee Gael*1142       endif
f09238ab8f Gael*1143       if ( gfoma(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj)
f1e8c1d7ee Gael*1144      &  .NE.0. ) then
                1145            anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
025a9bb173 antn*1146      &        etagcm(i,j,bi,bj)-gfoob(i,j,bi,bj)
52e9fa3441 Gael*1147      &        -mean_psMslaobs(i,j,bi,bj)
f1e8c1d7ee Gael*1148            anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
f09238ab8f Gael*1149      &        gfoob(i,j,bi,bj)
f1e8c1d7ee Gael*1150            anom_psMslaobs_NUM(i,j,bi,bj)=
                1151      &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
025a9bb173 antn*1152            if ( (2*jrec.EQ.nrecsave).OR.(2*jrec-1.EQ.nrecsave) )
d699b51f70 Gael*1153      &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0
f1e8c1d7ee Gael*1154       endif
                1155           enddo
                1156          enddo
                1157         enddo
                1158        enddo
                1159 
025a9bb173 antn*1160       enddo !do jrec=1,nrecsave
f1e8c1d7ee Gael*1161 
                1162         do bj = jtlo,jthi
                1163           do bi = itlo,ithi
                1164             do j = jmin,jmax
                1165               do i = imin,imax
                1166                if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
f8e779c983 antn*1167      &              ( maskc(i,j,1,bi,bj) .NE. 0. )
70fde697ea Gael*1168      &            ) then
f1e8c1d7ee Gael*1169                   anom_psMslaobs(i,j,bi,bj) =
                1170      &                 anom_psMslaobs(i,j,bi,bj) /
                1171      &                 anom_psMslaobs_NUM(i,j,bi,bj)
                1172                   anom_slaobs(i,j,bi,bj) =
                1173      &                 anom_slaobs(i,j,bi,bj) /
                1174      &                 anom_psMslaobs_NUM(i,j,bi,bj)
                1175                else
                1176                   anom_psMslaobs(i,j,bi,bj) = 0. _d 0
                1177                   anom_slaobs(i,j,bi,bj) = 0. _d 0
                1178                endif
                1179               enddo
                1180             enddo
                1181           enddo
                1182         enddo
                1183 
bc9bce397d Gael*1184 c PART 4.11: separate time variable offset
21dbf3b3cd Gael*1185 c ----------------------------------------
                1186 
                1187       slaoffset     = 0. _d 0
                1188       slaoffset_sum = 0. _d 0
bc9bce397d Gael*1189       slaoffset_weight = 0. _d 0
21dbf3b3cd Gael*1190 
0a1233e81a Ou W*1191       if (igen_lsc.GT.0) then
21dbf3b3cd Gael*1192         do bj = jtlo,jthi
                1193           do bi = itlo,ithi
                1194             do j = jmin,jmax
                1195               do i = imin,imax
                1196                if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
f8e779c983 antn*1197      &              ( maskc(i,j,1,bi,bj) .NE. 0. )
21dbf3b3cd Gael*1198      &            ) then
                1199                   slaoffset=slaoffset+RA(i,j,bi,bj)*
                1200      &                      anom_psMslaobs(i,j,bi,bj)
                1201                   slaoffset_sum=slaoffset_sum+RA(i,j,bi,bj)
025a9bb173 antn*1202 C Of interest is the total weight for an average of N
                1203 C indep sample (not the area weighted average of weight)
                1204 C Since slaoffset is itself area weighted, however, the
                1205 C total weight is only approx the simple weight sum:
bc9bce397d Gael*1206                   slaoffset_weight=slaoffset_weight+
                1207      &                      gencost_weight(i,j,bi,bj,igen_lsc)
21dbf3b3cd Gael*1208                endif
                1209               enddo
                1210             enddo
                1211           enddo
                1212         enddo
                1213 
f8e779c983 antn*1214       _GLOBAL_SUM_RL( slaoffset     , myThid )
                1215       _GLOBAL_SUM_RL( slaoffset_weight , myThid )
                1216       _GLOBAL_SUM_RL( slaoffset_sum , myThid )
21dbf3b3cd Gael*1217 
0a1233e81a Ou W*1218        if (igen_gmsl.GT.0) then
                1219         if (slaoffset_sum .EQ. 0.0) then
                1220          if (gencost_outputlevel(igen_gmsl).GT.0) then
                1221          _BEGIN_MASTER( myThid )
                1222          write(msgbuf,'(a)') ' sshv4: slaoffset_sum = zero!'
                1223          call print_message( msgbuf, standardmessageunit,
                1224      &                           SQUEEZE_RIGHT , myThid)
                1225          _END_MASTER( myThid )
                1226          endif
21dbf3b3cd Gael*1227 c        stop   '  ... stopped in cost_ssh.'
0a1233e81a Ou W*1228         else
                1229          if (gencost_outputlevel(igen_gmsl).GT.0) then
                1230          _BEGIN_MASTER( myThid )
                1231          write(msgbuf,'(a,3d22.15)') ' sshv4:slaoffset,sum,weight=',
                1232      &          slaoffset,slaoffset_sum,slaoffset_weight
                1233          call print_message( msgbuf, standardmessageunit,
                1234      &                           SQUEEZE_RIGHT , myThid)
                1235          _END_MASTER( myThid )
                1236          endif
cf137f827f Gael*1237         endif
0a1233e81a Ou W*1238 c      endif igen_gmsl .gt. 0
11c3150c71 Mart*1239        endif
21dbf3b3cd Gael*1240 
                1241 c--   Compute (average) slaoffset
0a1233e81a Ou W*1242        if (slaoffset_sum.NE.0. _d 0) then
                1243         slaoffset = slaoffset / slaoffset_sum
                1244        else
                1245         slaoffset = 0. _d 0
                1246        endif
21dbf3b3cd Gael*1247 
bc9bce397d Gael*1248 c--   Subtract slaoffset from anom_psMslaobs
21dbf3b3cd Gael*1249         do bj = jtlo,jthi
                1250           do bi = itlo,ithi
                1251             do j = jmin,jmax
                1252               do i = imin,imax
                1253                if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
f8e779c983 antn*1254      &              ( maskc(i,j,1,bi,bj) .NE. 0. )
21dbf3b3cd Gael*1255      &            ) then
                1256                   anom_psMslaobs(i,j,bi,bj) =
                1257      &                 anom_psMslaobs(i,j,bi,bj) - slaoffset
                1258                   anom_slaobs(i,j,bi,bj) =
                1259      &                 anom_slaobs(i,j,bi,bj) - slaoffset
                1260                else
                1261                   anom_psMslaobs(i,j,bi,bj) = 0. _d 0
                1262                   anom_slaobs(i,j,bi,bj) = 0. _d 0
                1263                endif
                1264               enddo
                1265             enddo
                1266           enddo
                1267         enddo
                1268 
f1e8c1d7ee Gael*1269 c PART 4.2: smooth anom_psMslaobs in space
                1270 c ----------------------------------------
                1271 
6b881d9117 Gael*1272       if (gencost_outputlevel(igen_lsc).GT.0) then
de57a2ec4b Mart*1273         write(fname4test,'(1a)') 'sladiff_raw'
f8e779c983 antn*1274         CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                1275      &                        anom_psMslaobs, irec, 1, myThid )
de57a2ec4b Mart*1276         write(fname4test,'(1a)') 'slaobs_raw'
f8e779c983 antn*1277         CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                1278      &                        anom_slaobs, irec, 1, myThid )
6b881d9117 Gael*1279       endif
f1e8c1d7ee Gael*1280 
adf79a5d72 Gael*1281 #ifdef ALLOW_SMOOTH
0a1233e81a Ou W*1282       if ( useSMOOTH.AND.(k2_lsc.GT.0) ) then
                1283         call smooth_hetero2d(anom_psMslaobs,maskInC,
3248fb8b5a Gael*1284      &     gencost_posproc_c(k2_lsc,igen_lsc),
f8e779c983 antn*1285      &     gencost_posproc_i(k2_lsc,igen_lsc),myThid)
0a1233e81a Ou W*1286         call smooth_hetero2d(anom_slaobs,maskInC,
3248fb8b5a Gael*1287      &     gencost_posproc_c(k2_lsc,igen_lsc),
f8e779c983 antn*1288      &     gencost_posproc_i(k2_lsc,igen_lsc),myThid)
0a1233e81a Ou W*1289 
                1290         if (gencost_outputlevel(igen_lsc).GT.0) then
de57a2ec4b Mart*1291          write(fname4test,'(1a)') 'sladiff_smooth'
0a1233e81a Ou W*1292          CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                1293      &                         anom_psMslaobs, irec, 1, myThid )
de57a2ec4b Mart*1294          write(fname4test,'(1a)') 'slaobs_smooth'
0a1233e81a Ou W*1295          CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                1296      &                         anom_slaobs, irec, 1, myThid )
                1297         endif
6b881d9117 Gael*1298       endif
0a1233e81a Ou W*1299 #endif
f1e8c1d7ee Gael*1300 
                1301 c PART 4.3: compute cost function term
                1302 c ------------------------------------
                1303 
                1304        do bj = jtlo,jthi
                1305         do bi = itlo,ithi
                1306          do j = jmin,jmax
                1307           do i = imin,imax
d699b51f70 Gael*1308           if ( (gencost_weight(i,j,bi,bj,igen_lsc).GT.0.).AND.
                1309      &         (anom_psMslaobs_MSK(i,j,bi,bj).GT.0.) ) then
                1310             junk = anom_psMslaobs(i,j,bi,bj)
                1311             anom_slaobs_SUB(i,j,bi,bj) = anom_slaobs(i,j,bi,bj)
                1312             anom_psMslaobs_SUB(i,j,bi,bj) = anom_psMslaobs(i,j,bi,bj)
                1313           else
                1314             junk = 0. _d 0
                1315             anom_slaobs_SUB(i,j,bi,bj)= 0. _d 0
                1316             anom_psMslaobs_SUB(i,j,bi,bj)= 0. _d 0
                1317           endif
61bb52ec0d Gael*1318           objf_gencost(bi,bj,igen_lsc) = objf_gencost(bi,bj,igen_lsc)
f1e8c1d7ee Gael*1319      &        + junk*junk*gencost_weight(i,j,bi,bj,igen_lsc)
                1320           if ( (gencost_weight(i,j,bi,bj,igen_lsc).GT.0.).AND.
d699b51f70 Gael*1321      &         (anom_psMslaobs_MSK(i,j,bi,bj).GT.0.) )
61bb52ec0d Gael*1322      &         num_gencost(bi,bj,igen_lsc) =
                1323      &         num_gencost(bi,bj,igen_lsc) + 1. _d 0
f1e8c1d7ee Gael*1324            enddo
                1325          enddo
                1326         enddo
                1327        enddo
                1328 
11c3150c71 Mart*1329        if (gencost_outputlevel(igen_lsc).GT.0) then
de57a2ec4b Mart*1330         write(fname4test,'(1a)') 'sladiff_subsample'
f8e779c983 antn*1331         CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                1332      &                        anom_psMslaobs_SUB, irec, 1, myThid )
de57a2ec4b Mart*1333         write(fname4test,'(1a)') 'slaobs_subsample'
f8e779c983 antn*1334         CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                1335      &                        anom_slaobs_SUB, irec, 1, myThid )
11c3150c71 Mart*1336        endif
                1337 
0a1233e81a Ou W*1338 c     endif igen_lsc .gt. 0
6b881d9117 Gael*1339       endif
d699b51f70 Gael*1340 
bc9bce397d Gael*1341 c PART 4.4: compute cost function term for global mean sea level
                1342 c --------------------------------------------------------------
                1343 
0a1233e81a Ou W*1344       IF ( MASTER_CPU_THREAD(myThid) .AND. ( igen_gmsl.GT.0 )
                1345      &     .AND. ( slaoffset_sum.GT.0.25 _d 0 * globalArea ) ) THEN
bc9bce397d Gael*1346           junk=slaoffset
cf137f827f Gael*1347           junkweight=1. _d 0
bc9bce397d Gael*1348           objf_gencost(1,1,igen_gmsl) = objf_gencost(1,1,igen_gmsl)
                1349      &        + junk*junk*junkweight
                1350           num_gencost(1,1,igen_gmsl) = num_gencost(1,1,igen_gmsl)
                1351      &        + 1. _d 0
                1352 c      print*,'igen_gmsl',igen_gmsl,
                1353       ENDIF
f1e8c1d7ee Gael*1354 
                1355 cgf =======================================================
                1356 cgf PART 5: compute raw SLA cost term
                1357 cgf =======================================================
                1358 
f9a3510c64 Gael*1359 c time associated with this ndaysrec mean
025a9bb173 antn*1360         krec = irec + (nrecsave/2)
f1e8c1d7ee Gael*1361 
101f75e5cd Gael*1362 #ifdef ALLOW_AUTODIFF
025a9bb173 antn*1363         call active_read_xy( fname, etagcm, krec, doglobalread,
f8e779c983 antn*1364      &                       ladinit, eccoiter, myThid,
025a9bb173 antn*1365      &                       gencost_dummy(igen_etagcm) )
101f75e5cd Gael*1366 #else
025a9bb173 antn*1367         CALL READ_REC_XY_RL( fname, etagcm, kRec, 1, myThid )
101f75e5cd Gael*1368 #endif
f1e8c1d7ee Gael*1369 
6b881d9117 Gael*1370         if (.NOT.useEtaMean)
025a9bb173 antn*1371      &  CALL REMOVE_MEAN_RL( 1, etagcm, maskInC, maskInC, rA, drF,
                1372      &                       'etagcm', zeroRL, 1, myThid )
0699f1d2f2 Gael*1373 
14be08a18d Gael*1374        do bj = jtlo,jthi
                1375         do bi = itlo,ithi
                1376          do j = jmin,jmax
                1377           do i = imin,imax
                1378               tpma(i,j,bi,bj) = 0. _d 0
                1379               tpob(i,j,bi,bj) = 0. _d 0
                1380               ersma(i,j,bi,bj) = 0. _d 0
                1381               ersob(i,j,bi,bj) = 0. _d 0
                1382               gfoma(i,j,bi,bj) = 0. _d 0
                1383               gfoob(i,j,bi,bj) = 0. _d 0
                1384           enddo
                1385          enddo
                1386         enddo
                1387        enddo
                1388 
f09238ab8f Gael*1389         if (using_tpj)
025a9bb173 antn*1390      &  call cost_sla_read( tpfi, tpsta, tpper, use_mon,
f09238ab8f Gael*1391      &                zeroRL, zeroRL,
                1392      &                tpob, tpma,
f8e779c983 antn*1393      &                krec, myThid )
f09238ab8f Gael*1394         if (using_ers)
025a9bb173 antn*1395      &  call cost_sla_read( ersfi, erssta, ersper, use_mon,
f09238ab8f Gael*1396      &                zeroRL, zeroRL,
                1397      &                ersob, ersma,
f8e779c983 antn*1398      &                krec, myThid )
f09238ab8f Gael*1399         if (using_gfo)
025a9bb173 antn*1400      &  call cost_sla_read( gfofi, gfosta, gfoper, use_mon,
f09238ab8f Gael*1401      &                zeroRL, zeroRL,
                1402      &                gfoob, gfoma,
f8e779c983 antn*1403      &                krec, myThid )
f1e8c1d7ee Gael*1404 
                1405        do bj = jtlo,jthi
                1406         do bi = itlo,ithi
                1407          do j = jmin,jmax
                1408           do i = imin,imax
                1409               anom_psMtpobs(i,j,bi,bj)  = 0. _d 0
                1410               anom_psMersobs(i,j,bi,bj) = 0. _d 0
                1411               anom_psMgfoobs(i,j,bi,bj) = 0. _d 0
                1412               anom_tpobs(i,j,bi,bj)  = 0. _d 0
                1413               anom_ersobs(i,j,bi,bj) = 0. _d 0
                1414               anom_gfoobs(i,j,bi,bj) = 0. _d 0
                1415           enddo
                1416          enddo
                1417         enddo
                1418        enddo
                1419 
                1420        do bj = jtlo,jthi
                1421         do bi = itlo,ithi
                1422          do j = jmin,jmax
                1423           do i = imin,imax
f09238ab8f Gael*1424       if ( tpma(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj).NE.0. )
f1e8c1d7ee Gael*1425      & then
                1426          anom_psMtpobs(i,j,bi,bj)=
025a9bb173 antn*1427      &       etagcm(i,j,bi,bj) - tpob(i,j,bi,bj)
f9a3510c64 Gael*1428      &       - mean_psMslaobs(i,j,bi,bj) - slaoffset
                1429      &       - anom_psMslaobs(i,j,bi,bj)
f09238ab8f Gael*1430          anom_tpobs(i,j,bi,bj)=tpob(i,j,bi,bj) - slaoffset
f9a3510c64 Gael*1431      &       - anom_slaobs(i,j,bi,bj)
f1e8c1d7ee Gael*1432       endif
f09238ab8f Gael*1433       if ( ersma(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj).NE.0. )
f1e8c1d7ee Gael*1434      & then
                1435          anom_psMersobs(i,j,bi,bj)=
025a9bb173 antn*1436      &       etagcm(i,j,bi,bj) - ersob(i,j,bi,bj)
f9a3510c64 Gael*1437      &       - mean_psMslaobs(i,j,bi,bj) - slaoffset
                1438      &       - anom_psMslaobs(i,j,bi,bj)
f09238ab8f Gael*1439          anom_ersobs(i,j,bi,bj)=ersob(i,j,bi,bj) - slaoffset
f9a3510c64 Gael*1440      &       - anom_slaobs(i,j,bi,bj)
f1e8c1d7ee Gael*1441       endif
f09238ab8f Gael*1442       if ( gfoma(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj).NE.0. )
f1e8c1d7ee Gael*1443      & then
                1444          anom_psMgfoobs(i,j,bi,bj)=
025a9bb173 antn*1445      &       etagcm(i,j,bi,bj) - gfoob(i,j,bi,bj)
f9a3510c64 Gael*1446      &       - mean_psMslaobs(i,j,bi,bj) - slaoffset
                1447      &       - anom_psMslaobs(i,j,bi,bj)
f09238ab8f Gael*1448          anom_gfoobs(i,j,bi,bj)=gfoob(i,j,bi,bj) - slaoffset
f9a3510c64 Gael*1449      &       - anom_slaobs(i,j,bi,bj)
f1e8c1d7ee Gael*1450       endif
                1451           enddo
                1452          enddo
                1453         enddo
                1454        enddo
                1455 
0a1233e81a Ou W*1456       if ( igen_tp .GT. 0 ) then
11c3150c71 Mart*1457        if (gencost_outputlevel(igen_tp).GT.0) then
de57a2ec4b Mart*1458         write(fname4test,'(1a)') 'sladiff_tp_raw'
f8e779c983 antn*1459         CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                1460      &                        anom_psMtpobs, irec, 1, myThid )
de57a2ec4b Mart*1461         write(fname4test,'(1a)') 'slaobs_tp_raw'
f8e779c983 antn*1462         CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                1463      &                        anom_tpobs, irec, 1, myThid )
11c3150c71 Mart*1464        endif
6b881d9117 Gael*1465       endif
                1466 
0a1233e81a Ou W*1467       if ( igen_ers .GT. 0 ) then
11c3150c71 Mart*1468        if (gencost_outputlevel(igen_ers).GT.0) then
de57a2ec4b Mart*1469         write(fname4test,'(1a)') 'sladiff_ers_raw'
f8e779c983 antn*1470         CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                1471      &                        anom_psMersobs, irec, 1, myThid )
de57a2ec4b Mart*1472         write(fname4test,'(1a)') 'slaobs_ers_raw'
f8e779c983 antn*1473         CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                1474      &                        anom_ersobs, irec, 1, myThid )
11c3150c71 Mart*1475        endif
6b881d9117 Gael*1476       endif
                1477 
0a1233e81a Ou W*1478       if ( igen_gfo .GT. 0 ) then
11c3150c71 Mart*1479        if (gencost_outputlevel(igen_gfo).GT.0) then
de57a2ec4b Mart*1480         write(fname4test,'(1a)') 'sladiff_gfo_raw'
f8e779c983 antn*1481         CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                1482      &                        anom_psMgfoobs, irec, 1, myThid )
de57a2ec4b Mart*1483         write(fname4test,'(1a)') 'slaobs_gfo_raw'
f8e779c983 antn*1484         CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                1485      &                        anom_gfoobs, irec, 1, myThid )
11c3150c71 Mart*1486        endif
6b881d9117 Gael*1487       endif
0a1233e81a Ou W*1488       if ( igen_tp .GT. 0 ) then
f1e8c1d7ee Gael*1489        do bj = jtlo,jthi
                1490         do bi = itlo,ithi
                1491          do j = jmin,jmax
                1492           do i = imin,imax
                1493 c--             The array psobs contains SSH anomalies.
52e9fa3441 Gael*1494                 junkweight = mean_psMslaobs_MSK(i,j,bi,bj)
f1e8c1d7ee Gael*1495      &                      *gencost_weight(i,j,bi,bj,igen_tp)
f09238ab8f Gael*1496      &                      *tpma(i,j,bi,bj)
f1e8c1d7ee Gael*1497                 junk = anom_psMtpobs(i,j,bi,bj)
61bb52ec0d Gael*1498                 objf_gencost(bi,bj,igen_tp) =
                1499      &              objf_gencost(bi,bj,igen_tp)+junk*junk*junkweight
0a1233e81a Ou W*1500                 if ( junkweight .NE. 0. )
61bb52ec0d Gael*1501      &              num_gencost(bi,bj,igen_tp) =
                1502      &              num_gencost(bi,bj,igen_tp) + 1. _d 0
11c3150c71 Mart*1503           enddo
                1504          enddo
                1505         enddo
                1506        enddo
                1507       endif
0a1233e81a Ou W*1508       if ( igen_ers .GT. 0 ) then
11c3150c71 Mart*1509        do bj = jtlo,jthi
                1510         do bi = itlo,ithi
                1511          do j = jmin,jmax
                1512           do i = imin,imax
f1e8c1d7ee Gael*1513 c--             The array ersobs contains SSH anomalies.
52e9fa3441 Gael*1514                 junkweight = mean_psMslaobs_MSK(i,j,bi,bj)
f1e8c1d7ee Gael*1515      &                      *gencost_weight(i,j,bi,bj,igen_ers)
f09238ab8f Gael*1516      &                      *ersma(i,j,bi,bj)
f1e8c1d7ee Gael*1517                 junk = anom_psMersobs(i,j,bi,bj)
61bb52ec0d Gael*1518                 objf_gencost(bi,bj,igen_ers) =
                1519      &               objf_gencost(bi,bj,igen_ers)+junk*junk*junkweight
0a1233e81a Ou W*1520                 if ( junkweight .NE. 0. )
61bb52ec0d Gael*1521      &              num_gencost(bi,bj,igen_ers) =
                1522      &              num_gencost(bi,bj,igen_ers) + 1. _d 0
11c3150c71 Mart*1523           enddo
                1524          enddo
                1525         enddo
                1526        enddo
                1527       endif
0a1233e81a Ou W*1528       if ( igen_gfo .GT. 0 ) then
11c3150c71 Mart*1529        do bj = jtlo,jthi
                1530         do bi = itlo,ithi
                1531          do j = jmin,jmax
                1532           do i = imin,imax
f1e8c1d7ee Gael*1533 c--             The array gfoobs contains SSH anomalies.
52e9fa3441 Gael*1534                 junkweight = mean_psMslaobs_MSK(i,j,bi,bj)
f1e8c1d7ee Gael*1535      &                      *gencost_weight(i,j,bi,bj,igen_gfo)
f09238ab8f Gael*1536      &                      *gfoma(i,j,bi,bj)
f1e8c1d7ee Gael*1537                 junk = anom_psMgfoobs(i,j,bi,bj)
61bb52ec0d Gael*1538                 objf_gencost(bi,bj,igen_gfo) =
                1539      &              objf_gencost(bi,bj,igen_gfo)+junk*junk*junkweight
0a1233e81a Ou W*1540                 if ( junkweight .NE. 0. )
61bb52ec0d Gael*1541      &              num_gencost(bi,bj,igen_gfo) =
                1542      &              num_gencost(bi,bj,igen_gfo) + 1. _d 0
11c3150c71 Mart*1543           enddo
f1e8c1d7ee Gael*1544          enddo
                1545         enddo
                1546        enddo
11c3150c71 Mart*1547       endif
f1e8c1d7ee Gael*1548 
025a9bb173 antn*1549       enddo !do irec = 1, nrec_total-nrecsave+1
f1e8c1d7ee Gael*1550 
                1551 #endif /* ifdef ALLOW_GENCOST_CONTRIBUTION */
                1552 
f8e779c983 antn*1553       RETURN
                1554       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!