/[MITgcm]/MITgcm/pkg/ecco/cost_gencost_sshv4.F
ViewVC logotype

Annotation of /MITgcm/pkg/ecco/cost_gencost_sshv4.F

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

Revision 1.2 - (hide annotations) (download)
Mon Oct 25 20:42:31 2010 UTC (13 years, 7 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint62o, checkpoint62n, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62x
Changes since 1.1: +2 -2 lines 
- cost_gencost_sshv4.F and cost_sshv4.F: compute smoothed
  35 day-average fields once a week only -- for speed up.
- ecco_cost_driver.F: necessary includes.
1 gforget 1.2 C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_gencost_sshv4.F,v 1.1 2010/09/07 21:20:39 gforget Exp $
2 gforget 1.1 C $Name: $
3    
4     #include "COST_CPPOPTIONS.h"
5    
6    
7     subroutine cost_gencost_sshv4(
8     I myiter,
9     I mytime,
10     I mythid
11     & )
12    
13     c ==================================================================
14     c SUBROUTINE cost_gencost_sshv4
15     c ==================================================================
16     c
17     c o Evaluate cost function contributions of sea surface height.
18     c
19     c started: Gael Forget, Oct-2009
20     c
21     c working assumption for the time mean dynamic topography (MDT) constraint:
22     c the various SLA data sets (tp, ers, gfo) have been consistenty cross-referenced,
23     c as done in the RADS data sets. We do not need to know the reference dynamic
24     c topography (or SSH/Geoid). Simply it is assumed that the biases
25     c between instruments have been taken care of. This is only a matter
26     c for the MDT constraint, not for SLA constraints (see below).
27     c
28     cgf 1) there are a few hardcoded numbers that will eventually be put in common
29     cgf blocks/namelists
30     cgf 2) there are a several refinements that should be considered, such as
31     cgf modulating weights with respect to numbers of samples
32     c
33     c ==================================================================
34     c SUBROUTINE cost_gencost_sshv4
35     c ==================================================================
36    
37     implicit none
38    
39     c == global variables ==
40    
41     #include "EEPARAMS.h"
42     #include "SIZE.h"
43     #include "PARAMS.h"
44     #include "GRID.h"
45    
46     #include "ecco_cost.h"
47     #include "ctrl.h"
48     #include "ctrl_dummy.h"
49     #include "optim.h"
50     #include "DYNVARS.h"
51     #ifdef ALLOW_PROFILES
52     #include "profiles.h"
53     #endif
54    
55     c == routine arguments ==
56    
57     integer myiter
58     _RL mytime
59     integer mythid
60    
61     #ifdef ALLOW_SSH_COST_CONTRIBUTION
62     #ifdef ALLOW_SMOOTH
63     #ifdef ALLOW_GENCOST_CONTRIBUTION
64     #ifdef ALLOW_GENCOST_FREEFORM
65     c == local variables ==
66    
67     integer bi,bj
68     integer i,j,k
69     integer itlo,ithi
70     integer jtlo,jthi
71     integer jmin,jmax
72     integer imin,imax
73     integer irec,jrec,krec
74     integer ilps
75     integer gwunit
76    
77     logical doglobalread
78     logical ladinit
79    
80     c mapping to gencost
81     integer igen_mdt, igen_lsc
82     integer igen_tp, igen_ers, igen_gfo
83    
84     _RL offset,fac
85     _RL offset_sum
86     _RL psmean ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
87     _RL diagnosfld ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
88    
89     c for PART 1: re-reference MDT (tpmean) to the inferred SLA reference field
90     _RL mean_slaobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
91     _RL mean_slaobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
92    
93     c for PART 2: compute time mean differences over the model period
94     _RL mean_slaobs2(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
95     _RL mean_psMssh_all(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
96     _RL mean_psMssh_all_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
97     _RL mean_psMssh_all_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
98    
99     _RL mean_psMtpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
100     _RL mean_psMtpobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
101     _RL mean_psMtpobs_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
102    
103     _RL mean_psMersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
104     _RL mean_psMersobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
105     _RL mean_psMersobs_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
106    
107     _RL mean_psMgfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
108     _RL mean_psMgfoobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
109     _RL mean_psMgfoobs_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
110    
111     c for PART 4/5: compute smooth/raw anomalies
112     _RL anom_psMslaobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
113     _RL anom_slaobs (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
114     _RL anom_psMslaobs_NUM (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
115    
116     _RL anom_psMtpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
117     _RL anom_psMtpobs_NUM (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
118     _RL anom_tpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
119    
120     _RL anom_psMersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
121     _RL anom_psMersobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
122     _RL anom_ersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
123    
124     _RL anom_psMgfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
125     _RL anom_psMgfoobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
126     _RL anom_gfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
127    
128     integer tpmean_y0,tpmean_y1,year,day
129     integer num_var
130    
131     _RL junk,junkweight
132    
133     integer ndaysave
134     _RL ndaysaveRL
135    
136     character*(80) fname
137     character*(80) fname4test
138     character*(MAX_LEN_MBUF) msgbuf
139    
140     c == external functions ==
141    
142     integer ilnblnk
143     external ilnblnk
144    
145     c == end of interface ==
146    
147     jtlo = mybylo(mythid)
148     jthi = mybyhi(mythid)
149     itlo = mybxlo(mythid)
150     ithi = mybxhi(mythid)
151     jmin = 1
152     jmax = sny
153     imin = 1
154     imax = snx
155    
156     c-- detect the relevant gencost indices
157     igen_mdt=0
158     igen_tp =0
159     igen_ers=0
160     igen_gfo=0
161     igen_lsc=0
162     do k=1,NGENCOST
163     if (gencost_name(k).EQ.'sshv4-mdt') igen_mdt=k
164     if (gencost_name(k).EQ.'sshv4-tp') igen_tp=k
165     if (gencost_name(k).EQ.'sshv4-ers') igen_ers=k
166     if (gencost_name(k).EQ.'sshv4-gfo') igen_gfo=k
167     if (gencost_name(k).EQ.'sshv4-lsc') igen_lsc=k
168     enddo
169    
170     c-- First, read tiled data.
171     doglobalread = .false.
172     ladinit = .false.
173    
174     write(fname(1:80),'(80a)') ' '
175     ilps=ilnblnk( psbarfile )
176     write(fname(1:80),'(2a,i10.10)')
177     & psbarfile(1:ilps),'.',optimcycle
178    
179    
180     cgf =======================================================
181     cgf PART 1:
182     cgf x Get the MDT (tpmean) ... compute the sample mean
183     cgf (mean_slaobs) of the SLA data (i.e. RADS for tp, ers, and gfo
184     cgf together) over the time interval of the MDT ... subtract
185     cgf mean_slaobs from tpmean.
186     cgf x At this point, tpmean is the inferred SLA reference field.
187     cgf x From there on, tpmean+sla will be directly comparable to
188     cgf the model SSH (psbar).
189     cgf =======================================================
190    
191     c-- Read mean field and mask
192     call cost_ReadTopexMean( mythid )
193    
194     c-- Compute mean_slaobs: sample mean SLA over the time period of tpmean.
195    
196     c pavlis and ecco/rio
197     tpmean_y0=1993
198     tpmean_y1=2004
199     c maximenko
200     c tpmean_y0=1992
201     c tpmean_y1=2002
202     c rio
203     c tpmean_y0=1993
204     c tpmean_y1=1999
205    
206     do bj = jtlo,jthi
207     do bi = itlo,ithi
208     do j = jmin,jmax
209     do i = imin,imax
210     mean_slaobs(i,j,bi,bj) = 0. _d 0
211     mean_slaobs_NUM(i,j,bi,bj) = 0. _d 0
212     enddo
213     enddo
214     enddo
215     enddo
216    
217     do year=tpmean_y0,tpmean_y1
218     do day=1,366
219     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
220     call cost_sla_read_yd( topexfile,
221     & tpobs, tpmask,
222     & year, day, mythid )
223     #endif
224     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
225     call cost_sla_read_yd( ersfile,
226     & ersobs, ersmask,
227     & year, day, mythid )
228     #endif
229     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
230     call cost_sla_read_yd( gfofile,
231     & gfoobs, gfomask,
232     & year, day, mythid )
233     #endif
234     do bj = jtlo,jthi
235     do bi = itlo,ithi
236     do j = jmin,jmax
237     do i = imin,imax
238     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
239     if ( tpmask(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)*
240     & gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then
241     mean_slaobs(i,j,bi,bj)= mean_slaobs(i,j,bi,bj)+
242     & tpobs(i,j,bi,bj)
243     mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
244     endif
245     #endif
246     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
247     if ( ersmask(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)*
248     & gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then
249     mean_slaobs(i,j,bi,bj)= mean_slaobs(i,j,bi,bj)+
250     & ersobs(i,j,bi,bj)
251     mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
252     endif
253     #endif
254     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
255     if ( gfomask(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)*
256     & gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then
257     mean_slaobs(i,j,bi,bj)= mean_slaobs(i,j,bi,bj)+
258     & gfoobs(i,j,bi,bj)
259     mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
260     endif
261     #endif
262     enddo
263     enddo
264     enddo
265     enddo
266    
267     enddo !do day=1,366
268     enddo !do year=tpmean_y0,tpmean_y1
269    
270     do bj = jtlo,jthi
271     do bi = itlo,ithi
272     do j = jmin,jmax
273     do i = imin,imax
274     if ( ( mean_slaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
275     & ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
276     & ( abs(YC(i,j,bi,bj)) .LE. 66. ).AND.
277     & ( tpmeanmask(i,j,bi,bj) .NE. 0. ) ) then
278     mean_slaobs(i,j,bi,bj) = mean_slaobs(i,j,bi,bj) /
279     & mean_slaobs_NUM(i,j,bi,bj)
280     else
281     mean_slaobs(i,j,bi,bj) = 0. _d 0
282     endif
283     enddo
284     enddo
285     enddo
286     enddo
287    
288    
289     c-- smooth mean_slaobs:
290    
291     write(fname4test(1:80),'(1a)') 'sla2mdt_raw'
292     call mdswritefield(fname4test,32,.false.,'RL',
293     & 1,mean_slaobs,1,1,mythid)
294    
295     call smooth_hetero2d(mean_slaobs,maskc,
296     & gencost_scalefile(igen_mdt),300,mythid)
297    
298     write(fname4test(1:80),'(1a)') 'sla2mdt_smooth'
299     call mdswritefield(fname4test,32,.false.,'RL',
300     & 1,mean_slaobs,1,1,mythid)
301    
302     c-- re-reference tpmean:
303     do bj = jtlo,jthi
304     do bi = itlo,ithi
305     do j = jmin,jmax
306     do i = imin,imax
307     if ( ( tpmeanmask(i,j,bi,bj) .NE. 0. ).AND.
308     & ( maskc(i,j,1,bi,bj) .NE. 0. ) ) then
309     tpmean(i,j,bi,bj) = tpmean(i,j,bi,bj)
310     & -mean_slaobs(i,j,bi,bj)
311     endif
312     enddo
313     enddo
314     enddo
315     enddo
316    
317    
318     cgf =======================================================
319     cgf PART 2: compute sample means of psbar-slaobs over the
320     cgf period that is covered by the model (i.e. psbar).
321     cgf x for all SLA data sets together: mean_psMssh_all, mean_psMssh_all_MSK,
322     cgf and offset will be used in PART 3 (MDT cost term).
323     cgf x for each SLA data individually. mean_psMtpobs, mean_psMtpobs_MS, etc.
324     cgf will be used in PART 4&5 (SLA cost terms).
325     cgf =======================================================
326    
327     do bj = jtlo,jthi
328     do bi = itlo,ithi
329     do j = jmin,jmax
330     do i = imin,imax
331    
332     psmean(i,j,bi,bj) = 0. _d 0
333     mean_psMtpobs(i,j,bi,bj) = 0. _d 0
334     mean_psMersobs(i,j,bi,bj) = 0. _d 0
335     mean_psMgfoobs(i,j,bi,bj) = 0. _d 0
336     mean_psMssh_all(i,j,bi,bj) = 0. _d 0
337     mean_slaobs2(i,j,bi,bj) = 0. _d 0
338    
339     mean_psMtpobs_NUM(i,j,bi,bj) = 0. _d 0
340     mean_psMersobs_NUM(i,j,bi,bj) = 0. _d 0
341     mean_psMgfoobs_NUM(i,j,bi,bj) = 0. _d 0
342     mean_psMssh_all_NUM(i,j,bi,bj) = 0. _d 0
343    
344     mean_psMtpobs_MSK(i,j,bi,bj) = 0. _d 0
345     mean_psMersobs_MSK(i,j,bi,bj) = 0. _d 0
346     mean_psMgfoobs_MSK(i,j,bi,bj) = 0. _d 0
347    
348     enddo
349     enddo
350     enddo
351     enddo
352     offset = 0. _d 0
353     offset_sum = 0. _d 0
354    
355    
356     do irec = 1, ndaysrec
357    
358     call active_read_xy( fname, psbar, irec, doglobalread,
359     & ladinit, optimcycle, mythid,
360     & xx_psbar_mean_dummy )
361    
362     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
363     call cost_sla_read( topexfile, topexstartdate, topexperiod,
364     & topexintercept, topexslope,
365     & tpobs, tpmask,
366     & irec, mythid )
367     #endif
368     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
369     call cost_sla_read( ersfile, ersstartdate, ersperiod,
370     & ersintercept, ersslope,
371     & ersobs, ersmask,
372     & irec, mythid )
373     #endif
374     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
375     call cost_sla_read( gfofile, gfostartdate, gfoperiod,
376     & gfointercept, gfoslope,
377     & gfoobs, gfomask,
378     & irec, mythid )
379     #endif
380    
381     do bj = jtlo,jthi
382     do bi = itlo,ithi
383     do j = jmin,jmax
384     do i = imin,imax
385     psmean(i,j,bi,bj) = psmean(i,j,bi,bj) +
386     & psbar(i,j,bi,bj) / float(ndaysrec)
387     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
388     if ( tpmask(i,j,bi,bj)*
389     & gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then
390     mean_slaobs2(i,j,bi,bj)=
391     & mean_slaobs2(i,j,bi,bj)+tpobs(i,j,bi,bj)
392     mean_psMtpobs(i,j,bi,bj) =
393     & mean_psMtpobs(i,j,bi,bj) +
394     & psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)
395     mean_psMtpobs_NUM(i,j,bi,bj) =
396     & mean_psMtpobs_NUM(i,j,bi,bj) + 1. _d 0
397     endif
398     #endif
399     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
400     if ( ersmask(i,j,bi,bj)*
401     & gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then
402     mean_slaobs2(i,j,bi,bj)=
403     & mean_slaobs2(i,j,bi,bj)+ersobs(i,j,bi,bj)
404     mean_psMersobs(i,j,bi,bj) =
405     & mean_psMersobs(i,j,bi,bj) +
406     & psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)
407     mean_psMersobs_NUM(i,j,bi,bj) =
408     & mean_psMersobs_NUM(i,j,bi,bj) + 1. _d 0
409     endif
410     #endif
411     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
412     if ( gfomask(i,j,bi,bj)*
413     & gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then
414     mean_slaobs2(i,j,bi,bj)=
415     & mean_slaobs2(i,j,bi,bj)+gfoobs(i,j,bi,bj)
416     mean_psMgfoobs(i,j,bi,bj) =
417     & mean_psMgfoobs(i,j,bi,bj) +
418     & psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)
419     mean_psMgfoobs_NUM(i,j,bi,bj) =
420     & mean_psMgfoobs_NUM(i,j,bi,bj) + 1. _d 0
421     endif
422     #endif
423     enddo
424     enddo
425     enddo
426     enddo
427    
428     c-- END loop over records for the first time.
429     enddo
430    
431     do bj = jtlo,jthi
432     do bi = itlo,ithi
433     do j = jmin,jmax
434     do i = imin,imax
435     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
436     if ( mean_psMtpobs_NUM(i,j,bi,bj) .NE. 0. ) then
437     mean_psMssh_all(i,j,bi,bj) =
438     & mean_psMssh_all(i,j,bi,bj) +
439     & mean_psMtpobs(i,j,bi,bj)
440     mean_psMssh_all_NUM(i,j,bi,bj) =
441     & mean_psMssh_all_NUM(i,j,bi,bj) +
442     & mean_psMtpobs_NUM(i,j,bi,bj)
443     mean_psMtpobs(i,j,bi,bj) =
444     & mean_psMtpobs(i,j,bi,bj) /
445     & mean_psMtpobs_NUM(i,j,bi,bj)
446     mean_psMtpobs_MSK(i,j,bi,bj) = 1. _d 0
447     endif
448     #endif
449     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
450     if ( mean_psMersobs_NUM(i,j,bi,bj) .NE. 0. ) then
451     mean_psMssh_all(i,j,bi,bj) =
452     & mean_psMssh_all(i,j,bi,bj) +
453     & mean_psMersobs(i,j,bi,bj)
454     mean_psMssh_all_NUM(i,j,bi,bj) =
455     & mean_psMssh_all_NUM(i,j,bi,bj) +
456     & mean_psMersobs_NUM(i,j,bi,bj)
457     mean_psMersobs(i,j,bi,bj) =
458     & mean_psMersobs(i,j,bi,bj) /
459     & mean_psMersobs_NUM(i,j,bi,bj)
460     mean_psMersobs_MSK(i,j,bi,bj) = 1. _d 0
461     endif
462     #endif
463     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
464     if ( mean_psMgfoobs_NUM(i,j,bi,bj) .NE. 0. ) then
465     mean_psMssh_all(i,j,bi,bj) =
466     & mean_psMssh_all(i,j,bi,bj) +
467     & mean_psMgfoobs(i,j,bi,bj)
468     mean_psMssh_all_NUM(i,j,bi,bj) =
469     & mean_psMssh_all_NUM(i,j,bi,bj) +
470     & mean_psMgfoobs_NUM(i,j,bi,bj)
471     mean_psMgfoobs(i,j,bi,bj) =
472     & mean_psMgfoobs(i,j,bi,bj) /
473     & mean_psMgfoobs_NUM(i,j,bi,bj)
474     mean_psMgfoobs_MSK(i,j,bi,bj) = 1. _d 0
475     endif
476     #endif
477     if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .NE. 0. ).AND.
478     & ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
479     & ( abs(YC(i,j,bi,bj)) .LE. 66. ).AND.
480     & ( tpmeanmask(i,j,bi,bj) .NE. 0. ) ) then
481     mean_slaobs2(i,j,bi,bj) =
482     & mean_slaobs2(i,j,bi,bj) /
483     & mean_psMssh_all_NUM(i,j,bi,bj)
484     mean_psMssh_all(i,j,bi,bj) =
485     & mean_psMssh_all(i,j,bi,bj) /
486     & mean_psMssh_all_NUM(i,j,bi,bj)-tpmean(i,j,bi,bj)
487     mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
488     offset=offset+mean_psMssh_all(i,j,bi,bj)*
489     & mean_psMssh_all_NUM(i,j,bi,bj)
490     offset_sum=offset_sum+mean_psMssh_all_NUM(i,j,bi,bj)
491     else
492     mean_slaobs2(i,j,bi,bj) = 0.d0
493     mean_psMssh_all(i,j,bi,bj) = 0. _d 0
494     mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0
495     endif
496     enddo
497     enddo
498     enddo
499     enddo
500    
501     c-- Do a global summation.
502     _GLOBAL_SUM_RL( offset , mythid )
503     _GLOBAL_SUM_RL( offset_sum , mythid )
504    
505     if (offset_sum .eq. 0.0) then
506     _BEGIN_MASTER( mythid )
507     write(msgbuf,'(a)') ' cost_ssh: offset_sum = zero!'
508     call print_message( msgbuf, standardmessageunit,
509     & SQUEEZE_RIGHT , mythid)
510     _END_MASTER( mythid )
511     c stop ' ... stopped in cost_ssh.'
512     else
513     _BEGIN_MASTER( mythid )
514     write(msgbuf,'(a,d22.15)')
515     & ' cost_ssh: offset_sum = ',offset_sum
516     call print_message( msgbuf, standardmessageunit,
517     & SQUEEZE_RIGHT , mythid)
518     _END_MASTER( mythid )
519     endif
520    
521     c-- Compute (average) offset
522     offset = offset / offset_sum
523    
524     c-- subtract offset from mean_psMssh_all and psmean
525     do bj = jtlo,jthi
526     do bi = itlo,ithi
527     do j = jmin,jmax
528     do i = imin,imax
529    
530     if ( (mean_psMssh_all_MSK(i,j,bi,bj) .NE. 0.) .AND.
531     & ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
532     & ( abs(YC(i,j,bi,bj)) .LE. 66. ) ) then
533     c use the re-referencing approach
534     mean_psMssh_all(i,j,bi,bj) =
535     & mean_psMssh_all(i,j,bi,bj) - offset
536     elseif ( ( tpmeanmask(i,j,bi,bj) .NE. 0. ) .AND.
537     & ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
538     & ( abs(YC(i,j,bi,bj)) .GT. 66. ) ) then
539     c use the simpler approach
540     mean_psMssh_all(i,j,bi,bj) =
541     & psmean(i,j,bi,bj) -tpmean(i,j,bi,bj) - offset
542     else
543     mean_psMssh_all(i,j,bi,bj) = 0. _d 0
544     endif
545    
546     if ( maskc(i,j,1,bi,bj) .NE. 0. )
547     & psmean(i,j,bi,bj)=psmean(i,j,bi,bj)-offset
548     enddo
549     enddo
550     enddo
551     enddo
552    
553     c-- smooth mean_psMssh_all
554     write(fname4test(1:80),'(1a)') 'mdtdiff_raw'
555     call mdswritefield(fname4test,32,.false.,'RL',
556     & 1,mean_psMssh_all,1,1,mythid)
557    
558     call smooth_hetero2d(mean_psMssh_all,maskc,
559     & gencost_scalefile(igen_mdt),300,mythid)
560    
561     write(fname4test(1:80),'(1a)') 'mdtdiff_smooth'
562     call mdswritefield(fname4test,32,.false.,'RL',
563     & 1,mean_psMssh_all,1,1,mythid)
564    
565     write(fname4test(1:80),'(1a)') 'sla2model_raw'
566     call mdswritefield(fname4test,32,.false.,'RL',
567     & 1,mean_slaobs2,1,1,mythid)
568    
569     call smooth_hetero2d(mean_slaobs2,maskc,
570     & gencost_scalefile(igen_mdt),300,mythid)
571    
572     write(fname4test(1:80),'(1a)') 'sla2model_smooth'
573     call mdswritefield(fname4test,32,.false.,'RL',
574     & 1,mean_slaobs2,1,1,mythid)
575    
576     cgf at this point:
577     cgf 1) mean_psMssh_all is the sample mean <psbar-slaobs-tpmean-offset>,
578     cgf to which a smoothing filter has been applied.
579     cgf 2) mean_psMtpobs is the (unsmoothed) sample mean <psbar-tpobs>.
580     cgf And similarly for ers and gfo, each treated separately.
581    
582     #ifdef ALLOW_PROFILES
583     do bj = jtlo,jthi
584     do bi = itlo,ithi
585     do j = 1,sny
586     do i = 1,snx
587     prof_etan_mean(i,j,bi,bj)=psmean(i,j,bi,bj)
588     enddo
589     enddo
590     enddo
591     enddo
592     _EXCH_XY_RL( prof_etan_mean, mythid )
593     #endif
594    
595    
596     cgf =======================================================
597     cgf PART 3: compute MDT cost term
598     cgf =======================================================
599    
600    
601     #ifdef ALLOW_SSH_MEAN_COST_CONTRIBUTION
602    
603     do bj = jtlo,jthi
604     do bi = itlo,ithi
605     do j = jmin,jmax
606     do i = imin,imax
607     junk = mean_psMssh_all(i,j,bi,bj)
608     junkweight = gencost_weight(i,j,bi,bj,igen_mdt)
609     & *tpmeanmask(i,j,bi,bj)
610     objf_gencost(igen_mdt,bi,bj) = objf_gencost(igen_mdt,bi,bj)
611     & + junk*junk*junkweight
612     if ( junkweight .ne. 0. ) num_gencost(igen_mdt,bi,bj) =
613     & num_gencost(igen_mdt,bi,bj) + 1. _d 0
614     diagnosfld(i,j,bi,bj) = junk*junk*junkweight
615     enddo
616     enddo
617     enddo
618     enddo
619    
620     CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean', ' ', diagnosfld,
621     & optimcycle, mythid )
622    
623     #endif /* ALLOW_SSH_MEAN_COST_CONTRIBUTION */
624    
625    
626    
627     cgf =======================================================
628     cgf PART 4: compute smooth SLA cost term
629     cgf =======================================================
630    
631    
632     ndaysave=35
633     ndaysaveRL=ndaysave
634    
635 gforget 1.2 do irec = 1, ndaysrec-ndaysave+1, 7
636 gforget 1.1
637     do bj = jtlo,jthi
638     do bi = itlo,ithi
639     do j = jmin,jmax
640     do i = imin,imax
641     anom_psMslaobs(i,j,bi,bj) = 0. _d 0
642     anom_slaobs(i,j,bi,bj) = 0. _d 0
643     anom_psMtpobs(i,j,bi,bj) = 0. _d 0
644     anom_psMersobs(i,j,bi,bj) = 0. _d 0
645     anom_psMgfoobs(i,j,bi,bj) = 0. _d 0
646     anom_psMslaobs_NUM(i,j,bi,bj) = 0. _d 0
647     anom_psMtpobs_NUM(i,j,bi,bj) = 0. _d 0
648     anom_psMersobs_NUM(i,j,bi,bj) = 0. _d 0
649     anom_psMgfoobs_NUM(i,j,bi,bj) = 0. _d 0
650     enddo
651     enddo
652     enddo
653     enddo
654    
655     c PART 4.1: compute running sample average over ndaysave
656     c ------------------------------------------------------
657    
658     do jrec=1,ndaysave
659    
660     krec=irec+jrec-1
661    
662     call active_read_xy( fname, psbar, krec, doglobalread,
663     & ladinit, optimcycle, mythid,
664     & xx_psbar_mean_dummy )
665    
666     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
667     call cost_sla_read( topexfile, topexstartdate, topexperiod,
668     & topexintercept, topexslope,
669     & tpobs, tpmask,
670     & krec, mythid )
671     #endif
672     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
673     call cost_sla_read( ersfile, ersstartdate, ersperiod,
674     & ersintercept, ersslope,
675     & ersobs, ersmask,
676     & krec, mythid )
677     #endif
678     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
679     call cost_sla_read( gfofile, gfostartdate, gfoperiod,
680     & gfointercept, gfoslope,
681     & gfoobs, gfomask,
682     & krec, mythid )
683     #endif
684    
685     do bj = jtlo,jthi
686     do bi = itlo,ithi
687     do j = jmin,jmax
688     do i = imin,imax
689     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
690     if ( tpmask(i,j,bi,bj)*mean_psMtpobs_MSK(i,j,bi,bj)
691     & .NE.0. ) then
692     anom_psMtpobs(i,j,bi,bj)= anom_psMtpobs(i,j,bi,bj)+
693     & psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)
694     & -mean_psMtpobs(i,j,bi,bj)
695     anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
696     & psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)
697     & -mean_psMtpobs(i,j,bi,bj)
698     anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
699     & tpobs(i,j,bi,bj)
700     anom_psMtpobs_NUM(i,j,bi,bj)=
701     & anom_psMtpobs_NUM(i,j,bi,bj)+1. _d 0
702     anom_psMslaobs_NUM(i,j,bi,bj)=
703     & anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
704     endif
705     #endif
706     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
707     if ( ersmask(i,j,bi,bj)*mean_psMersobs_MSK(i,j,bi,bj)
708     & .NE.0. ) then
709     anom_psMersobs(i,j,bi,bj)= anom_psMersobs(i,j,bi,bj)+
710     & psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)
711     & -mean_psMersobs(i,j,bi,bj)
712     anom_psMersobs_NUM(i,j,bi,bj)=
713     & anom_psMersobs_NUM(i,j,bi,bj)+1. _d 0
714     anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
715     & psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)
716     & -mean_psMersobs(i,j,bi,bj)
717     anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
718     & ersobs(i,j,bi,bj)
719     anom_psMslaobs_NUM(i,j,bi,bj)=
720     & anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
721     endif
722     #endif
723     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
724     if ( gfomask(i,j,bi,bj)*mean_psMgfoobs_MSK(i,j,bi,bj)
725     & .NE.0. ) then
726     anom_psMgfoobs(i,j,bi,bj)= anom_psMgfoobs(i,j,bi,bj)+
727     & psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)
728     & -mean_psMgfoobs(i,j,bi,bj)
729     anom_psMgfoobs_NUM(i,j,bi,bj)=
730     & anom_psMgfoobs_NUM(i,j,bi,bj)+1. _d 0
731     anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
732     & psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)
733     & -mean_psMgfoobs(i,j,bi,bj)
734     anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
735     & gfoobs(i,j,bi,bj)
736     anom_psMslaobs_NUM(i,j,bi,bj)=
737     & anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
738     endif
739     #endif
740     enddo
741     enddo
742     enddo
743     enddo
744    
745     enddo !do jrec=1,ndaysave
746    
747     do bj = jtlo,jthi
748     do bi = itlo,ithi
749     do j = jmin,jmax
750     do i = imin,imax
751     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
752     if ( anom_psMtpobs_NUM(i,j,bi,bj) .NE. 0. ) then
753     anom_psMtpobs(i,j,bi,bj) =
754     & anom_psMtpobs(i,j,bi,bj) /
755     & anom_psMtpobs_NUM(i,j,bi,bj)
756     endif
757     #endif
758     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
759     if ( anom_psMersobs_NUM(i,j,bi,bj) .NE. 0. ) then
760     anom_psMersobs(i,j,bi,bj) =
761     & anom_psMersobs(i,j,bi,bj) /
762     & anom_psMersobs_NUM(i,j,bi,bj)
763     endif
764     #endif
765     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
766     if ( anom_psMgfoobs_NUM(i,j,bi,bj) .NE. 0. ) then
767     anom_psMgfoobs(i,j,bi,bj) =
768     & anom_psMgfoobs(i,j,bi,bj) /
769     & anom_psMgfoobs_NUM(i,j,bi,bj)
770     endif
771     #endif
772     if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
773     & ( maskc(i,j,1,bi,bj) .NE. 0. ) ) then
774     anom_psMslaobs(i,j,bi,bj) =
775     & anom_psMslaobs(i,j,bi,bj) /
776     & anom_psMslaobs_NUM(i,j,bi,bj)
777     anom_slaobs(i,j,bi,bj) =
778     & anom_slaobs(i,j,bi,bj) /
779     & anom_psMslaobs_NUM(i,j,bi,bj)
780     else
781     anom_psMslaobs(i,j,bi,bj) = 0. _d 0
782     anom_slaobs(i,j,bi,bj) = 0. _d 0
783     endif
784     enddo
785     enddo
786     enddo
787     enddo
788    
789     c PART 4.2: smooth anom_psMslaobs in space
790     c ----------------------------------------
791    
792     #ifdef ALLOW_GENCOST_SSHV4_OUTPUT
793     write(fname4test(1:80),'(1a)') 'sladiff_raw'
794     call mdswritefield(fname4test,32,.false.,'RL',
795     & 1,anom_psMslaobs,irec,1,mythid)
796    
797     write(fname4test(1:80),'(1a)') 'slaobs_raw'
798     call mdswritefield(fname4test,32,.false.,'RL',
799     & 1,anom_slaobs,irec,1,mythid)
800     #endif
801    
802     call smooth_hetero2d(anom_psMslaobs,maskc,
803     & gencost_scalefile(igen_lsc),300,mythid)
804    
805     #ifdef ALLOW_GENCOST_SSHV4_OUTPUT
806     call smooth_hetero2d(anom_slaobs,maskc,
807     & gencost_scalefile(igen_lsc),300,mythid)
808    
809     write(fname4test(1:80),'(1a)') 'sladiff_smooth'
810     call mdswritefield(fname4test,32,.false.,'RL',
811     & 1,anom_psMslaobs,irec,1,mythid)
812    
813     write(fname4test(1:80),'(1a)') 'slaobs_smooth'
814     call mdswritefield(fname4test,32,.false.,'RL',
815     & 1,anom_slaobs,irec,1,mythid)
816     #endif
817    
818     c PART 4.3: compute cost function term
819     c ------------------------------------
820    
821     do bj = jtlo,jthi
822     do bi = itlo,ithi
823     do j = jmin,jmax
824     do i = imin,imax
825     # if (defined (ALLOW_SSH_GFOANOM_COST_CONTRIBUTION) || \
826     defined (ALLOW_SSH_TPANOM_COST_CONTRIBUTION) || \
827     defined (ALLOW_SSH_ERSANOM_COST_CONTRIBUTION))
828     junk = anom_psMslaobs(i,j,bi,bj)
829     objf_gencost(igen_lsc,bi,bj) = objf_gencost(igen_lsc,bi,bj)
830     & + junk*junk*gencost_weight(i,j,bi,bj,igen_lsc)
831     & *maskc(i,j,1,bi,bj)/ndaysaveRL
832     if ( (gencost_weight(i,j,bi,bj,igen_lsc).GT.0.).AND.
833     & (anom_psMslaobs_NUM(i,j,bi,bj).GT.0.).AND.
834     & (maskc(i,j,1,bi,bj) .ne. 0.) )
835     & num_gencost(igen_lsc,bi,bj) =
836     & num_gencost(igen_lsc,bi,bj) + 1. _d 0 /ndaysaveRL
837     #endif
838     enddo
839     enddo
840     enddo
841     enddo
842    
843     enddo
844    
845    
846     cgf =======================================================
847     cgf PART 5: compute raw SLA cost term
848     cgf =======================================================
849    
850    
851     do irec = 1, ndaysrec
852    
853     call active_read_xy( fname, psbar, irec, doglobalread,
854     & ladinit, optimcycle, mythid,
855     & xx_psbar_mean_dummy )
856    
857     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
858     call cost_readtopex( irec, mythid )
859     #endif
860     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
861     call cost_readers( irec, mythid )
862     #endif
863     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
864     call cost_readgfo( irec, mythid )
865     #endif
866    
867     do bj = jtlo,jthi
868     do bi = itlo,ithi
869     do j = jmin,jmax
870     do i = imin,imax
871     anom_psMtpobs(i,j,bi,bj) = 0. _d 0
872     anom_psMersobs(i,j,bi,bj) = 0. _d 0
873     anom_psMgfoobs(i,j,bi,bj) = 0. _d 0
874     anom_tpobs(i,j,bi,bj) = 0. _d 0
875     anom_ersobs(i,j,bi,bj) = 0. _d 0
876     anom_gfoobs(i,j,bi,bj) = 0. _d 0
877     enddo
878     enddo
879     enddo
880     enddo
881    
882     do bj = jtlo,jthi
883     do bi = itlo,ithi
884     do j = jmin,jmax
885     do i = imin,imax
886     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
887     if ( tpmask(i,j,bi,bj)*mean_psMtpobs_MSK(i,j,bi,bj).NE.0. )
888     & then
889     anom_psMtpobs(i,j,bi,bj)=
890     & psbar(i,j,bi,bj) - tpobs(i,j,bi,bj)
891     & - mean_psMtpobs(i,j,bi,bj)
892     anom_tpobs(i,j,bi,bj)=tpobs(i,j,bi,bj)
893     endif
894     #endif
895     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
896     if ( ersmask(i,j,bi,bj)*mean_psMersobs_MSK(i,j,bi,bj).NE.0. )
897     & then
898     anom_psMersobs(i,j,bi,bj)=
899     & psbar(i,j,bi,bj) - ersobs(i,j,bi,bj)
900     & - mean_psMersobs(i,j,bi,bj)
901     anom_ersobs(i,j,bi,bj)=ersobs(i,j,bi,bj)
902     endif
903     #endif
904     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
905     if ( gfomask(i,j,bi,bj)*mean_psMgfoobs_MSK(i,j,bi,bj).NE.0. )
906     & then
907     anom_psMgfoobs(i,j,bi,bj)=
908     & psbar(i,j,bi,bj) - gfoobs(i,j,bi,bj)
909     & - mean_psMgfoobs(i,j,bi,bj)
910     anom_gfoobs(i,j,bi,bj)=gfoobs(i,j,bi,bj)
911     endif
912     #endif
913     enddo
914     enddo
915     enddo
916     enddo
917    
918     #ifdef ALLOW_GENCOST_SSHV4_OUTPUT
919     write(fname4test(1:80),'(1a)') 'sladiff_tp_raw'
920     call mdswritefield(fname4test,32,.false.,'RL',
921     & 1,anom_psMtpobs,irec,1,mythid)
922     write(fname4test(1:80),'(1a)') 'sladiff_ers_raw'
923     call mdswritefield(fname4test,32,.false.,'RL',
924     & 1,anom_psMersobs,irec,1,mythid)
925     write(fname4test(1:80),'(1a)') 'sladiff_gfo_raw'
926     call mdswritefield(fname4test,32,.false.,'RL',
927     & 1,anom_psMgfoobs,irec,1,mythid)
928    
929     write(fname4test(1:80),'(1a)') 'slaobs_tp_raw'
930     call mdswritefield(fname4test,32,.false.,'RL',
931     & 1,anom_tpobs,irec,1,mythid)
932     write(fname4test(1:80),'(1a)') 'slaobs_ers_raw'
933     call mdswritefield(fname4test,32,.false.,'RL',
934     & 1,anom_ersobs,irec,1,mythid)
935     write(fname4test(1:80),'(1a)') 'slaobs_gfo_raw'
936     call mdswritefield(fname4test,32,.false.,'RL',
937     & 1,anom_gfoobs,irec,1,mythid)
938     #endif
939    
940     do bj = jtlo,jthi
941     do bi = itlo,ithi
942     do j = jmin,jmax
943     do i = imin,imax
944     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
945     c-- The array psobs contains SSH anomalies.
946     junkweight = mean_psMtpobs_MSK(i,j,bi,bj)
947     & *gencost_weight(i,j,bi,bj,igen_tp)
948     & *tpmask(i,j,bi,bj)
949     & *cosphi(i,j,bi,bj)
950     junk = anom_psMtpobs(i,j,bi,bj)
951     objf_gencost(igen_tp,bi,bj) =
952     & objf_gencost(igen_tp,bi,bj)+junk*junk*junkweight
953     if ( junkweight .ne. 0. )
954     & num_gencost(igen_tp,bi,bj) =
955     & num_gencost(igen_tp,bi,bj) + 1. _d 0
956     #endif
957     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
958     c-- The array ersobs contains SSH anomalies.
959     junkweight = mean_psMersobs_MSK(i,j,bi,bj)
960     & *gencost_weight(i,j,bi,bj,igen_ers)
961     & *ersmask(i,j,bi,bj)
962     & *cosphi(i,j,bi,bj)
963     junk = anom_psMersobs(i,j,bi,bj)
964     objf_gencost(igen_ers,bi,bj) =
965     & objf_gencost(igen_ers,bi,bj)+junk*junk*junkweight
966     if ( junkweight .ne. 0. )
967     & num_gencost(igen_ers,bi,bj) =
968     & num_gencost(igen_ers,bi,bj) + 1. _d 0
969     #endif
970     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
971     c-- The array gfoobs contains SSH anomalies.
972     junkweight = mean_psMgfoobs_MSK(i,j,bi,bj)
973     & *gencost_weight(i,j,bi,bj,igen_gfo)
974     & *gfomask(i,j,bi,bj)
975     & *cosphi(i,j,bi,bj)
976     junk = anom_psMgfoobs(i,j,bi,bj)
977     objf_gencost(igen_gfo,bi,bj) =
978     & objf_gencost(igen_gfo,bi,bj)+junk*junk*junkweight
979     if ( junkweight .ne. 0. )
980     & num_gencost(igen_gfo,bi,bj) =
981     & num_gencost(igen_gfo,bi,bj) + 1. _d 0
982     #endif
983     enddo
984     enddo
985     enddo
986     enddo
987    
988     enddo
989    
990    
991     #endif /* ifdef ALLOW_GENCOST_FREEFORM */
992     #endif /* ifdef ALLOW_GENCOST_CONTRIBUTION */
993     #endif /* ifdef ALLOW_SMOOTH */
994     #endif /* ifdef ALLOW_SSH_COST_CONTRIBUTION */
995    
996     end
  ViewVC Help
Powered by ViewVC 1.1.22