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

Annotation 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 - (hide 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 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

  ViewVC Help
Powered by ViewVC 1.1.22