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

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

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


Revision 1.38 - (hide annotations) (download)
Mon Oct 6 22:57:36 2014 UTC (9 years, 7 months ago) by gforget
Branch: MAIN
Changes since 1.37: +7 -4 lines
- cost_gencost_boxmean.F, cost_gencost_bpv4.F, ecco_cost_init_fixed.F,
  ecco_cost_weights.F : treat cases of missing files at run time.
- cost_gencost_sshv4.F : rename etabar as etaday to reflect
  the required averaging frequency, move the test for etaday (now
  used instead of psbar) definition to ecco_check.F
- ecco_check.F : if no gencost term properly defines etaday (e.g.
  due to missing files) then turn off sshv4 cost all together.

1 gforget 1.38 C $Header: /u/gcmpack/MITgcm/pkg/ecco/ecco_cost_init_fixed.F,v 1.37 2014/10/03 15:01:07 gforget Exp $
2 jmc 1.11 C $Name: $
3 heimbach 1.1
4 jmc 1.28 #include "ECCO_OPTIONS.h"
5 heimbach 1.2 #include "AD_CONFIG.h"
6 heimbach 1.1
7     subroutine ecco_cost_init_fixed( mythid )
8    
9     c ==================================================================
10     c SUBROUTINE ecco_cost_init_fixed
11     c ==================================================================
12     c
13     c o Set contributions to the cost function and the cost function
14     c itself to zero. The cost function and the individual contribu-
15     c tions are defined in the header file "ecco_cost.h".
16     c
17     c started: Christian Eckert eckert@mit.edu 30-Jun-1999
18     c
19     c changed: Christian Eckert eckert@mit.edu 25-Feb-2000
20     c
21     c - Restructured the code in order to create a package
22     c for the MITgcmUV.
23     c
24     c changed: Ralf Giering 18-Jan-2001
25     c
26     c - move namelist reading to cost_readparms.F
27     c
28     c ==================================================================
29     c SUBROUTINE ecco_cost_init_fixed
30     c ==================================================================
31    
32     implicit none
33    
34     c == global variables ==
35    
36     #include "EEPARAMS.h"
37     #include "SIZE.h"
38     #include "GRID.h"
39     #include "PARAMS.h"
40    
41 heimbach 1.3 #ifdef ALLOW_CAL
42 heimbach 1.1 #include "cal.h"
43 heimbach 1.3 #endif
44 gforget 1.34 #ifdef ALLOW_ECCO
45 gforget 1.35 c# ifdef ALLOW_OLD_ESTIM_CODES
46     # include "ecco_cost.h"
47     c# else
48     c# include "ecco.h"
49     c# endif
50 gforget 1.34 #endif
51     #ifdef ALLOW_CTRL
52     # include "optim.h"
53     #endif
54 heimbach 1.1
55     c == routine arguments ==
56    
57     integer mythid
58    
59     c == local variables ==
60    
61 gforget 1.35 integer k
62     logical exst
63     c#ifdef ALLOW_OLD_ESTIM_CODES
64     _RL dummy
65 heimbach 1.1 integer tempDate1(4)
66     integer tempDate2(4)
67 heimbach 1.9 integer gwunit
68 gforget 1.38 integer il, ilo,ihi
69 gforget 1.35 character*(max_len_mbuf) msgbuf
70     integer irec
71 heimbach 1.9 _RL missingObsFlag
72     PARAMETER ( missingObsFlag = 1. _d 23 )
73 gforget 1.35 c#endif
74     #ifdef ALLOW_GENCOST_CONTRIBUTION
75     integer bi,bj
76     integer i,j,kk,k2
77     #endif
78 heimbach 1.1
79     c == external functions ==
80    
81 heimbach 1.8 integer cal_IntYears
82     external cal_IntYears
83 heimbach 1.1 integer cal_IntMonths
84     external cal_IntMonths
85     integer cal_IntDays
86     external cal_IntDays
87 heimbach 1.9 integer ifnblnk
88     external ifnblnk
89     integer ilnblnk
90     external ilnblnk
91 heimbach 1.1
92     c == end of interface ==
93    
94 gforget 1.34 #ifdef ALLOW_CTRL
95     eccoiter=optimcycle
96     #else
97     eccoiter=0
98     #endif
99    
100 heimbach 1.1 #ifdef ALLOW_CAL
101    
102     c-- The number of monthly and daily averages generated by the
103     c-- current model integration.
104 heimbach 1.8 nyearsrec = cal_IntYears( mythid )
105 heimbach 1.1 nmonsrec = cal_IntMonths( mythid )
106     ndaysrec = cal_IntDays( mythid )
107    
108     _BEGIN_MASTER( myThid )
109    
110 gforget 1.35 c#ifdef ALLOW_OLD_ESTIM_CODES
111    
112 heimbach 1.1 c-- Get the complete dates of the ...
113     c-- ... TMI data.
114 mlosch 1.18 if ( tmidatfile .ne. ' ' )
115     & call cal_FullDate( tmistartdate1, tmistartdate2,
116     & tmistartdate, mythid )
117 heimbach 1.1 c-- ... SST data.
118 mlosch 1.18 if ( sstdatfile .ne. ' ' )
119     & call cal_FullDate( sststartdate1, sststartdate2,
120     & sststartdate, mythid )
121 heimbach 1.1 c-- ... SSS data.
122 mlosch 1.18 if ( sssdatfile .ne. ' ' )
123     & call cal_FullDate( sssstartdate1, sssstartdate2,
124     & sssstartdate, mythid )
125 heimbach 1.10 c-- ... BP data.
126 mlosch 1.18 if ( bpdatfile .ne. ' ' )
127     & call cal_FullDate( bpstartdate1, bpstartdate2,
128     & bpstartdate, mythid )
129 heimbach 1.26 c-- ... IES data.
130     if ( iesdatfile .ne. ' ' )
131     & call cal_FullDate( iesstartdate1, iesstartdate2,
132     & iesstartdate, mythid )
133 atn 1.33 #ifdef ALLOW_SSH_MEAN_COST_CONTRIBUTION
134     c-- ... mdt data.
135     if ( mdtdatfile .ne. ' ' )
136     & call cal_FullDate( mdtstartdate1, mdtstartdate2,
137     & mdtstartdate, mythid )
138     c-- ... mdt data.
139     if ( mdtdatfile .ne. ' ' )
140     & call cal_FullDate( mdtenddate1, mdtenddate2,
141     & mdtenddate, mythid )
142     #endif /* ALLOW_SSH_MEAN_COST_CONTRIBUTION */
143 heimbach 1.1 c-- ... T/P data.
144 mlosch 1.18 if ( topexfile .ne. ' ' )
145     & call cal_FullDate( topexstartdate1, topexstartdate2,
146 heimbach 1.1 & topexstartdate, mythid )
147     c-- ... ERS data.
148 mlosch 1.18 if ( ersfile .ne. ' ' )
149     & call cal_FullDate( ersstartdate1, ersstartdate2,
150     & ersstartdate, mythid )
151 heimbach 1.4 c-- ... GFO data.
152 mlosch 1.18 if ( gfofile .ne. ' ' )
153     & call cal_FullDate( gfostartdate1, gfostartdate2,
154     & gfostartdate, mythid )
155 heimbach 1.1 c-- ... SCAT data.
156 mlosch 1.18 if ( scatxdatfile .ne. ' ' )
157     & call cal_FullDate( scatstartdate1, scatstartdate2,
158     & scatxstartdate, mythid )
159     if ( scatydatfile .ne. ' ' )
160     & call cal_FullDate( scatstartdate1, scatstartdate2,
161     & scatystartdate, mythid )
162 heimbach 1.1 c-- ... ARGO data.
163 mlosch 1.18 if ( argotfile .ne. ' ' )
164     & call cal_FullDate( argotstartdate1, argotstartdate2,
165 heimbach 1.1 & argotstartdate, mythid )
166 mlosch 1.18 if ( argosfile .ne. ' ' )
167     & call cal_FullDate( argosstartdate1, argotstartdate2,
168 heimbach 1.1 & argosstartdate, mythid )
169 gforget 1.35 c#endif /* ALLOW_OLD_ESTIM_CODES */
170 heimbach 1.1
171 gforget 1.23 #ifdef ALLOW_GENCOST_CONTRIBUTION
172     do k = 1, NGENCOST
173 gforget 1.34
174 gforget 1.37 c-- skip averaging when several cost terms use the
175     c same barfile or when barfile is undefined
176     gencost_barskip(k)=.FALSE.
177     if ( gencost_barfile(k).EQ.' ' )
178     & gencost_barskip(k)=.TRUE.
179     do k2 = 1,k-1
180     if ( gencost_barfile(k2).EQ.gencost_barfile(k) )
181     & gencost_barskip(k)=.TRUE.
182     enddo
183    
184 gforget 1.35 c-- set time averaging parameters
185 gforget 1.34 if ( (using_gencost(k)).AND.( (gencost_flag(k).GE.1).OR.
186     & (gencost_avgperiod(k).NE.' ') ) ) then
187 gforget 1.23 if ( gencost_avgperiod(k) .EQ. 'day' .OR.
188     & gencost_avgperiod(k) .EQ. 'DAY' ) then
189     gencost_nrec(k) = ndaysrec
190     gencost_period(k) = 86400.
191     else if ( gencost_avgperiod(k) .EQ. 'month' .OR.
192     & gencost_avgperiod(k) .EQ. 'MONTH' ) then
193     gencost_nrec(k) =nmonsrec
194     gencost_period(k) = 0.
195     else if ( gencost_avgperiod(k) .EQ. 'year' .OR.
196     & gencost_avgperiod(k) .EQ. 'YEAR' ) then
197 heimbach 1.24 STOP
198     & 'ecco_cost_init_fixed: yearly data not yet implemented'
199 gforget 1.23 else
200 jmc 1.28 STOP
201 heimbach 1.24 & 'ecco_cost_init_fixed: gencost_avgperiod wrongly specified'
202 gforget 1.23 endif
203 gforget 1.36 if ( gencost_nrecperiod(k).EQ.0)
204     & gencost_nrecperiod(k)=gencost_nrec(k)
205 heimbach 1.24 endif
206 gforget 1.34
207 gforget 1.35 c-- set observation start/enddate
208 gforget 1.34 if (gencost_startdate1(k).GT.0) then
209     call cal_FullDate(
210     & gencost_startdate1(k), gencost_startdate2(k),
211     & gencost_startdate(1,k), mythid )
212     else
213     call cal_CopyDate(modelStartDate,
214     & gencost_startdate(1,k),mythid)
215     endif
216    
217     if (gencost_enddate1(k).GT.0) then
218     call cal_FullDate(
219     & gencost_enddate1(k), gencost_enddate2(k),
220     & gencost_enddate(1,k), mythid )
221     else
222     call cal_CopyDate(modelEndDate,
223     & gencost_enddate(1,k),mythid)
224     endif
225    
226 gforget 1.35 c-- set weights
227     if ( gencost_errfile(k) .NE. ' ' ) then
228    
229     kk=gencost_pointer3d(k)
230 gforget 1.38 il = ilnblnk(gencost_errfile(k))
231     inquire( file=gencost_errfile(k)(1:il), exist=exst )
232     if (.NOT.exst) using_gencost(k)=.FALSE.
233 gforget 1.35
234     if ( (.NOT. gencost_timevaryweight(k)).AND.
235 gforget 1.38 & (.NOT.gencost_is3d(k)).AND.(exst) ) then
236 gforget 1.35
237     call mdsreadfield( gencost_errfile(k), cost_iprec,
238     & cost_yftype, 1, gencost_weight(:,:,1,1,k), 1, mythid )
239    
240     DO bj=myByLo(myThid),myByHi(myThid)
241     DO bi=myBxLo(myThid),myBxHi(myThid)
242     DO j = 1-Oly,sNy+Oly
243     DO i = 1-Olx,sNx+Olx
244     c-- initialize weight to 0 if needed
245     if( gencost_timevaryweight(k)) then
246     gencost_weight(i,j,bi,bj,k) = 0. _d 0
247     endif
248     c-- Test for missing values.
249     if (gencost_weight(i,j,bi,bj,k) .lt. -9900.) then
250     gencost_weight(i,j,bi,bj,k) = 0. _d 0
251     endif
252     c-- Convert to weight
253     if (gencost_weight(i,j,bi,bj,k) .ne. 0.) then
254     gencost_weight(i,j,bi,bj,k) =
255     & 1./gencost_weight(i,j,bi,bj,k)/
256     & gencost_weight(i,j,bi,bj,k)
257     endif
258     enddo
259     enddo
260     enddo
261     enddo
262    
263     #ifdef ALLOW_GENCOST3D
264     elseif ( (.NOT. gencost_timevaryweight(k)).AND.
265 gforget 1.38 & (gencost_is3d(k)).AND.(kk.LE.NGENCOST3D).AND.(exst) ) then
266 gforget 1.35
267     call mdsreadfield( gencost_errfile(k), cost_iprec,
268     & cost_yftype, nr, gencost_wei3d(:,:,1,1,1,kk), 1, mythid )
269    
270     DO bj=myByLo(myThid),myByHi(myThid)
271     DO bi=myBxLo(myThid),myBxHi(myThid)
272     DO j = 1-Oly,sNy+Oly
273     DO i = 1-Olx,sNx+Olx
274     DO k2 = 1,nr
275     c-- initialize weight to 0 if needed
276     if( gencost_timevaryweight(k)) then
277     gencost_wei3d(i,j,k2,bi,bj,kk) = 0. _d 0
278     endif
279     c-- Test for missing values.
280     if (gencost_wei3d(i,j,k2,bi,bj,kk) .lt. -9900.) then
281     gencost_wei3d(i,j,k2,bi,bj,kk) = 0. _d 0
282     endif
283     c-- Convert to weight
284     if (gencost_wei3d(i,j,k2,bi,bj,kk) .ne. 0.) then
285     gencost_wei3d(i,j,k2,bi,bj,kk) =
286     & 1./gencost_wei3d(i,j,k2,bi,bj,kk)/
287     & gencost_wei3d(i,j,k2,bi,bj,kk)
288     endif
289     enddo
290     enddo
291     enddo
292     enddo
293     enddo
294     #endif /* ALLOW_GENCOST3D */
295    
296     endif !if ( (.NOT. gencost_timevaryweight(k)).AND.
297     endif !if ( gencost_errfile(k) .NE. ' ' ) then
298     enddo !do k = 1, NGENCOST
299 gforget 1.23 #endif /* ALLOW_GENCOST_CONTRIBUTION */
300    
301 heimbach 1.1 _END_MASTER( mythid )
302    
303     #endif /* ALLOW_CAL */
304    
305 heimbach 1.14 call ecco_check( myThid )
306    
307 heimbach 1.1 c-- Get the weights that are to be used for the individual cost
308     c-- function contributions.
309     call ecco_cost_weights( mythid )
310    
311     c-- Initialise adjoint of monthly mean files calculated
312     c-- in cost_averagesfields (and their ad...).
313     cph(
314     cph The following init. shoud not be applied if in the middle
315     cph of a divided adjoint run
316     cph)
317     #ifndef ALLOW_TANGENTLINEAR_RUN
318     cph!!! and I think it needs to be seen by TAF
319     cph!!! for repeated TLM runs
320     cph!!!
321     inquire( file='costfinal', exist=exst )
322     if ( .NOT. exst) then
323     call ecco_cost_init_barfiles( mythid )
324     endif
325     #endif
326    
327 gforget 1.35 c#ifdef ALLOW_OLD_ESTIM_CODES
328    
329 heimbach 1.9 #ifdef ALLOW_TRANSPORT_COST_CONTRIBUTION
330 heimbach 1.12 do irec = 1, ndaysrec
331     wtransp(irec) = 0. _d 0
332     transpobs(irec) = 0. _d 0
333     enddo
334    
335     if ( costTranspDataFile .NE. ' ' ) then
336 heimbach 1.9 _BEGIN_MASTER(myThid)
337     ilo = ifnblnk(costTranspDataFile)
338     ihi = ilnblnk(costTranspDataFile)
339     CALL OPEN_COPY_DATA_FILE(
340 jmc 1.11 I costTranspDataFile(ilo:ihi),
341 heimbach 1.9 I 'ECCO_COST_INIT_FIXED',
342     O gwunit,
343     I myThid )
344     do irec = 1, ndaysrec
345     c-- read daily transport time series
346     c-- 1st: transport in m/s
347     c-- 2nd: date in YYYYMMDD
348     c-- 3rd: uncertainty in m/s
349     read(gwunit,*) transpobs(irec), dummy, wtransp(irec)
350     c-- convert std.dev. to weight
351     if ( wtransp(irec) .NE. 0. )
352 heimbach 1.12 & wtransp(irec) =1.0/(wtransp(irec)*wtransp(irec))
353 heimbach 1.9 c-- set weight to zero for missing values
354     if ( transpobs(irec) .EQ. missingObsFlag )
355 heimbach 1.12 & wtransp(irec) = 0. _d 0
356 heimbach 1.9 enddo
357     _END_MASTER(myThid)
358     _BARRIER
359 heimbach 1.12 endif
360 heimbach 1.9 #endif /* ALLOW_TRANSPORT_COST_CONTRIBUTION */
361    
362 gforget 1.29 #ifdef ALLOW_NEW_SSH_COST
363 heimbach 1.12
364     c-- Read flags for picking SSH time averages
365     do irec = 1, ndaysrec
366     tpTimeMask(irec) = 1. _d 0
367     ersTimeMask(irec) = 1. _d 0
368     gfoTimeMask(irec) = 1. _d 0
369     enddo
370     c
371 heimbach 1.13 _BEGIN_MASTER(myThid)
372 heimbach 1.12 c
373     #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
374     if ( tpTimeMaskFile .NE. ' ' ) then
375     ilo = ifnblnk(tpTimeMaskFile)
376     ihi = ilnblnk(tpTimeMaskFile)
377     CALL OPEN_COPY_DATA_FILE(
378     I tpTimeMaskFile(ilo:ihi),
379     I 'cost_ssh tp',
380     O gwunit,
381     I myThid )
382     do irec = 1, ndaysrec
383     read(gwunit,*) tpTimeMask(irec)
384     enddo
385     endif
386     #endif
387     c
388     #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
389     if ( ersTimeMaskFile .NE. ' ' ) then
390     ilo = ifnblnk(ersTimeMaskFile)
391     ihi = ilnblnk(ersTimeMaskFile)
392     CALL OPEN_COPY_DATA_FILE(
393     I ersTimeMaskFile(ilo:ihi),
394     I 'cost_ssh ers',
395     O gwunit,
396     I myThid )
397     do irec = 1, ndaysrec
398     read(gwunit,*) ersTimeMask(irec)
399     enddo
400     endif
401     #endif
402     c
403     #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
404     if ( gfoTimeMaskFile .NE. ' ' ) then
405     ilo = ifnblnk(gfoTimeMaskFile)
406     ihi = ilnblnk(gfoTimeMaskFile)
407     CALL OPEN_COPY_DATA_FILE(
408     I gfoTimeMaskFile(ilo:ihi),
409     I 'cost_ssh gfo',
410     O gwunit,
411     I myThid )
412     do irec = 1, ndaysrec
413     read(gwunit,*) gfoTimeMask(irec)
414     enddo
415     endif
416     #endif
417     c
418 heimbach 1.15 do irec = 1, ndaysrec
419 jmc 1.28 if (
420 heimbach 1.15 & ( tpTimeMask(irec).NE.0. .AND. tpTimeMask(irec).NE.1. ) .OR.
421     & ( ersTimeMask(irec).NE.0. .AND. ersTimeMask(irec).NE.1. ) .OR.
422     & ( ersTimeMask(irec).NE.0. .AND. ersTimeMask(irec).NE.1. ) )
423     & then
424 jmc 1.17 WRITE(msgBuf,'(2A,I10)')
425 heimbach 1.15 & 'ecco_cost_init_fixed: (SSH)TimeMask not 0. or 1. ',
426     & 'for irec (=day) ', irec
427 gforget 1.31 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
428 heimbach 1.15 & SQUEEZE_RIGHT , myThid )
429     CALL PRINT_ERROR( msgBuf , myThid )
430     STOP 'ABNORMAL END: S/R ECCO_COST_INIT_FIXED'
431     endif
432     enddo
433     c
434 heimbach 1.12 _END_MASTER(myThid)
435     _BARRIER
436 gforget 1.29 #endif /* ALLOW_NEW_SSH_COST */
437 heimbach 1.12
438 gforget 1.35 c#endif /* ALLOW_OLD_ESTIM_CODES */
439    
440 jmc 1.20 c-- Summarize the cost function setup.
441 heimbach 1.1 _BEGIN_MASTER( mythid )
442 gforget 1.36 call ecco_summary( mythid )
443 heimbach 1.1 call ecco_cost_summary( mythid )
444     _END_MASTER( mythid )
445    
446     _BARRIER
447    
448     end

  ViewVC Help
Powered by ViewVC 1.1.22