| 1 | mmazloff | 1.1 | C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_ssh_new.F,v 1.5 2009/04/28 18:13:28 jmc Exp $ | 
| 2 |  |  | C $Name:  $ | 
| 3 |  |  |  | 
| 4 |  |  | #include "COST_CPPOPTIONS.h" | 
| 5 |  |  |  | 
| 6 |  |  |  | 
| 7 |  |  | subroutine cost_ssh_new( | 
| 8 |  |  | I                     myiter, | 
| 9 |  |  | I                     mytime, | 
| 10 |  |  | I                     mythid | 
| 11 |  |  | &                   ) | 
| 12 |  |  |  | 
| 13 |  |  | c     ================================================================== | 
| 14 |  |  | c     SUBROUTINE cost_ssh | 
| 15 |  |  | c     ================================================================== | 
| 16 |  |  | c | 
| 17 |  |  | c     o Evaluate cost function contribution of sea surface height. | 
| 18 |  |  | c       using of geoid error covariances requires regular model grid | 
| 19 |  |  | c | 
| 20 |  |  | c     started: Detlef Stammer, Ralf Giering Jul-1996 | 
| 21 |  |  | c | 
| 22 |  |  | c     changed: Christian Eckert eckert@mit.edu 25-Feb-2000 | 
| 23 |  |  | c | 
| 24 |  |  | c              - Restructured the code in order to create a package | 
| 25 |  |  | c                for the MITgcmUV. | 
| 26 |  |  | c | 
| 27 |  |  | c     changed: Ralf Giering Ralf.Giering@FastOpt.de 12-Jun-2001 | 
| 28 |  |  | c | 
| 29 |  |  | c              - totally rewrite for parallel processing | 
| 30 |  |  | c | 
| 31 |  |  | c     ================================================================== | 
| 32 |  |  | c     SUBROUTINE cost_ssh | 
| 33 |  |  | c     ================================================================== | 
| 34 |  |  |  | 
| 35 |  |  | implicit none | 
| 36 |  |  |  | 
| 37 |  |  | c     == global variables == | 
| 38 |  |  |  | 
| 39 |  |  | #include "EEPARAMS.h" | 
| 40 |  |  | #include "SIZE.h" | 
| 41 |  |  | #include "PARAMS.h" | 
| 42 |  |  |  | 
| 43 |  |  | #include "ecco_cost.h" | 
| 44 |  |  | #include "ctrl.h" | 
| 45 |  |  | #include "ctrl_dummy.h" | 
| 46 |  |  | #include "optim.h" | 
| 47 |  |  | #include "DYNVARS.h" | 
| 48 |  |  | #ifdef ALLOW_PROFILES | 
| 49 |  |  | #include "profiles.h" | 
| 50 |  |  | #endif | 
| 51 |  |  |  | 
| 52 |  |  | c     == routine arguments == | 
| 53 |  |  |  | 
| 54 |  |  | integer myiter | 
| 55 |  |  | _RL     mytime | 
| 56 |  |  | integer mythid | 
| 57 |  |  |  | 
| 58 |  |  | #ifdef ALLOW_SSH_COST_CONTRIBUTION | 
| 59 |  |  | c     == local variables == | 
| 60 |  |  |  | 
| 61 |  |  | integer bi,bj | 
| 62 |  |  | integer i,j | 
| 63 |  |  | integer itlo,ithi | 
| 64 |  |  | integer jtlo,jthi | 
| 65 |  |  | integer jmin,jmax | 
| 66 |  |  | integer imin,imax | 
| 67 |  |  | integer irec | 
| 68 |  |  | integer ilps | 
| 69 |  |  | integer gwunit | 
| 70 |  |  |  | 
| 71 |  |  | logical doglobalread | 
| 72 |  |  | logical ladinit | 
| 73 |  |  |  | 
| 74 |  |  | _RL offset | 
| 75 |  |  | _RL costmean | 
| 76 |  |  | _RL offset_sum | 
| 77 |  |  | _RL psmean    ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 78 |  |  | _RL psmeantp  ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 79 |  |  | _RL psmeaners ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 80 |  |  | _RL psmeangfo ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 81 |  |  | _RL sumtp  ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 82 |  |  | _RL sumers ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 83 |  |  | _RL sumgfo ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 84 |  |  |  | 
| 85 |  |  | _RL wwwtp1  ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 86 |  |  | _RL wwwers1 ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 87 |  |  | _RL wwwgfo1 ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 88 |  |  | _RL wwwtp2  ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 89 |  |  | _RL wwwers2 ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 90 |  |  | _RL wwwgfo2 ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 91 |  |  | _RL junk,diff | 
| 92 |  |  | CMM( | 
| 93 |  |  | _RL obsmeantp  ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 94 |  |  | _RL obsmeaners ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 95 |  |  | _RL obsmeangfo ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 96 |  |  | _RL tpmean_loc ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy ) | 
| 97 |  |  | _RL mean_slaobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) | 
| 98 |  |  | _RL mean_slaobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) | 
| 99 |  |  | _RL mean_psMssh_all(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) | 
| 100 |  |  | _RL mean_psMssh_all_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) | 
| 101 |  |  | _RL mean_psMssh_all_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) | 
| 102 |  |  | _RL diagnosfld(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) | 
| 103 |  |  | _RL sumc, num_hmean_mm | 
| 104 |  |  | _RL factor, spval | 
| 105 |  |  | integer IDOT, k | 
| 106 |  |  | character*(80) fname4test | 
| 107 |  |  | CMM) | 
| 108 |  |  | character*(80) fname | 
| 109 |  |  | character*(MAX_LEN_MBUF) msgbuf | 
| 110 |  |  |  | 
| 111 |  |  | c     == external functions == | 
| 112 |  |  |  | 
| 113 |  |  | integer  ilnblnk | 
| 114 |  |  | external ilnblnk | 
| 115 |  |  |  | 
| 116 |  |  | c     == end of interface == | 
| 117 |  |  |  | 
| 118 |  |  | jtlo = mybylo(mythid) | 
| 119 |  |  | jthi = mybyhi(mythid) | 
| 120 |  |  | itlo = mybxlo(mythid) | 
| 121 |  |  | ithi = mybxhi(mythid) | 
| 122 |  |  | jmin = 1 | 
| 123 |  |  | jmax = sny | 
| 124 |  |  | imin = 1 | 
| 125 |  |  | imax = snx | 
| 126 |  |  |  | 
| 127 |  |  | c--   Initialise local variables. | 
| 128 |  |  | costmean   = 0. _d 0 | 
| 129 |  |  |  | 
| 130 |  |  | c--   First, read tiled data. | 
| 131 |  |  | doglobalread = .false. | 
| 132 |  |  | ladinit      = .false. | 
| 133 |  |  |  | 
| 134 |  |  | write(fname(1:80),'(80a)') ' ' | 
| 135 |  |  | ilps=ilnblnk( psbarfile ) | 
| 136 |  |  | write(fname(1:80),'(2a,i10.10)') | 
| 137 |  |  | &     psbarfile(1:ilps),'.',optimcycle | 
| 138 |  |  |  | 
| 139 |  |  | c--   ============ | 
| 140 |  |  | c--   Mean values. | 
| 141 |  |  | c--   ============ | 
| 142 |  |  |  | 
| 143 |  |  | do bj = jtlo,jthi | 
| 144 |  |  | do bi = itlo,ithi | 
| 145 |  |  | do j = jmin,jmax | 
| 146 |  |  | do i = imin,imax | 
| 147 |  |  | psmean(i,j,bi,bj)    = 0. _d 0 | 
| 148 |  |  | psmeantp(i,j,bi,bj)  = 0. _d 0 | 
| 149 |  |  | psmeaners(i,j,bi,bj) = 0. _d 0 | 
| 150 |  |  | psmeangfo(i,j,bi,bj) = 0. _d 0 | 
| 151 |  |  | CMM(  initialize | 
| 152 |  |  | obsmeantp(i,j,bi,bj) = 0. _d 0 | 
| 153 |  |  | obsmeaners(i,j,bi,bj) = 0. _d 0 | 
| 154 |  |  | obsmeangfo(i,j,bi,bj) = 0. _d 0 | 
| 155 |  |  | mean_slaobs(i,j,bi,bj)  = 0. _d 0 | 
| 156 |  |  | mean_slaobs_NUM(i,j,bi,bj)  = 0. _d 0 | 
| 157 |  |  | mean_psMssh_all(i,j,bi,bj) = 0. _d 0 | 
| 158 |  |  | mean_psMssh_all_NUM(i,j,bi,bj) = 0. _d 0 | 
| 159 |  |  | mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0 | 
| 160 |  |  | diagnosfld(i,j,bi,bj)  = 0. _d 0 | 
| 161 |  |  | tpmean_loc(i,j,bi,bj) = 0. _d 0 | 
| 162 |  |  | CMM) | 
| 163 |  |  | sumtp(i,j,bi,bj)  = 0. _d 0 | 
| 164 |  |  | sumers(i,j,bi,bj) = 0. _d 0 | 
| 165 |  |  | sumgfo(i,j,bi,bj) = 0. _d 0 | 
| 166 |  |  | wwwtp1(i,j,bi,bj)  = 0. _d 0 | 
| 167 |  |  | wwwers1(i,j,bi,bj) = 0. _d 0 | 
| 168 |  |  | wwwgfo1(i,j,bi,bj) = 0. _d 0 | 
| 169 |  |  | wwwtp2(i,j,bi,bj)  = 0. _d 0 | 
| 170 |  |  | wwwers2(i,j,bi,bj) = 0. _d 0 | 
| 171 |  |  | wwwgfo2(i,j,bi,bj) = 0. _d 0 | 
| 172 |  |  | enddo | 
| 173 |  |  | enddo | 
| 174 |  |  | enddo | 
| 175 |  |  | enddo | 
| 176 |  |  |  | 
| 177 |  |  | c--   Read mean field and generate mask | 
| 178 |  |  | call cost_ReadTopexMean( mythid ) | 
| 179 |  |  |  | 
| 180 |  |  | CMM(  HAD TO MOVE OFFSET CALC TO GET IT INTO RECORD LOOP | 
| 181 |  |  | c--   Compute and remove offset for current tile and sum over all | 
| 182 |  |  | c--   tiles of this instance. | 
| 183 |  |  | offset     = 0. _d 0 | 
| 184 |  |  | offset_sum = 0. _d 0 | 
| 185 |  |  | CMM) | 
| 186 |  |  | c--   Loop over records for the first time. | 
| 187 |  |  | do irec = 1, ndaysrec | 
| 188 |  |  |  | 
| 189 |  |  | c--     Compute the mean over all psbar records. | 
| 190 |  |  | call active_read_xy( fname, psbar, irec, doglobalread, | 
| 191 |  |  | &                       ladinit, optimcycle, mythid, | 
| 192 |  |  | &                       xx_psbar_mean_dummy ) | 
| 193 |  |  |  | 
| 194 |  |  | #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION | 
| 195 |  |  | cph these if would be more efficient, but need recomputations | 
| 196 |  |  | c        if ( tpTimeMask(irec) .NE. 0. ) | 
| 197 |  |  | call cost_readtopex( irec, mythid ) | 
| 198 |  |  | #endif | 
| 199 |  |  | #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION | 
| 200 |  |  | cph these if would be more efficient, but need recomputations | 
| 201 |  |  | c        if ( ersTimeMask(irec) .NE. 0. ) | 
| 202 |  |  | call cost_readers( irec, mythid ) | 
| 203 |  |  | #endif | 
| 204 |  |  | #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION | 
| 205 |  |  | cph these if would be more efficient, but need recomputations | 
| 206 |  |  | c        if ( gfoTimeMask(irec) .NE. 0. ) | 
| 207 |  |  | call cost_readgfo( irec, mythid ) | 
| 208 |  |  | #endif | 
| 209 |  |  |  | 
| 210 |  |  | do bj = jtlo,jthi | 
| 211 |  |  | do bi = itlo,ithi | 
| 212 |  |  | do j = jmin,jmax | 
| 213 |  |  | do i = imin,imax | 
| 214 |  |  | psmean(i,j,bi,bj) = psmean(i,j,bi,bj) + | 
| 215 |  |  | &                psbar(i,j,bi,bj) / float(ndaysrec) | 
| 216 |  |  | #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION | 
| 217 |  |  | if ( tpmask(i,j,bi,bj)*wtp(i,j,bi,bj) | 
| 218 |  |  | &               *tpTimeMask(irec) .NE. 0. ) then | 
| 219 |  |  | psmeantp(i,j,bi,bj) = psmeantp(i,j,bi,bj) + | 
| 220 |  |  | &                 psbar(i,j,bi,bj) | 
| 221 |  |  | obsmeantp(i,j,bi,bj) = obsmeantp(i,j,bi,bj) + | 
| 222 |  |  | &                 tpobs(i,j,bi,bj) | 
| 223 |  |  | sumtp(i,j,bi,bj) = sumtp(i,j,bi,bj) + 1. _d 0 | 
| 224 |  |  | endif | 
| 225 |  |  | #endif | 
| 226 |  |  | #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION | 
| 227 |  |  | if ( ersmask(i,j,bi,bj)*wers(i,j,bi,bj) | 
| 228 |  |  | &               *ersTimeMask(irec) .NE. 0. ) then | 
| 229 |  |  | psmeaners(i,j,bi,bj) = psmeaners(i,j,bi,bj) + | 
| 230 |  |  | &                 psbar(i,j,bi,bj) | 
| 231 |  |  | obsmeaners(i,j,bi,bj) = obsmeaners(i,j,bi,bj) + | 
| 232 |  |  | &                 ersobs(i,j,bi,bj) | 
| 233 |  |  | sumers(i,j,bi,bj) = sumers(i,j,bi,bj) + 1. _d 0 | 
| 234 |  |  | endif | 
| 235 |  |  | #endif | 
| 236 |  |  | #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION | 
| 237 |  |  | if ( gfomask(i,j,bi,bj)*wgfo(i,j,bi,bj) | 
| 238 |  |  | &               *gfoTimeMask(irec) .NE. 0. ) then | 
| 239 |  |  | psmeangfo(i,j,bi,bj) = psmeangfo(i,j,bi,bj) + | 
| 240 |  |  | &                 psbar(i,j,bi,bj) | 
| 241 |  |  | obsmeangfo(i,j,bi,bj) = obsmeangfo(i,j,bi,bj) + | 
| 242 |  |  | &                 gfoobs(i,j,bi,bj) | 
| 243 |  |  | sumgfo(i,j,bi,bj) = sumgfo(i,j,bi,bj) + 1. _d 0 | 
| 244 |  |  | endif | 
| 245 |  |  | #endif | 
| 246 |  |  | enddo | 
| 247 |  |  | enddo | 
| 248 |  |  | enddo | 
| 249 |  |  | enddo | 
| 250 |  |  |  | 
| 251 |  |  | c--   END loop over records for the first time. | 
| 252 |  |  | enddo | 
| 253 |  |  |  | 
| 254 |  |  | do bj = jtlo,jthi | 
| 255 |  |  | do bi = itlo,ithi | 
| 256 |  |  | do j = jmin,jmax | 
| 257 |  |  | do i = imin,imax | 
| 258 |  |  | CMM(  average up all SLA | 
| 259 |  |  | mean_slaobs(i,j,bi,bj) = obsmeantp(i,j,bi,bj) + | 
| 260 |  |  | &               obsmeaners(i,j,bi,bj) + obsmeangfo(i,j,bi,bj) | 
| 261 |  |  | mean_slaobs_NUM(i,j,bi,bj) = sumtp(i,j,bi,bj) + | 
| 262 |  |  | &               sumers(i,j,bi,bj) + sumgfo(i,j,bi,bj) | 
| 263 |  |  | if ( ( mean_slaobs_NUM(i,j,bi,bj) .NE. 0. ).AND. | 
| 264 |  |  | &              ( tpmeanmask(i,j,bi,bj) .NE. 0. ) ) then | 
| 265 |  |  | mean_slaobs(i,j,bi,bj) = mean_slaobs(i,j,bi,bj) / | 
| 266 |  |  | &                 mean_slaobs_NUM(i,j,bi,bj) | 
| 267 |  |  | else | 
| 268 |  |  | mean_slaobs(i,j,bi,bj) = 0. _d 0 | 
| 269 |  |  | endif | 
| 270 |  |  | #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION | 
| 271 |  |  | if ( sumtp(i,j,bi,bj) .NE. 0. ) then | 
| 272 |  |  | mean_psMssh_all(i,j,bi,bj) = | 
| 273 |  |  | &                 mean_psMssh_all(i,j,bi,bj) + | 
| 274 |  |  | &                 psmeantp(i,j,bi,bj) - obsmeantp(i,j,bi,bj) | 
| 275 |  |  | mean_psMssh_all_NUM(i,j,bi,bj) = | 
| 276 |  |  | &                 mean_psMssh_all_NUM(i,j,bi,bj) + | 
| 277 |  |  | &                 sumtp(i,j,bi,bj) | 
| 278 |  |  | psmeantp(i,j,bi,bj) = psmeantp(i,j,bi,bj) / | 
| 279 |  |  | &                 sumtp(i,j,bi,bj) | 
| 280 |  |  | obsmeantp(i,j,bi,bj) = obsmeantp(i,j,bi,bj) / | 
| 281 |  |  | &                 sumtp(i,j,bi,bj) | 
| 282 |  |  | wwwtp1(i,j,bi,bj) = 1. _d 0 | 
| 283 |  |  | endif | 
| 284 |  |  | #endif | 
| 285 |  |  | #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION | 
| 286 |  |  | if ( sumers(i,j,bi,bj) .NE. 0. ) then | 
| 287 |  |  | mean_psMssh_all(i,j,bi,bj) = | 
| 288 |  |  | &                 mean_psMssh_all(i,j,bi,bj) + | 
| 289 |  |  | &                 psmeaners(i,j,bi,bj) - obsmeaners(i,j,bi,bj) | 
| 290 |  |  | mean_psMssh_all_NUM(i,j,bi,bj) = | 
| 291 |  |  | &                 mean_psMssh_all_NUM(i,j,bi,bj) + | 
| 292 |  |  | &                 sumers(i,j,bi,bj) | 
| 293 |  |  | psmeaners(i,j,bi,bj) = psmeaners(i,j,bi,bj) / | 
| 294 |  |  | &                 sumers(i,j,bi,bj) | 
| 295 |  |  | obsmeaners(i,j,bi,bj) = obsmeaners(i,j,bi,bj) / | 
| 296 |  |  | &                 sumers(i,j,bi,bj) | 
| 297 |  |  | wwwers1(i,j,bi,bj) = 1. _d 0 | 
| 298 |  |  | endif | 
| 299 |  |  | #endif | 
| 300 |  |  | #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION | 
| 301 |  |  | if ( sumgfo(i,j,bi,bj) .NE. 0. ) then | 
| 302 |  |  | mean_psMssh_all(i,j,bi,bj) = | 
| 303 |  |  | &                 mean_psMssh_all(i,j,bi,bj) + | 
| 304 |  |  | &                 psmeangfo(i,j,bi,bj)-obsmeangfo(i,j,bi,bj) | 
| 305 |  |  | mean_psMssh_all_NUM(i,j,bi,bj) = | 
| 306 |  |  | &                 mean_psMssh_all_NUM(i,j,bi,bj) + | 
| 307 |  |  | &                 sumgfo(i,j,bi,bj) | 
| 308 |  |  | psmeangfo(i,j,bi,bj) = psmeangfo(i,j,bi,bj) / | 
| 309 |  |  | &                 sumgfo(i,j,bi,bj) | 
| 310 |  |  | obsmeangfo(i,j,bi,bj) = obsmeangfo(i,j,bi,bj) / | 
| 311 |  |  | &                 sumgfo(i,j,bi,bj) | 
| 312 |  |  | wwwgfo1(i,j,bi,bj) = 1. _d 0 | 
| 313 |  |  | endif | 
| 314 |  |  | #endif | 
| 315 |  |  | if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .NE. 0. ).AND. | 
| 316 |  |  | &              ( tpmeanmask(i,j,bi,bj) .NE. 0. ) ) then | 
| 317 |  |  | mean_psMssh_all(i,j,bi,bj) = | 
| 318 |  |  | &                 mean_psMssh_all(i,j,bi,bj) / | 
| 319 |  |  | &                 mean_psMssh_all_NUM(i,j,bi,bj)-tpmean(i,j,bi,bj) | 
| 320 |  |  | &                 -mean_slaobs(i,j,bi,bj) | 
| 321 |  |  | mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0 | 
| 322 |  |  | offset=offset+mean_psMssh_all(i,j,bi,bj)* | 
| 323 |  |  | &                 mean_psMssh_all_NUM(i,j,bi,bj) | 
| 324 |  |  | offset_sum=offset_sum+mean_psMssh_all_NUM(i,j,bi,bj) | 
| 325 |  |  | else | 
| 326 |  |  | mean_psMssh_all(i,j,bi,bj) = 0. _d 0 | 
| 327 |  |  | mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0 | 
| 328 |  |  | endif | 
| 329 |  |  | CMM(  OFFSET IS < DOT + obs - model > | 
| 330 |  |  | enddo | 
| 331 |  |  | enddo | 
| 332 |  |  | enddo | 
| 333 |  |  | enddo | 
| 334 |  |  | CMM( output SLA mean for referencing | 
| 335 |  |  | write(fname4test(1:80),'(1a)') 'sla2mdt_raw' | 
| 336 |  |  | call mdswritefield(fname4test,32,.false.,'RL', | 
| 337 |  |  | & 1,mean_slaobs,1,1,mythid) | 
| 338 |  |  | CMM) | 
| 339 |  |  |  | 
| 340 |  |  | c--   Do a global summation. | 
| 341 |  |  | _GLOBAL_SUM_RL( offset     , mythid ) | 
| 342 |  |  | _GLOBAL_SUM_RL( offset_sum , mythid ) | 
| 343 |  |  |  | 
| 344 |  |  | if (offset_sum .eq. 0.0) then | 
| 345 |  |  | _BEGIN_MASTER( mythid ) | 
| 346 |  |  | write(msgbuf,'(a)') ' cost_ssh: offset_sum = zero!' | 
| 347 |  |  | call print_message( msgbuf, standardmessageunit, | 
| 348 |  |  | &                          SQUEEZE_RIGHT , mythid) | 
| 349 |  |  | _END_MASTER( mythid ) | 
| 350 |  |  | stop   '  ... stopped in cost_ssh.' | 
| 351 |  |  | else | 
| 352 |  |  | _BEGIN_MASTER( mythid ) | 
| 353 |  |  | write(msgbuf,'(a,d22.15)') | 
| 354 |  |  | &        ' cost_ssh: offset_sum = ',offset_sum | 
| 355 |  |  | call print_message( msgbuf, standardmessageunit, | 
| 356 |  |  | &                          SQUEEZE_RIGHT , mythid) | 
| 357 |  |  | write(msgbuf,'(a,d22.15)') | 
| 358 |  |  | &          ' cost_ssh: offset_tot = ',offset | 
| 359 |  |  | call print_message( msgbuf, standardmessageunit, | 
| 360 |  |  | &                          SQUEEZE_RIGHT , mythid) | 
| 361 |  |  | _END_MASTER( mythid ) | 
| 362 |  |  | endif | 
| 363 |  |  |  | 
| 364 |  |  | offset = offset / offset_sum | 
| 365 |  |  |  | 
| 366 |  |  | #ifdef ALLOW_PROFILES | 
| 367 |  |  | do bj = jtlo,jthi | 
| 368 |  |  | do bi = itlo,ithi | 
| 369 |  |  | do j = 1,sny | 
| 370 |  |  | do i = 1,snx | 
| 371 |  |  | prof_etan_mean(i,j,bi,bj)=offset+psmean(i,j,bi,bj) | 
| 372 |  |  | enddo | 
| 373 |  |  | enddo | 
| 374 |  |  | enddo | 
| 375 |  |  | enddo | 
| 376 |  |  | _EXCH_XY_RL( prof_etan_mean, mythid ) | 
| 377 |  |  | #endif | 
| 378 |  |  |  | 
| 379 |  |  | #ifdef ALLOW_SSH_MEAN_COST_CONTRIBUTION | 
| 380 |  |  | #ifndef ALLOW_SSH_TOT | 
| 381 |  |  | c--   ========== | 
| 382 |  |  | c--      Mean | 
| 383 |  |  | c--   ========== | 
| 384 |  |  | c--   compute mean ssh difference cost contribution | 
| 385 |  |  | c--      Mean | 
| 386 |  |  | CMM   mean_psMssh_all  = <model - obs> - (DOT + <obs>) | 
| 387 |  |  | do bj = jtlo,jthi | 
| 388 |  |  | do bi = itlo,ithi | 
| 389 |  |  | do j = jmin,jmax | 
| 390 |  |  | do i = imin,imax | 
| 391 |  |  | if ( psmean(i,j,bi,bj) .ne. 0. ) then | 
| 392 |  |  | diagnosfld(i,j,bi,bj) = (mean_psMssh_all(i,j,bi,bj) -  offset) | 
| 393 |  |  | &     *mean_psMssh_all_MSK(i,j,bi,bj) | 
| 394 |  |  | diff = mean_psMssh_all(i,j,bi,bj) -  offset | 
| 395 |  |  | sumc = sumc + diff*diff | 
| 396 |  |  | &          * wp(i,j,bi,bj)*mean_psMssh_all_MSK(i,j,bi,bj) | 
| 397 |  |  | if ( wp(i,j,bi,bj)*mean_psMssh_all_MSK(i,j,bi,bj) .ne. 0. ) | 
| 398 |  |  | &             num_hmean_mm = num_hmean_mm + 1. _d 0 | 
| 399 |  |  | CMM         diagnosfld(i,j,bi,bj) = diff*diff | 
| 400 |  |  | CMM     &      * wp(i,j,bi,bj)*mean_psMssh_all_MSK(i,j,bi,bj) | 
| 401 |  |  | endif | 
| 402 |  |  | enddo | 
| 403 |  |  | enddo | 
| 404 |  |  | enddo | 
| 405 |  |  | enddo | 
| 406 |  |  |  | 
| 407 |  |  | CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean1', ' ', diagnosfld, | 
| 408 |  |  | &                           optimcycle, mythid ) | 
| 409 |  |  |  | 
| 410 |  |  | objf_hmean = objf_hmean + sumc | 
| 411 |  |  | num_hmean = num_hmean + num_hmean_mm | 
| 412 |  |  | CMM) | 
| 413 |  |  | CMM(((((HERE ADD OTHER DOT FOR DOT PROJECT | 
| 414 |  |  | CMM     THIS IS ALL HARDCODED! | 
| 415 |  |  | CMM     DO THIS for 4 DOTS -- first read the field | 
| 416 |  |  | CMM first get orig cost | 
| 417 |  |  | c--   Do a global summation. | 
| 418 |  |  | _GLOBAL_SUM_RL( sumc         , mythid ) | 
| 419 |  |  | _GLOBAL_SUM_RL( num_hmean_mm , mythid ) | 
| 420 |  |  | write(standardmessageunit,'(A,D22.15)') | 
| 421 |  |  | & 'CMM DOT tmp ref: offset ', offset | 
| 422 |  |  | write(standardmessageunit,'(A,D22.15)') | 
| 423 |  |  | & 'CMM DOT tmp ref: offset N ', offset_sum | 
| 424 |  |  | write(standardmessageunit,'(A,D22.15)') | 
| 425 |  |  | & 'CMM DOT tmp ref: cost integral ', sumc | 
| 426 |  |  | write(standardmessageunit,'(A,D22.15)') | 
| 427 |  |  | & 'CMM DOT tmp ref: cost N ', num_hmean_mm | 
| 428 |  |  | CMM  now to | 
| 429 |  |  | factor =  0.01 | 
| 430 |  |  | spval  = -9990. | 
| 431 |  |  | write(standardmessageunit,'(A)') | 
| 432 |  |  | & 'CMM ssh_cost - reading DOT in order: EGM08 AVISO GRACE MAXI' | 
| 433 |  |  | do IDOT = 1,4 | 
| 434 |  |  | CMM reinitialize due to recomps | 
| 435 |  |  | do bj = jtlo,jthi | 
| 436 |  |  | do bi = itlo,ithi | 
| 437 |  |  | do j = jmin,jmax | 
| 438 |  |  | do i = imin,imax | 
| 439 |  |  | tpmean_loc(i,j,bi,bj) = 0. _d 0 | 
| 440 |  |  | enddo | 
| 441 |  |  | enddo | 
| 442 |  |  | enddo | 
| 443 |  |  | enddo | 
| 444 |  |  | sumc = 0. | 
| 445 |  |  | num_hmean_mm = 0. | 
| 446 |  |  | if     (IDOT.eq.1) then | 
| 447 |  |  | call mdsreadfield( 'DOT_egm08.bin', cost_iprec,cost_yftype | 
| 448 |  |  | &                  ,1,tpmean_loc, 1, mythid ) | 
| 449 |  |  | elseif (IDOT.eq.2) then | 
| 450 |  |  | call mdsreadfield( 'DOT_aviso.bin', cost_iprec,cost_yftype, | 
| 451 |  |  | &                   1,tpmean_loc, 1, mythid ) | 
| 452 |  |  | elseif (IDOT.eq.3) then | 
| 453 |  |  | call mdsreadfield( 'DOT_grace2.bin', cost_iprec,cost_yftype, | 
| 454 |  |  | &                   1,tpmean_loc, 1, mythid ) | 
| 455 |  |  | elseif (IDOT.eq.4) then | 
| 456 |  |  | call mdsreadfield( 'DOT_maxi.bin', cost_iprec,cost_yftype, | 
| 457 |  |  | &                   1,tpmean_loc, 1, mythid ) | 
| 458 |  |  | endif | 
| 459 |  |  | CMM  NOW MAKE MASK | 
| 460 |  |  | do bj = jtlo,jthi | 
| 461 |  |  | do bi = itlo,ithi | 
| 462 |  |  | k = 1 | 
| 463 |  |  | do j = jmin,jmax | 
| 464 |  |  | do i = imin,imax | 
| 465 |  |  | mean_psMssh_all_MSK(i,j,bi,bj) =  tpmeanmask(i,j,bi,bj) | 
| 466 |  |  | if (tpmean_loc(i,j,bi,bj) .lt. spval) then | 
| 467 |  |  | mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0 | 
| 468 |  |  | endif | 
| 469 |  |  | if (tpmean_loc(i,j,bi,bj) .eq. 0. _d 0 ) then | 
| 470 |  |  | mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0 | 
| 471 |  |  | endif | 
| 472 |  |  | tpmean_loc(i,j,bi,bj) = tpmean_loc(i,j,bi,bj)* | 
| 473 |  |  | &                            mean_psMssh_all_MSK(i,j,bi,bj)* | 
| 474 |  |  | &                             factor | 
| 475 |  |  | enddo | 
| 476 |  |  | enddo | 
| 477 |  |  | enddo | 
| 478 |  |  | enddo | 
| 479 |  |  | CMM NOW GET OFFSET | 
| 480 |  |  | c--   Compute and remove offset for current tile and sum over all | 
| 481 |  |  | c--   tiles of this instance. | 
| 482 |  |  | offset     = 0. _d 0 | 
| 483 |  |  | offset_sum = 0. _d 0 | 
| 484 |  |  | c--   Sum over this thread tiles. | 
| 485 |  |  | do bj = jtlo,jthi | 
| 486 |  |  | do bi = itlo,ithi | 
| 487 |  |  | do j = 1,sny | 
| 488 |  |  | do i = 1,snx | 
| 489 |  |  | offset     = offset + | 
| 490 |  |  | &                mean_psMssh_all_MSK(i,j,bi,bj)*cosphi(i,j,bi,bj)* | 
| 491 |  |  | &                     (tpmean_loc(i,j,bi,bj) - psmean(i,j,bi,bj)) | 
| 492 |  |  | offset_sum = offset_sum + | 
| 493 |  |  | &                mean_psMssh_all_MSK(i,j,bi,bj)*cosphi(i,j,bi,bj) | 
| 494 |  |  | enddo | 
| 495 |  |  | enddo | 
| 496 |  |  | enddo | 
| 497 |  |  | enddo | 
| 498 |  |  |  | 
| 499 |  |  | c--   Do a global summation. | 
| 500 |  |  | _GLOBAL_SUM_RL( offset     , mythid ) | 
| 501 |  |  | _GLOBAL_SUM_RL( offset_sum , mythid ) | 
| 502 |  |  |  | 
| 503 |  |  | if (offset_sum .eq. 0.0) then | 
| 504 |  |  | _BEGIN_MASTER( mythid ) | 
| 505 |  |  | write(msgbuf,'(a)') ' cost_ssh: offset_sum = zero!' | 
| 506 |  |  | call print_message( msgbuf, standardmessageunit, | 
| 507 |  |  | &                          SQUEEZE_RIGHT , mythid) | 
| 508 |  |  | _END_MASTER( mythid ) | 
| 509 |  |  | stop   '  ... stopped in cost_ssh.' | 
| 510 |  |  | else | 
| 511 |  |  | _BEGIN_MASTER( mythid ) | 
| 512 |  |  | write(msgbuf,'(a,d22.15)') | 
| 513 |  |  | &          ' cost_ssh: offset = ',offset | 
| 514 |  |  | call print_message( msgbuf, standardmessageunit, | 
| 515 |  |  | &                          SQUEEZE_RIGHT , mythid) | 
| 516 |  |  | write(msgbuf,'(a,d22.15)') | 
| 517 |  |  | &          ' cost_ssh: offset_sum = ',offset_sum | 
| 518 |  |  | call print_message( msgbuf, standardmessageunit, | 
| 519 |  |  | &                          SQUEEZE_RIGHT , mythid) | 
| 520 |  |  |  | 
| 521 |  |  | _END_MASTER( mythid ) | 
| 522 |  |  | endif | 
| 523 |  |  | offset = offset / offset_sum | 
| 524 |  |  | do bj = jtlo,jthi | 
| 525 |  |  | do bi = itlo,ithi | 
| 526 |  |  | do j = jmin,jmax | 
| 527 |  |  | do i = imin,imax | 
| 528 |  |  | diagnosfld(i,j,bi,bj)=(psmean(i,j,bi,bj)-tpmean_loc(i,j,bi,bj) | 
| 529 |  |  | &              + offset)*mean_psMssh_all_MSK(i,j,bi,bj) | 
| 530 |  |  | diff = psmean(i,j,bi,bj) - tpmean_loc(i,j,bi,bj) + offset | 
| 531 |  |  | sumc = sumc + diff*diff | 
| 532 |  |  | &           *wp(i,j,bi,bj)*mean_psMssh_all_MSK(i,j,bi,bj) | 
| 533 |  |  | if ( wp(i,j,bi,bj)*mean_psMssh_all_MSK(i,j,bi,bj) .ne. 0. ) | 
| 534 |  |  | &             num_hmean_mm = num_hmean_mm + 1. _d 0 | 
| 535 |  |  | CMM              diagnosfld(i,j,bi,bj) = diff*diff | 
| 536 |  |  | CMM     &           *wp(i,j,bi,bj)*mean_psMssh_all_MSK(i,j,bi,bj) | 
| 537 |  |  | enddo | 
| 538 |  |  | enddo | 
| 539 |  |  | enddo | 
| 540 |  |  | enddo | 
| 541 |  |  |  | 
| 542 |  |  | if     (IDOT.eq.1) then | 
| 543 |  |  | CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean_mm_A', ' ', diagnosfld, | 
| 544 |  |  | &                           optimcycle, mythid ) | 
| 545 |  |  | elseif (IDOT.eq.2) then | 
| 546 |  |  | CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean_mm_B', ' ', diagnosfld, | 
| 547 |  |  | &                           optimcycle, mythid ) | 
| 548 |  |  | elseif (IDOT.eq.3) then | 
| 549 |  |  | CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean_mm_C', ' ', diagnosfld, | 
| 550 |  |  | &                           optimcycle, mythid ) | 
| 551 |  |  | elseif (IDOT.eq.4) then | 
| 552 |  |  | CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean_mm_D', ' ', diagnosfld, | 
| 553 |  |  | &                           optimcycle, mythid ) | 
| 554 |  |  | endif | 
| 555 |  |  |  | 
| 556 |  |  | CMM ADD TO COST FUNCTION | 
| 557 |  |  | objf_hmean = objf_hmean + sumc | 
| 558 |  |  | num_hmean = num_hmean + num_hmean_mm | 
| 559 |  |  |  | 
| 560 |  |  | c--   Do a global summation. | 
| 561 |  |  | _GLOBAL_SUM_RL( sumc     , mythid ) | 
| 562 |  |  | _GLOBAL_SUM_RL( num_hmean_mm , mythid ) | 
| 563 |  |  |  | 
| 564 |  |  | write(standardmessageunit,'(A,I6)') | 
| 565 |  |  | & 'CMM DOT sta ref number: ', IDOT | 
| 566 |  |  | write(standardmessageunit,'(A,D22.15)') | 
| 567 |  |  | & 'CMM DOT sta ref: offset ', offset | 
| 568 |  |  | write(standardmessageunit,'(A,D22.15)') | 
| 569 |  |  | & 'CMM DOT sta ref: offset N ', offset_sum | 
| 570 |  |  | write(standardmessageunit,'(A,D22.15)') | 
| 571 |  |  | & 'CMM DOT sta ref: cost integral ', sumc | 
| 572 |  |  | write(standardmessageunit,'(A,D22.15)') | 
| 573 |  |  | & 'CMM DOT sta ref: cost N ', num_hmean_mm | 
| 574 |  |  |  | 
| 575 |  |  | enddo | 
| 576 |  |  | CMM end of loop over DOT products | 
| 577 |  |  |  | 
| 578 |  |  |  | 
| 579 |  |  | #endif /* ALLOW_SSH_TOT */ | 
| 580 |  |  | #endif /* ALLOW_SSH_MEAN_COST_CONTRIBUTION */ | 
| 581 |  |  |  | 
| 582 |  |  | c--   ========== | 
| 583 |  |  | c--   Anomalies. | 
| 584 |  |  | c--   ========== | 
| 585 |  |  |  | 
| 586 |  |  | c--   Loop over records for the second time. | 
| 587 |  |  | do irec = 1, ndaysrec | 
| 588 |  |  |  | 
| 589 |  |  | call active_read_xy( fname, psbar, irec, doglobalread, | 
| 590 |  |  | &                       ladinit, optimcycle, mythid, | 
| 591 |  |  | &                       xx_psbar_mean_dummy ) | 
| 592 |  |  |  | 
| 593 |  |  | #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION | 
| 594 |  |  | call cost_readtopex( irec, mythid ) | 
| 595 |  |  | #endif | 
| 596 |  |  | #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION | 
| 597 |  |  | call cost_readers( irec, mythid ) | 
| 598 |  |  | #endif | 
| 599 |  |  | #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION | 
| 600 |  |  | call cost_readgfo( irec, mythid ) | 
| 601 |  |  | #endif | 
| 602 |  |  |  | 
| 603 |  |  | do bj = jtlo,jthi | 
| 604 |  |  | do bi = itlo,ithi | 
| 605 |  |  |  | 
| 606 |  |  | #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION | 
| 607 |  |  | do j = jmin,jmax | 
| 608 |  |  | do i = imin,imax | 
| 609 |  |  | c--             The array psobs contains SSH anomalies. | 
| 610 |  |  | wwwtp2(i,j,bi,bj) = wwwtp1(i,j,bi,bj) | 
| 611 |  |  | &                             *wtp(i,j,bi,bj) | 
| 612 |  |  | &                             *tpmask(i,j,bi,bj) | 
| 613 |  |  | &                             *cosphi(i,j,bi,bj) | 
| 614 |  |  | #ifndef ALLOW_SSH_TOT | 
| 615 |  |  | CMM                junk = ( psbar(i,j,bi,bj) | 
| 616 |  |  | CMM     &                   - psmeantp(i,j,bi,bj) - tpobs(i,j,bi,bj) ) | 
| 617 |  |  | junk =  ( psbar(i,j,bi,bj) - psmeantp(i,j,bi,bj)  ) | 
| 618 |  |  | &                - ( tpobs(i,j,bi,bj) - obsmeantp(i,j,bi,bj) ) | 
| 619 |  |  | objf_tp(bi,bj) = objf_tp(bi,bj) | 
| 620 |  |  | &              + junk*junk*wwwtp2(i,j,bi,bj) | 
| 621 |  |  | if ( wwwtp2(i,j,bi,bj) .ne. 0. ) | 
| 622 |  |  | &               num_tp(bi,bj) = num_tp(bi,bj) + 1. _d 0 | 
| 623 |  |  | #else | 
| 624 |  |  | if (tpmeanmask(i,j,bi,bj)*tpmask(i,j,bi,bj) | 
| 625 |  |  | &                  *wp(i,j,bi,bj)*wwwtp2(i,j,bi,bj) .ne.0.) then | 
| 626 |  |  | junk       = ( psbar(i,j,bi,bj) - | 
| 627 |  |  | &                 (tpobs(i,j,bi,bj)+tpmean(i,j,bi,bj)-offset) ) | 
| 628 |  |  | objf_tp(bi,bj) = objf_tp(bi,bj) | 
| 629 |  |  | &              +junk*junk/( 1/wp(i,j,bi,bj)+1/wwwtp2(i,j,bi,bj) ) | 
| 630 |  |  | num_tp(bi,bj) = num_tp(bi,bj) + 1. _d 0 | 
| 631 |  |  | endif | 
| 632 |  |  | #endif /* ALLOW_SSH_TOT */ | 
| 633 |  |  | enddo | 
| 634 |  |  | enddo | 
| 635 |  |  | #endif | 
| 636 |  |  |  | 
| 637 |  |  | #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION | 
| 638 |  |  | do j = jmin,jmax | 
| 639 |  |  | do i = imin,imax | 
| 640 |  |  | c--             The array ersobs contains SSH anomalies. | 
| 641 |  |  | wwwers2(i,j,bi,bj) = wwwers1(i,j,bi,bj) | 
| 642 |  |  | &                              *wers(i,j,bi,bj) | 
| 643 |  |  | &                              *ersmask(i,j,bi,bj) | 
| 644 |  |  | &                              *cosphi(i,j,bi,bj) | 
| 645 |  |  | #ifndef ALLOW_SSH_TOT | 
| 646 |  |  | CMM                junk = ( psbar(i,j,bi,bj) | 
| 647 |  |  | CMM     &                   - psmeaners(i,j,bi,bj) - ersobs(i,j,bi,bj) ) | 
| 648 |  |  | junk =  ( psbar(i,j,bi,bj) - psmeaners(i,j,bi,bj)  ) | 
| 649 |  |  | &                - ( ersobs(i,j,bi,bj) - obsmeaners(i,j,bi,bj) ) | 
| 650 |  |  | objf_ers(bi,bj) = objf_ers(bi,bj) | 
| 651 |  |  | &              + junk*junk*wwwers2(i,j,bi,bj) | 
| 652 |  |  | if ( wwwers2(i,j,bi,bj) .ne. 0. ) | 
| 653 |  |  | &               num_ers(bi,bj) = num_ers(bi,bj) + 1. _d 0 | 
| 654 |  |  | #else | 
| 655 |  |  | if (tpmeanmask(i,j,bi,bj)*ersmask(i,j,bi,bj) | 
| 656 |  |  | &                  *wp(i,j,bi,bj)*wwwers2(i,j,bi,bj) .ne.0.) then | 
| 657 |  |  | junk       = ( psbar(i,j,bi,bj) - | 
| 658 |  |  | &                 (ersobs(i,j,bi,bj)+tpmean(i,j,bi,bj)-offset) ) | 
| 659 |  |  | objf_ers(bi,bj) = objf_ers(bi,bj) | 
| 660 |  |  | &              +junk*junk/( 1/wp(i,j,bi,bj)+1/wwwers2(i,j,bi,bj) ) | 
| 661 |  |  | num_ers(bi,bj) = num_ers(bi,bj) + 1. _d 0 | 
| 662 |  |  | endif | 
| 663 |  |  | #endif /* ALLOW_SSH_TOT */ | 
| 664 |  |  | enddo | 
| 665 |  |  | enddo | 
| 666 |  |  | #endif | 
| 667 |  |  |  | 
| 668 |  |  | #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION | 
| 669 |  |  | do j = jmin,jmax | 
| 670 |  |  | do i = imin,imax | 
| 671 |  |  | c--             The array gfoobs contains SSH anomalies. | 
| 672 |  |  | wwwgfo2(i,j,bi,bj) = wwwgfo1(i,j,bi,bj) | 
| 673 |  |  | &                        *wgfo(i,j,bi,bj) | 
| 674 |  |  | &                        *gfomask(i,j,bi,bj) | 
| 675 |  |  | &                        *cosphi(i,j,bi,bj) | 
| 676 |  |  | #ifndef ALLOW_SSH_TOT | 
| 677 |  |  | CMM                junk = ( psbar(i,j,bi,bj) | 
| 678 |  |  | CMM     &                   - psmeangfo(i,j,bi,bj) - gfoobs(i,j,bi,bj) ) | 
| 679 |  |  | junk =  ( psbar(i,j,bi,bj) - psmeangfo(i,j,bi,bj)  ) | 
| 680 |  |  | &                - ( gfoobs(i,j,bi,bj) - obsmeangfo(i,j,bi,bj) ) | 
| 681 |  |  | objf_gfo(bi,bj) = objf_gfo(bi,bj) | 
| 682 |  |  | &              + junk*junk*wwwgfo2(i,j,bi,bj) | 
| 683 |  |  | if ( wwwgfo2(i,j,bi,bj) .ne. 0. ) | 
| 684 |  |  | &               num_gfo(bi,bj) = num_gfo(bi,bj) + 1. _d 0 | 
| 685 |  |  | #else | 
| 686 |  |  | if (tpmeanmask(i,j,bi,bj)*gfomask(i,j,bi,bj) | 
| 687 |  |  | &                  *wp(i,j,bi,bj)*wwwgfo2(i,j,bi,bj) .ne.0.) then | 
| 688 |  |  | junk       = ( psbar(i,j,bi,bj) - | 
| 689 |  |  | &                 (gfoobs(i,j,bi,bj)+tpmean(i,j,bi,bj)-offset) ) | 
| 690 |  |  | objf_gfo(bi,bj) = objf_gfo(bi,bj) | 
| 691 |  |  | &              +junk*junk/( 1/wp(i,j,bi,bj)+1/wwwgfo2(i,j,bi,bj) ) | 
| 692 |  |  | num_gfo(bi,bj) = num_gfo(bi,bj) + 1. _d 0 | 
| 693 |  |  | endif | 
| 694 |  |  | #endif /* ALLOW_SSH_TOT */ | 
| 695 |  |  | enddo | 
| 696 |  |  | enddo | 
| 697 |  |  | #endif | 
| 698 |  |  |  | 
| 699 |  |  | enddo | 
| 700 |  |  | enddo | 
| 701 |  |  |  | 
| 702 |  |  | enddo | 
| 703 |  |  | c--   End of second loop over records. | 
| 704 |  |  |  | 
| 705 |  |  | do bj = jtlo,jthi | 
| 706 |  |  | do bi = itlo,ithi | 
| 707 |  |  | objf_h(bi,bj) = objf_h(bi,bj) + | 
| 708 |  |  | &        objf_tp(bi,bj) + objf_ers(bi,bj) + objf_gfo(bi,bj) | 
| 709 |  |  | num_h(bi,bj) = num_h(bi,bj) + | 
| 710 |  |  | &        num_tp(bi,bj) + num_ers(bi,bj) + num_gfo(bi,bj) | 
| 711 |  |  | enddo | 
| 712 |  |  | enddo | 
| 713 |  |  |  | 
| 714 |  |  |  | 
| 715 |  |  | #endif /* ifdef ALLOW_SSH_COST_CONTRIBUTION */ | 
| 716 |  |  |  | 
| 717 |  |  | end |