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

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

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


Revision 1.9 - (hide annotations) (download)
Sun Jan 26 06:10:10 2014 UTC (10 years, 4 months ago) by atn
Branch: MAIN
CVS Tags: checkpoint64u, checkpoint64t
Changes since 1.8: +6 -3 lines
Replace hardcoded value 300 for eccov4 with general runtime param smooth2Ddiffnbt

1 atn 1.9 C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_gencost_sstv4.F,v 1.8 2012/10/02 01:40:03 gforget Exp $
2 gforget 1.1 C $Name: $
3    
4 jmc 1.7 #include "ECCO_OPTIONS.h"
5 gforget 1.1
6    
7     subroutine cost_gencost_sstv4(
8     I myiter,
9     I mytime,
10     I mythid
11     & )
12    
13     c ==================================================================
14     c SUBROUTINE cost_gencost_sstv4
15     c ==================================================================
16     c
17     c o Evaluate cost function contributions of sea surface temperature.
18     c (Daily Pointwise and then Large Scale)
19     c
20     c started: Gael Forget, Oct-2009
21     c
22     c ==================================================================
23     c SUBROUTINE cost_gencost_sstv4
24     c ==================================================================
25    
26     implicit none
27    
28     c == global variables ==
29    
30     #include "EEPARAMS.h"
31     #include "SIZE.h"
32     #include "PARAMS.h"
33     #include "GRID.h"
34    
35     #include "ecco_cost.h"
36 heimbach 1.6 #include "CTRL_SIZE.h"
37 gforget 1.1 #include "ctrl.h"
38     #include "ctrl_dummy.h"
39     #include "optim.h"
40     #include "DYNVARS.h"
41     #include "cal.h"
42 atn 1.9 #ifdef ALLOW_SMOOTH
43     #include "SMOOTH.h"
44     #endif
45 gforget 1.1
46     c == routine arguments ==
47    
48     integer myiter
49     _RL mytime
50     integer mythid
51    
52     #ifdef ALLOW_DAILYSST_COST_CONTRIBUTION
53     #ifdef ALLOW_SMOOTH
54     #ifdef ALLOW_GENCOST_CONTRIBUTION
55     c == local variables ==
56    
57     integer bi,bj
58     integer i,j,k
59     integer itlo,ithi
60     integer jtlo,jthi
61     integer jmin,jmax
62     integer imin,imax
63     integer irec,jrec,krec
64     integer ilps
65     integer gwunit
66    
67     logical doglobalread
68     logical ladinit
69    
70     _RL anom_sst(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
71     _RL obs_sst (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
72     _RL nb_sst (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
73 jmc 1.4 _RL msk_sst (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
74 gforget 1.1 _RL tmp_sst (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
75     _RL spval
76    
77     integer num_var
78 jmc 1.4
79 gforget 1.1 _RL junk,junkweight
80    
81     integer ndaysave
82     _RL ndaysaveRL
83    
84     character*(80) fname
85     character*(80) fname2
86 gforget 1.3
87     #ifdef ALLOW_ECCO_DEBUG
88     character*(MAX_LEN_MBUF) msgBuf
89     INTEGER ioUnit
90     #endif
91     logical exst
92 gforget 1.1
93     _RL daytime
94     _RL diffsecs
95     integer il, localrec
96     integer dayiter
97     integer daydate(4)
98     integer difftime(4)
99 jmc 1.4 integer tempDate_1
100 gforget 1.1 integer middate(4)
101 gforget 1.3 integer locstartdate(4)
102 gforget 1.1 integer yday, ymod
103     integer md, dd, sd, ld, wd
104     integer mody, modm
105    
106     integer igen_amsre, igen_amsre_lsc
107    
108     c == external functions ==
109    
110     integer ilnblnk
111     external ilnblnk
112    
113     c == end of interface ==
114    
115     jtlo = mybylo(mythid)
116     jthi = mybyhi(mythid)
117     itlo = mybxlo(mythid)
118     ithi = mybxhi(mythid)
119     jmin = 1
120     jmax = sny
121     imin = 1
122     imax = snx
123    
124 gforget 1.3 #ifdef ALLOW_ECCO_DEBUG
125     ioUnit=standardMessageUnit
126     #endif
127    
128     call cal_FullDate(19920101,0,locstartdate,mythid)
129    
130 gforget 1.1 c-- detect the relevant gencost indices
131     igen_amsre=0
132     igen_amsre_lsc=0
133     do k=1,NGENCOST
134     if (gencost_name(k).EQ.'sstv4-amsre') igen_amsre=k
135     if (gencost_name(k).EQ.'sstv4-amsre-lsc') igen_amsre_lsc=k
136     enddo
137    
138     c-- First, read tiled data.
139     doglobalread = .false.
140     ladinit = .false.
141    
142     write(fname(1:80),'(80a)') ' '
143     ilps=ilnblnk( sstbarfile )
144     write(fname(1:80),'(2a,i10.10)')
145     & sstbarfile(1:ilps),'.',optimcycle
146    
147     spval = -9990.
148    
149     cgf =======================================================
150     cgf PART 1: compute smooth SST cost term
151     cgf =======================================================
152    
153    
154     ndaysave=7
155     ndaysaveRL=ndaysave
156    
157     do irec = 1, ndaysrec-ndaysave+1, 7
158    
159     do bj = jtlo,jthi
160     do bi = itlo,ithi
161     do j = jmin,jmax
162     do i = imin,imax
163     anom_sst(i,j,bi,bj) = 0. _d 0
164     obs_sst(i,j,bi,bj) = 0. _d 0
165     nb_sst(i,j,bi,bj) = 0. _d 0
166     msk_sst(i,j,bi,bj) = 0. _d 0
167     enddo
168     enddo
169     enddo
170     enddo
171    
172     c PART 1.1: compute running sample average over ndaysave
173     c ------------------------------------------------------
174    
175     do jrec=1,ndaysave
176    
177     krec=irec+jrec-1
178    
179     c get modeled sst:
180     call active_read_xy( fname, sstbar, krec, doglobalread,
181     & ladinit, optimcycle, mythid,
182     & xx_sstbar_mean_dummy )
183    
184     c get observed sst:
185 gforget 1.3 daytime = FLOAT(secondsperday*(krec-1)) + modelstart
186     dayiter = hoursperday*(krec-1)+modeliter0
187 gforget 1.1 call cal_getdate( dayiter, daytime, daydate, mythid )
188     call cal_convdate( daydate,yday,md,dd,sd,ld,wd,mythid )
189 gforget 1.3 ymod = locstartdate(1)/10000
190     if ( ymod .GE. yday ) then
191 gforget 1.5 middate(1)=1
192     call cal_FullDate(locstartdate(1),0,middate,mythid)
193 gforget 1.1 else
194 gforget 1.5 middate(1)=1
195 jmc 1.4 tempDate_1 = yday*10000+100+1
196     call cal_FullDate( tempDate_1, 0, middate, mythid)
197 gforget 1.1 endif
198     call cal_TimePassed( middate, daydate, difftime, mythid )
199     call cal_ToSeconds( difftime, diffsecs, mythid )
200 gforget 1.3 c localrec = floor(diffsecs/86400.) + 1
201 gforget 1.1 localrec = int(diffsecs/86400.) + 1
202    
203     il=ilnblnk(gencost_datafile(igen_amsre))
204     write(fname2(1:80),'(2a,i4)')
205     & gencost_datafile(igen_amsre)(1:il), '_', yday
206 gforget 1.3 inquire( file=fname2, exist=exst )
207    
208     #ifdef ALLOW_ECCO_DEBUG
209     WRITE(msgBuf,'(A,I4,A,I4,A,I10,A,1PE15.2)') 'sstv4 reading ',
210     & yday,' ',ymod,' ',localrec,' ',diffsecs
211     CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
212     c
213     CALL CAL_PRINTDATE(middate,mythid)
214     CALL CAL_PRINTDATE(daydate,mythid)
215     CALL CAL_PRINTDATE(difftime,mythid)
216 jmc 1.4 #endif
217 gforget 1.1
218 gforget 1.5 if ( ( localrec .GT. 0 ).AND.(diffsecs .GT. 0.d0) ) then
219 gforget 1.1 call mdsreadfield( fname2, cost_iprec, cost_yftype, 1,
220     & tmp_sst, localrec, mythid )
221     else
222     do bj = jtlo,jthi
223     do bi = itlo,ithi
224     do j = jmin,jmax
225     do i = imin,imax
226     tmp_sst(i,j,bi,bj) = spval
227     enddo
228     enddo
229     enddo
230     enddo
231     endif
232    
233     c accumulate obs and misfit:
234     do bj = jtlo,jthi
235     do bi = itlo,ithi
236     do j = jmin,jmax
237     do i = imin,imax
238     if ( (tmp_sst(i,j,bi,bj).GT.spval).AND.
239     & (maskc(i,j,1,bi,bj).EQ.1.) ) then
240     anom_sst(i,j,bi,bj)= anom_sst(i,j,bi,bj)+
241     & sstbar(i,j,bi,bj)-tmp_sst(i,j,bi,bj)
242     obs_sst(i,j,bi,bj)= obs_sst(i,j,bi,bj)+
243 jmc 1.4 & tmp_sst(i,j,bi,bj)
244 gforget 1.1 nb_sst(i,j,bi,bj)=nb_sst(i,j,bi,bj)+1. _d 0
245     endif
246     enddo
247     enddo
248     enddo
249     enddo
250    
251     enddo !do jrec=1,ndaysave
252    
253     c average obs and misfit:
254     do bj = jtlo,jthi
255     do bi = itlo,ithi
256     do j = jmin,jmax
257     do i = imin,imax
258     if ( nb_sst(i,j,bi,bj) .NE. 0. ) then
259 jmc 1.4 obs_sst(i,j,bi,bj) =
260 gforget 1.1 & obs_sst(i,j,bi,bj)/nb_sst(i,j,bi,bj)
261 jmc 1.4 anom_sst(i,j,bi,bj) =
262 gforget 1.1 & anom_sst(i,j,bi,bj)/nb_sst(i,j,bi,bj)
263     msk_sst(i,j,bi,bj) = 1. _d 0
264     endif
265     enddo
266     enddo
267     enddo
268     enddo
269    
270    
271     c PART 1.2: smooth anom_sst in space
272     c ----------------------------------------
273    
274     #ifdef ALLOW_GENCOST_SSTV4_OUTPUT
275     write(fname2(1:80),'(1a)') 'sstdiff_raw'
276     call mdswritefield(fname2,32,.false.,'RL',
277     & 1,anom_sst,irec,1,mythid)
278    
279     write(fname2(1:80),'(1a)') 'sstobs_raw'
280     call mdswritefield(fname2,32,.false.,'RL',
281     & 1,obs_sst,irec,1,mythid)
282     #endif
283    
284     call smooth_hetero2d(anom_sst,maskc,
285 atn 1.9 & gencost_scalefile(igen_amsre_lsc),smooth2Ddiffnbt,mythid)
286 gforget 1.1
287     #ifdef ALLOW_GENCOST_SSTV4_OUTPUT
288     call smooth_hetero2d(obs_sst,maskc,
289 atn 1.9 & gencost_scalefile(igen_amsre_lsc),smooth2Ddiffnbt,mythid)
290 gforget 1.1
291     write(fname2(1:80),'(1a)') 'sstdiff_smooth'
292     call mdswritefield(fname2,32,.false.,'RL',
293 jmc 1.4 & 1,anom_sst,irec,1,mythid)
294 gforget 1.1
295     write(fname2(1:80),'(1a)') 'sstobs_smooth'
296     call mdswritefield(fname2,32,.false.,'RL',
297     & 1,obs_sst,irec,1,mythid)
298 jmc 1.4 #endif
299 gforget 1.1
300     c PART 1.3: compute cost function term
301     c ------------------------------------
302    
303     do bj = jtlo,jthi
304     do bi = itlo,ithi
305     do j = jmin,jmax
306     do i = imin,imax
307     junk = anom_sst(i,j,bi,bj)
308     junkweight = gencost_weight(i,j,bi,bj,igen_amsre_lsc)*
309     & maskc(i,j,1,bi,bj)
310 jmc 1.4 objf_gencost(igen_amsre_lsc,bi,bj) =
311 gforget 1.1 & objf_gencost(igen_amsre_lsc,bi,bj)
312     & +junk*junk*junkweight/ndaysaveRL
313     if ( (junkweight.GT.0.).AND.(nb_sst(i,j,bi,bj).GT.0.) )
314 jmc 1.4 & num_gencost(igen_amsre_lsc,bi,bj) =
315 gforget 1.1 & num_gencost(igen_amsre_lsc,bi,bj) + 1. _d 0 /ndaysaveRL
316     enddo
317     enddo
318     enddo
319     enddo
320    
321     enddo
322    
323     cgf =======================================================
324     cgf PART 2: compute raw SST cost term
325     cgf =======================================================
326    
327     do irec = 1, ndaysrec
328    
329     c get modeled sst:
330 gforget 1.2 call active_read_xy( fname, sstbar, irec, doglobalread,
331 gforget 1.1 & ladinit, optimcycle, mythid,
332 gforget 1.2 & xx_sstbar_mean_dummy )
333 gforget 1.1
334     c get observed sst:
335 gforget 1.3 daytime = FLOAT(secondsperday*(irec-1)) + modelstart
336     dayiter = hoursperday*(irec-1)+modeliter0
337 gforget 1.1 call cal_getdate( dayiter, daytime, daydate, mythid )
338     call cal_convdate( daydate,yday,md,dd,sd,ld,wd,mythid )
339 gforget 1.3 ymod = locstartdate(1)/10000
340     if ( ymod .GE. yday ) then
341 gforget 1.5 middate(1)=1
342     call cal_FullDate(locstartdate(1),0,middate,mythid)
343 gforget 1.1 else
344 gforget 1.5 middate(1)=1
345     tempDate_1 = yday*10000+100+1
346     call cal_FullDate( tempDate_1, 0, middate, mythid)
347 gforget 1.1 endif
348     call cal_TimePassed( middate, daydate, difftime, mythid )
349 jmc 1.4 call cal_ToSeconds( difftime, diffsecs, mythid )
350 gforget 1.3 c localrec = floor(diffsecs/86400.) + 1
351 gforget 1.1 localrec = int(diffsecs/86400.) + 1
352    
353     il=ilnblnk(gencost_datafile(igen_amsre))
354     write(fname2(1:80),'(2a,i4)')
355     & gencost_datafile(igen_amsre)(1:il), '_', yday
356 gforget 1.3 inquire( file=fname2, exist=exst )
357    
358     #ifdef ALLOW_ECCO_DEBUG
359     WRITE(msgBuf,'(A,I4,A,I4,A,I10,A,1PE15.2)') 'sstv4 reading ',
360     & yday,' ',ymod,' ',localrec,' ',diffsecs
361     CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
362     c
363     CALL CAL_PRINTDATE(middate,mythid)
364     CALL CAL_PRINTDATE(daydate,mythid)
365     CALL CAL_PRINTDATE(difftime,mythid)
366     #endif
367 gforget 1.1
368 gforget 1.5 if ( ( localrec .GT. 0 ).AND.(diffsecs .GT. 0.d0) ) then
369 gforget 1.1 call mdsreadfield( fname2, cost_iprec, cost_yftype, 1,
370     & tmp_sst, localrec, mythid )
371     else
372     do bj = jtlo,jthi
373     do bi = itlo,ithi
374     do j = jmin,jmax
375     do i = imin,imax
376     tmp_sst(i,j,bi,bj) = spval
377     enddo
378     enddo
379     enddo
380     enddo
381     endif
382    
383     c compute misfit:
384     do bj = jtlo,jthi
385     do bi = itlo,ithi
386     do j = jmin,jmax
387     do i = imin,imax
388     if ( (tmp_sst(i,j,bi,bj).GT.spval).AND.
389     & (maskc(i,j,1,bi,bj).EQ.1.) ) then
390     anom_sst(i,j,bi,bj) =
391     & sstbar(i,j,bi,bj)-tmp_sst(i,j,bi,bj)
392     msk_sst(i,j,bi,bj) = 1. _d 0
393     else
394     anom_sst(i,j,bi,bj) = 0. _d 0
395     msk_sst(i,j,bi,bj) = 0. _d 0
396     endif
397     enddo
398     enddo
399     enddo
400     enddo
401    
402     #ifdef ALLOW_GENCOST_SSTV4_OUTPUT
403     write(fname2(1:80),'(1a)') 'sstdiff_point'
404     call mdswritefield(fname2,32,.false.,'RL',
405     & 1,anom_sst,irec,1,mythid)
406 jmc 1.4 #endif
407 gforget 1.1
408     c compute cost:
409    
410     do bj = jtlo,jthi
411     do bi = itlo,ithi
412     do j = jmin,jmax
413     do i = imin,imax
414     junk = anom_sst(i,j,bi,bj)
415     junkweight = gencost_weight(i,j,bi,bj,igen_amsre)*
416     & maskc(i,j,1,bi,bj)*msk_sst(i,j,bi,bj)
417     objf_gencost(igen_amsre,bi,bj) =
418     & objf_gencost(igen_amsre,bi,bj)+junk*junk*junkweight
419     if (junkweight.GT.0.)
420     & num_gencost(igen_amsre,bi,bj) =
421     & num_gencost(igen_amsre,bi,bj) + 1. _d 0
422     enddo
423     enddo
424     enddo
425     enddo
426    
427     enddo
428    
429     #endif /* ifdef ALLOW_GENCOST_CONTRIBUTION */
430     #endif /* ifdef ALLOW_SMOOTH */
431     #endif /* ifdef ALLOW_DAILYSST_COST_CONTRIBUTION */
432    
433     end

  ViewVC Help
Powered by ViewVC 1.1.22