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

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

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


Revision 1.47 - (show annotations) (download)
Mon Apr 3 23:16:38 2017 UTC (7 years ago) by ou.wang
Branch: MAIN
CVS Tags: checkpoint66g, checkpoint66f, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, HEAD
Changes since 1.46: +38 -1 lines
Add and update some files related to ECCO v4 release 3. The main changes are
 adding cost_gencost_glbmean.F to compute the costs for global mean time-series of OBP
 & SSH, splitting the profile cost to sample-mean and -anomaly parts,
 updating the cost formula for the proxy sea-ice costs, adding initial velocity
 as part of the control variables, and applying the correction of global
 mean steric height change to OBP.
- cost_gencost_glbmean.F:    compute cost for global mean OBP & SSH time-series
- cost_profiles.F:           add a CPP option (ALLOW_PROFILES_SAMPLESPLIT_COST)
                              that if defined splits the profile cost to
                              sample-mean and sample-anomaly costs.
- cost_gencost_seaicev4.F:   add a CPP option (SEAICECOST_JPL) that if defined
                             reads in the weights for the proxy sea-ice costs.
- cost_gencost_curstomize.F: if ALLOW_PSBAR_STERIC is defined, correct OBP with
                              the global mean steric ssh change to be consistent
                              with SSH.
- ecco_phys.F:               include sterGloH in a common block to be used in
                              cost_gencost_customize.F to correct OBP when
                              ALLOW_PSBAR_STERIC is defined.
- ctrl_map_ini_genarr.F:     add CPP options to include initial UV as the control
                              variables.
- cost_gencost_all.F:        changes related to cost for the global mean time-series.
- cost_gencost_bpv4.F:       minor change related to cost for the global mean time-series.
- ecco_readparms.F:          changes related to cost for the global mean time-series.
- ecco_ad_diff.list:         include the new added cost_gencost_glbmean.F.
- ctrl_get_gen.F:            check if the first seven characters of xx_tauu_file &
                              xx_tauv_filei matches "xx_tauu" or "xx_tauv".
- ecco_cost_final.F:         add modification related to the sample-mean profile cost.
- ecco_cost_init_fixed.F:    initialize the data and weights for the global mean costs
                              for OBP & SSH, if ALLOW_GENCOST_1D is defined.
- profiles_readvector.F:     read in climatology field in the profile file, and mask
                              out the weight where climatology data is flagged. Active
                              if ALLOW_PROFILES_CLIMMASK is defined.
- PROFILES_SIZE.h:           changes related to the split of profile cost to sample-mean
                              and sample-anomaly.
- profiles.h:                changes related to the split of profile cost to sample-mean
                              and sample-anomaly.
- ecco.h:                    changes related to the global mean costs
- profiles_readparms.F:      changes related to the profile cost splitting and masking
                              by climatology data.
- profiles_init_fixed.F:     changes related to the profile cost splitting.
- profiles_init_varia.F:     changes related to the profile cost splitting.

