/[MITgcm]/MITgcm_contrib/SOSE/code_ad/cost_ssh_new.F
ViewVC logotype

Contents of /MITgcm_contrib/SOSE/code_ad/cost_ssh_new.F

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


Revision 1.1 - (show annotations) (download)
Fri Apr 23 19:55:11 2010 UTC (15 years, 3 months ago) by mmazloff
Branch: MAIN
CVS Tags: HEAD
original files

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

  ViewVC Help
Powered by ViewVC 1.1.22