1 C $Header: /u/gcmpack/MITgcm/pkg/ecco/ecco_cost_init_fixed.F,v 1.46 2016/09/20 17:15:13 gforget Exp $
2 C $Name: $
3
4 #include "ECCO_OPTIONS.h"
5 #include "AD_CONFIG.h"
6
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 #ifdef ALLOW_CAL
42 #include "cal.h"
43 #endif
44 #ifdef ALLOW_ECCO
45 # ifdef ECCO_CTRL_DEPRECATED
46 # include "ecco_cost.h"
47 # else
48 # include "ecco.h"
49 # endif
50 #endif
51 #ifdef ALLOW_CTRL
52 # include "optim.h"
53 #endif
54
55 c == routine arguments ==
56
57 integer mythid
58
59 c == local variables ==
60
61 integer k
62 logical exst
63 _RL dummy
64 integer tempDate1(4)
65 integer tempDate2(4)
66 integer gwunit
67 integer il, ilo,ihi
68 character*(max_len_mbuf) msgbuf
69 integer irec
70 _RL missingObsFlag
71 PARAMETER ( missingObsFlag = 1. _d 23 )
72 #ifdef ALLOW_GENCOST_CONTRIBUTION
73 integer bi,bj
74 integer i,j,kk,k2
75 #endif
76
77 c == external functions ==
78
79 integer cal_IntYears
80 external cal_IntYears
81 integer cal_IntMonths
82 external cal_IntMonths
83 integer cal_IntDays
84 external cal_IntDays
85 integer ifnblnk
86 external ifnblnk
87 integer ilnblnk
88 external ilnblnk
89
90 c == end of interface ==
91
92 #ifdef ALLOW_CTRL
93 eccoiter=optimcycle
94 #else
95 eccoiter=0
96 #endif
97
98 #ifdef ALLOW_CAL
99
100 c-- The number of monthly and daily averages generated by the
101 c-- current model integration.
102 nyearsrec = cal_IntYears( mythid )
103 nmonsrec = cal_IntMonths( mythid )
104 ndaysrec = cal_IntDays( mythid )
105
106 _BEGIN_MASTER( myThid )
107
108 #ifdef ECCO_CTRL_DEPRECATED
109
110 c-- Get the complete dates of the ...
111 c-- ... TMI data.
112 if ( tmidatfile .ne. ' ' )
113 & call cal_FullDate( tmistartdate1, tmistartdate2,
114 & tmistartdate, mythid )
115 c-- ... SST data.
116 if ( sstdatfile .ne. ' ' )
117 & call cal_FullDate( sststartdate1, sststartdate2,
118 & sststartdate, mythid )
119 c-- ... SSS data.
120 if ( sssdatfile .ne. ' ' )
121 & call cal_FullDate( sssstartdate1, sssstartdate2,
122 & sssstartdate, mythid )
123 c-- ... BP data.
124 if ( bpdatfile .ne. ' ' )
125 & call cal_FullDate( bpstartdate1, bpstartdate2,
126 & bpstartdate, mythid )
127 c-- ... IES data.
128 if ( iesdatfile .ne. ' ' )
129 & call cal_FullDate( iesstartdate1, iesstartdate2,
130 & iesstartdate, mythid )
131 #ifdef ALLOW_SSH_MEAN_COST_CONTRIBUTION
132 c-- ... mdt data.
133 if ( mdtdatfile .ne. ' ' )
134 & call cal_FullDate( mdtstartdate1, mdtstartdate2,
135 & mdtstartdate, mythid )
136 c-- ... mdt data.
137 if ( mdtdatfile .ne. ' ' )
138 & call cal_FullDate( mdtenddate1, mdtenddate2,
139 & mdtenddate, mythid )
140 #endif /* ALLOW_SSH_MEAN_COST_CONTRIBUTION */
141 c-- ... T/P data.
142 if ( topexfile .ne. ' ' )
143 & call cal_FullDate( topexstartdate1, topexstartdate2,
144 & topexstartdate, mythid )
145 c-- ... ERS data.
146 if ( ersfile .ne. ' ' )
147 & call cal_FullDate( ersstartdate1, ersstartdate2,
148 & ersstartdate, mythid )
149 c-- ... GFO data.
150 if ( gfofile .ne. ' ' )
151 & call cal_FullDate( gfostartdate1, gfostartdate2,
152 & gfostartdate, mythid )
153 c-- ... SCAT data.
154 if ( scatxdatfile .ne. ' ' )
155 & call cal_FullDate( scatstartdate1, scatstartdate2,
156 & scatxstartdate, mythid )
157 if ( scatydatfile .ne. ' ' )
158 & call cal_FullDate( scatstartdate1, scatstartdate2,
159 & scatystartdate, mythid )
160 c-- ... ARGO data.
161 if ( argotfile .ne. ' ' )
162 & call cal_FullDate( argotstartdate1, argotstartdate2,
163 & argotstartdate, mythid )
164 if ( argosfile .ne. ' ' )
165 & call cal_FullDate( argosstartdate1, argotstartdate2,
166 & argosstartdate, mythid )
167 #endif /* ECCO_CTRL_DEPRECATED */
168
169 #ifdef ALLOW_GENCOST_CONTRIBUTION
170 do k = 1, NGENCOST
171
172 c-- skip averaging when several cost terms use the
173 c same barfile or when barfile is undefined
174 gencost_barskip(k)=.FALSE.
175 if ( gencost_barfile(k).EQ.' ' )
176 & gencost_barskip(k)=.TRUE.
177 do k2 = 1,k-1
178 if ( gencost_barfile(k2).EQ.gencost_barfile(k) )
179 & gencost_barskip(k)=.TRUE.
180 enddo
181
182 c-- set time averaging parameters
183 if ( (using_gencost(k)).AND.( (gencost_flag(k).GE.1).OR.
184 & (gencost_avgperiod(k).NE.' ') ) ) then
185 if ( gencost_avgperiod(k) .EQ. 'day' .OR.
186 & gencost_avgperiod(k) .EQ. 'DAY' ) then
187 gencost_nrec(k) = ndaysrec
188 gencost_period(k) = 86400.
189 else if ( gencost_avgperiod(k) .EQ. 'month' .OR.
190 & gencost_avgperiod(k) .EQ. 'MONTH' ) then
191 gencost_nrec(k) =nmonsrec
192 gencost_period(k) = 0.
193 else if ( gencost_avgperiod(k) .EQ. 'step' .OR.
194 & gencost_avgperiod(k) .EQ. 'STEP' ) then
195 gencost_nrec(k) =nTimeSteps+1
196 gencost_period(k) = dTtracerLev(1)
197 else if ( gencost_avgperiod(k) .EQ. 'const' .OR.
198 & gencost_avgperiod(k) .EQ. 'CONST' ) then
199 gencost_nrec(k) =1
200 gencost_period(k) = dTtracerLev(1)
201 else if ( gencost_avgperiod(k) .EQ. 'year' .OR.
202 & gencost_avgperiod(k) .EQ. 'YEAR' ) then
203 STOP
204 & 'ecco_cost_init_fixed: yearly data not yet implemented'
205 else
206 STOP
207 & 'ecco_cost_init_fixed: gencost_avgperiod wrongly specified'
208 endif
209 endif
210
211 c-- set observation start/enddate
212 if (gencost_startdate1(k).GT.0) then
213 call cal_FullDate(
214 & gencost_startdate1(k), gencost_startdate2(k),
215 & gencost_startdate(1,k), mythid )
216 else
217 call cal_CopyDate(modelStartDate,
218 & gencost_startdate(1,k),mythid)
219 gencost_startdate1(k)=startdate_1
220 gencost_startdate2(k)=startdate_2
221 endif
222
223 if (gencost_enddate1(k).GT.0) then
224 call cal_FullDate(
225 & gencost_enddate1(k), gencost_enddate2(k),
226 & gencost_enddate(1,k), mythid )
227 else
228 call cal_CopyDate(modelEndDate,
229 & gencost_enddate(1,k),mythid)
230 endif
231
232 #ifdef ALLOW_GENCOST_1D
233 if ( (gencost_name(k).EQ.'gmbp') .OR.
234 & (gencost_name(k).EQ.'gmsl')) then
235 if(gencost_nrec(k).GT.N1DDATA)then
236 WRITE(msgBuf,'(2A,2i8)') 'ecco_cost_init_fixed: ',
237 & 'Increase N1DDATA', N1DDATA, gencost_nrec(k)
238 CALL PRINT_ERROR( msgBuf, myThid)
239 endif
240
241 do irec = 1, N1DDATA
242 gencost_1DDATA(irec,k) = 0. _d 0
243 enddo
244
245 if(gencost_wei1d(k).NE.0. _d 0)then
246 gencost_wei1d(k) = 1. _d 0 / gencost_wei1d(k)
247 & /gencost_wei1d(k)
248 endif
249
250 if ( gencost_datafile(k) .NE. ' ' ) then
251 ilo = ifnblnk(gencost_datafile(k))
252 ihi = ilnblnk(gencost_datafile(k))
253
254 CALL OPEN_COPY_DATA_FILE(
255 I gencost_datafile(k)(ilo:ihi),
256 I 'ECCO_COST_INIT_FIXED: ',
257 O gwunit,
258 I myThid )
259 do irec = 1, gencost_nrec(k)
260 read(gwunit,*) gencost_1DDATA(irec,k)
261 enddo
262 close(gwunit)
263 _BARRIER
264 endif
265
266 endif
267 #endif /* ALLOW_GENCOST_1D */
268
269 enddo !do k = 1, NGENCOST
270 #endif /* ALLOW_GENCOST_CONTRIBUTION */
271
272 _END_MASTER( mythid )
273
274 #endif /* ALLOW_CAL */
275
276 call ecco_check( myThid )
277
278 c-- Get the weights that are to be used for the individual cost
279 c-- function contributions.
280 call ecco_cost_weights( mythid )
281
282 c-- Initialise adjoint of monthly mean files calculated
283 c-- in cost_averagesfields (and their ad...).
284 cph(
285 cph The following init. shoud not be applied if in the middle
286 cph of a divided adjoint run
287 cph)
288 #ifndef ALLOW_TANGENTLINEAR_RUN
289 cph!!! and I think it needs to be seen by TAF
290 cph!!! for repeated TLM runs
291 cph!!!
292 inquire( file='costfinal', exist=exst )
293 if ( .NOT. exst) then
294 call ecco_cost_init_barfiles( mythid )
295 endif
296 #endif
297
298 #ifdef ECCO_CTRL_DEPRECATED
299
300 #ifdef ALLOW_TRANSPORT_COST_CONTRIBUTION
301 do irec = 1, ndaysrec
302 wtransp(irec) = 0. _d 0
303 transpobs(irec) = 0. _d 0
304 enddo
305
306 if ( costTranspDataFile .NE. ' ' ) then
307 _BEGIN_MASTER(myThid)
308 ilo = ifnblnk(costTranspDataFile)
309 ihi = ilnblnk(costTranspDataFile)
310 CALL OPEN_COPY_DATA_FILE(
311 I costTranspDataFile(ilo:ihi),
312 I 'ECCO_COST_INIT_FIXED',
313 O gwunit,
314 I myThid )
315 do irec = 1, ndaysrec
316 c-- read daily transport time series
317 c-- 1st: transport in m/s
318 c-- 2nd: date in YYYYMMDD
319 c-- 3rd: uncertainty in m/s
320 read(gwunit,*) transpobs(irec), dummy, wtransp(irec)
321 c-- convert std.dev. to weight
322 if ( wtransp(irec) .NE. 0. )
323 & wtransp(irec) =1.0/(wtransp(irec)*wtransp(irec))
324 c-- set weight to zero for missing values
325 if ( transpobs(irec) .EQ. missingObsFlag )
326 & wtransp(irec) = 0. _d 0
327 enddo
328 _END_MASTER(myThid)
329 _BARRIER
330 endif
331 #endif /* ALLOW_TRANSPORT_COST_CONTRIBUTION */
332
333 #ifdef ALLOW_NEW_SSH_COST
334
335 c-- Read flags for picking SSH time averages
336 do irec = 1, ndaysrec
337 tpTimeMask(irec) = 1. _d 0
338 ersTimeMask(irec) = 1. _d 0
339 gfoTimeMask(irec) = 1. _d 0
340 enddo
341 c
342 _BEGIN_MASTER(myThid)
343 c
344 #ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
345 if ( tpTimeMaskFile .NE. ' ' ) then
346 ilo = ifnblnk(tpTimeMaskFile)
347 ihi = ilnblnk(tpTimeMaskFile)
348 CALL OPEN_COPY_DATA_FILE(
349 I tpTimeMaskFile(ilo:ihi),
350 I 'cost_ssh tp',
351 O gwunit,
352 I myThid )
353 do irec = 1, ndaysrec
354 read(gwunit,*) tpTimeMask(irec)
355 enddo
356 endif
357 #endif
358 c
359 #ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
360 if ( ersTimeMaskFile .NE. ' ' ) then
361 ilo = ifnblnk(ersTimeMaskFile)
362 ihi = ilnblnk(ersTimeMaskFile)
363 CALL OPEN_COPY_DATA_FILE(
364 I ersTimeMaskFile(ilo:ihi),
365 I 'cost_ssh ers',
366 O gwunit,
367 I myThid )
368 do irec = 1, ndaysrec
369 read(gwunit,*) ersTimeMask(irec)
370 enddo
371 endif
372 #endif
373 c
374 #ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
375 if ( gfoTimeMaskFile .NE. ' ' ) then
376 ilo = ifnblnk(gfoTimeMaskFile)
377 ihi = ilnblnk(gfoTimeMaskFile)
378 CALL OPEN_COPY_DATA_FILE(
379 I gfoTimeMaskFile(ilo:ihi),
380 I 'cost_ssh gfo',
381 O gwunit,
382 I myThid )
383 do irec = 1, ndaysrec
384 read(gwunit,*) gfoTimeMask(irec)
385 enddo
386 endif
387 #endif
388 c
389 do irec = 1, ndaysrec
390 if (
391 & ( tpTimeMask(irec).NE.0. .AND. tpTimeMask(irec).NE.1. ) .OR.
392 & ( ersTimeMask(irec).NE.0. .AND. ersTimeMask(irec).NE.1. ) .OR.
393 & ( ersTimeMask(irec).NE.0. .AND. ersTimeMask(irec).NE.1. ) )
394 & then
395 WRITE(msgBuf,'(2A,I10)')
396 & 'ecco_cost_init_fixed: (SSH)TimeMask not 0. or 1. ',
397 & 'for irec (=day) ', irec
398 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
399 & SQUEEZE_RIGHT , myThid )
400 CALL PRINT_ERROR( msgBuf , myThid )
401 STOP 'ABNORMAL END: S/R ECCO_COST_INIT_FIXED'
402 endif
403 enddo
404 c
405 _END_MASTER(myThid)
406 _BARRIER
407 #endif /* ALLOW_NEW_SSH_COST */
408
409 #endif /* ECCO_CTRL_DEPRECATED */
410
411 c-- Summarize the cost function setup.
412 _BEGIN_MASTER( mythid )
413 call ecco_summary( mythid )
414 call ecco_cost_summary( mythid )
415 _END_MASTER( mythid )
416
417 _BARRIER
418
419 end

  ViewVC Help
Powered by ViewVC 1.1.22