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

Diff of /MITgcm/pkg/ecco/cost_generic.F

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

revision 1.23 by gforget, Sat Oct 4 20:28:20 2014 UTC revision 1.36 by gforget, Thu Oct 29 15:08:26 2015 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "ECCO_OPTIONS.h"  #include "ECCO_OPTIONS.h"
5    
6    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7    CBOP
8    C     !ROUTINE: cost_generic
9    C     !INTERFACE:
10        subroutine cost_generic(        subroutine cost_generic(
11       &     nnzbar, localbarfile, localbar, xx_localbar_mean_dummy,       &     nnzbar, localbarfile, dummy,
12       &     nnzobs, localobsfile, mult_local,       &     nnzobs, localobsfile, localerrfile,
13       &     nrecloc, nrecobs,       &     mult_local, nrecloc, nrecobs,
14       &     localstartdate, localperiod,       &     localstartdate, localperiod,
15       &     ylocmask, localweight,       &     ylocmask, spminloc, spmaxloc, spzeroloc,
16       &     spminloc, spmaxloc, spzeroloc,       &     preproc, preproc_c, preproc_i, preproc_r,
17       &     preproc, posproc, scalefile,       &     posproc, posproc_c, posproc_i, posproc_r,
18       &     outlev, outname,       &     outlev, outname,
19       &     objf_local, num_local,       &     objf_local, num_local,
20       &     myiter, mytime, mythid )       &     myiter, mytime, mythid )
21    
22  c     ==================================================================  C     !DESCRIPTION: \bv
23  c     SUBROUTINE cost_generic  C     Generic routine for evaluating time-dependent
 c     ==================================================================  
 c  
 c     o Generic routine for evaluating time-dependent  
24  c       cost function contribution  c       cost function contribution
25  c  C     \ev
 c     ==================================================================  
 c     SUBROUTINE cost_generic  
 c     ==================================================================  
26    
27    C     !USES:
28        implicit none        implicit none
29    
30  c     == global variables ==  c     == global variables ==
# Line 47  c     == global variables == Line 45  c     == global variables ==
45    
46  c     == routine arguments ==  c     == routine arguments ==
47    
       integer nnzbar, nnzobs  
       integer nrecloc, nrecobs  
48        integer myiter        integer myiter
49        integer mythid        integer mythid
50          integer nnzbar, nnzobs
51          integer nrecloc, nrecobs
52        integer localstartdate(4)        integer localstartdate(4)
53        integer outlev        integer outlev
54          integer preproc_i(NGENPPROC)
55          integer posproc_i(NGENPPROC)
56    
57        _RL localbar   (1-olx:snx+olx,1-oly:sny+oly,nnzbar,nsx,nsy)        _RL objf_local(nsx,nsy)
58        _RL localweight(1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)        _RL num_local(nsx,nsy)
59        _RL xx_localbar_mean_dummy        _RL dummy
60        _RL mult_local        _RL mult_local
61        _RL mytime        _RL mytime
62        _RL localperiod        _RL localperiod
63        _RL spminloc        _RL spminloc
64        _RL spmaxloc        _RL spmaxloc
65        _RL spzeroloc        _RL spzeroloc
66          _RL preproc_r(NGENPPROC)
67          _RL posproc_r(NGENPPROC)
68    
69          character*(1) ylocmask
70          character*(MAX_LEN_FNAM) localbarfile
71          character*(MAX_LEN_FNAM) localobsfile
72          character*(MAX_LEN_FNAM) localerrfile
73          character*(MAX_LEN_FNAM) preproc(NGENPPROC)
74          character*(MAX_LEN_FNAM) preproc_c(NGENPPROC)
75          character*(MAX_LEN_FNAM) posproc(NGENPPROC)
76          character*(MAX_LEN_FNAM) posproc_c(NGENPPROC)
77          character*(MAX_LEN_FNAM) outname
78    
79    #ifdef ALLOW_ECCO
80    
81    c     == local variables ==
82    
83          integer bi,bj,k2
84          integer itlo,ithi
85          integer jtlo,jthi
86          logical domean, doanom
87    
88          _RL localdifmean1  (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
89          _RL localdifmean2  (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
90    
91    CEOP
92    
93          jtlo = mybylo(mythid)
94          jthi = mybyhi(mythid)
95          itlo = mybxlo(mythid)
96          ithi = mybxhi(mythid)
97    
98    c--   Initialise local variables.
99    
100          do bj = jtlo,jthi
101            do bi = itlo,ithi
102              objf_local(bi,bj) = 0. _d 0
103              num_local(bi,bj) = 0. _d 0
104            enddo
105          enddo
106    
107          call ecco_zero(localdifmean1,Nr,zeroRL,myThid)
108          call ecco_zero(localdifmean2,Nr,zeroRL,myThid)
109    
110          domean=.FALSE.
111          doanom=.FALSE.
112          do k2 = 1, NGENPPROC
113              if (preproc(k2).EQ.'mean') domean=.TRUE.
114              if (preproc(k2).EQ.'anom') doanom=.TRUE.
115          enddo
116    
117    C Extra time loop to compute time-mean fields and costs
118          if ( (.NOT. ( localobsfile.EQ.' ' ) )
119         &   .AND. ( domean .OR. doanom ) ) then
120            call cost_genloop(
121         &     localdifmean1,localdifmean2,.FALSE.,
122         &     nnzbar, localbarfile, dummy,
123         &     nnzobs, localobsfile, localerrfile,
124         &     mult_local, nrecloc, nrecobs,
125         &     localstartdate, localperiod,
126         &     ylocmask, spminloc, spmaxloc, spzeroloc,
127         &     preproc, preproc_c, preproc_i, preproc_r,
128         &     posproc, posproc_c, posproc_i, posproc_r,
129         &     outlev, outname,
130         &     objf_local, num_local,
131         &     myiter, mytime, mythid )
132          endif
133    
134          call ecco_zero(localdifmean1,Nr,zeroRL,myThid)
135    
136          if ((.NOT.(localobsfile.EQ.' ')).AND.(.NOT.domean)) then
137            call cost_genloop(
138         &     localdifmean2,localdifmean1,.TRUE.,
139         &     nnzbar, localbarfile, dummy,
140         &     nnzobs, localobsfile, localerrfile,
141         &     mult_local, nrecloc, nrecobs,
142         &     localstartdate, localperiod,
143         &     ylocmask, spminloc, spmaxloc, spzeroloc,
144         &     preproc, preproc_c, preproc_i, preproc_r,
145         &     posproc, posproc_c, posproc_i, posproc_r,
146         &     outlev, outname,
147         &     objf_local, num_local,
148         &     myiter, mytime, mythid )
149          endif
150    
151    #endif /* ALLOW_ECCO */
152    
153          return
154          end
155    
156    C--------------
157    
158          subroutine cost_genloop(
159         &     localdifmeanIn,localdifmeanOut, addVariaCost,
160         &     nnzbar, localbarfile, dummy,
161         &     nnzobs, localobsfile, localerrfile,
162         &     mult_local, nrecloc, nrecobs,
163         &     localstartdate, localperiod,
164         &     ylocmask, spminloc, spmaxloc, spzeroloc,
165         &     preproc, preproc_c, preproc_i, preproc_r,
166         &     posproc, posproc_c, posproc_i, posproc_r,
167         &     outlev, outname,
168         &     objf_local, num_local,
169         &     myiter, mytime, mythid )
170    
171    C     !DESCRIPTION: \bv
172    C     Generic routine for evaluating time-dependent
173    c       cost function contribution
174    C     \ev
175    
176    C     !USES:
177          implicit none
178    
179    c     == global variables ==
180    
181    #include "EEPARAMS.h"
182    #include "SIZE.h"
183    #include "PARAMS.h"
184    #include "GRID.h"
185    #ifdef ALLOW_CAL
186    # include "cal.h"
187    #endif
188    #ifdef ALLOW_ECCO
189    # include "ecco.h"
190    #endif
191    #ifdef ALLOW_SEAICE
192    # include "SEAICE_COST.h"
193    #endif
194    
195    c     == routine arguments ==
196    
197          integer myiter
198          integer mythid
199          integer nnzbar, nnzobs
200          integer nrecloc, nrecobs
201          integer localstartdate(4)
202          integer outlev
203          integer preproc_i(NGENPPROC)
204          integer posproc_i(NGENPPROC)
205    
206        _RL objf_local(nsx,nsy)        _RL objf_local(nsx,nsy)
207        _RL num_local(nsx,nsy)        _RL num_local(nsx,nsy)
208          _RL dummy
209          _RL mult_local
210          _RL mytime
211          _RL localperiod
212          _RL spminloc
213          _RL spmaxloc
214          _RL spzeroloc
215          _RL preproc_r(NGENPPROC)
216          _RL posproc_r(NGENPPROC)
217    
218        character*(1) ylocmask        character*(1) ylocmask
219        character*(MAX_LEN_FNAM) localbarfile        character*(MAX_LEN_FNAM) localbarfile
220        character*(MAX_LEN_FNAM) localobsfile        character*(MAX_LEN_FNAM) localobsfile
221        character*(16) preproc        character*(MAX_LEN_FNAM) localerrfile
222        character*(16) posproc        character*(MAX_LEN_FNAM) preproc(NGENPPROC)
223        character*(MAX_LEN_FNAM) scalefile        character*(MAX_LEN_FNAM) preproc_c(NGENPPROC)
224          character*(MAX_LEN_FNAM) posproc(NGENPPROC)
225          character*(MAX_LEN_FNAM) posproc_c(NGENPPROC)
226        character*(MAX_LEN_FNAM) outname        character*(MAX_LEN_FNAM) outname
227    
228          logical addVariaCost
229          _RL localdifmeanIn   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
230          _RL localdifmeanOut  (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
231    
232  #ifdef ALLOW_ECCO  #ifdef ALLOW_ECCO
233    
234  c     == local variables ==  c     == local variables ==
235    
236        integer bi,bj        integer bi,bj
       integer i,j,k  
237        integer itlo,ithi        integer itlo,ithi
238        integer jtlo,jthi        integer jtlo,jthi
239        integer jmin,jmax        integer irec, jrec
240        integer imin,imax        integer il, k2
       integer irec  
       integer  il  
241        integer localrec, obsrec        integer localrec, obsrec
242          integer nrecloop, nrecclim, k2smooth
243          logical domean, doanom, dovarwei, doclim, dosmooth
244    
245        logical doglobalread        _RL localmask  (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
       logical ladinit  
246    
247        _RL localwww        _RL localbar   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
248        _RL localcost        _RL localweight(1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
249        _RL junk        _RL localtmp   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
250          _RL localobs   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
251        _RL localmask  (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)        _RL localdif   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
252        _RL localobs   (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)        _RL difmask    (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
253        _RL localdif   (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)  
254        _RL difmask  (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)        _RL localdifmsk   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
255          _RL localdifsum   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
256          _RL localdifnum   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
257    
258        character*(128) fname1, fname2, fname3        character*(128) fname1, fname2, fname3
259          character*200 msgbuf
260    
261        logical exst        logical exst
262    
# Line 109  c     == external functions == Line 265  c     == external functions ==
265        integer  ilnblnk        integer  ilnblnk
266        external ilnblnk        external ilnblnk
267    
268  c     == end of interface ==  CEOP
   
       jtlo = mybylo(mythid)  
       jthi = mybyhi(mythid)  
       itlo = mybxlo(mythid)  
       ithi = mybxhi(mythid)  
       jmin = 1  
       jmax = sny  
       imin = 1  
       imax = snx  
   
 c--   Initialise local variables.  
269    
270        localwww = 0. _d 0        call ecco_zero(localbar,Nr,zeroRL,myThid)
271          call ecco_zero(localweight,Nr,zeroRL,myThid)
272        do bj = jtlo,jthi        call ecco_zero(localtmp,Nr,zeroRL,myThid)
273          do bi = itlo,ithi        call ecco_zero(localmask,Nr,zeroRL,myThid)
274            objf_local(bi,bj) = 0. _d 0  
275            num_local(bi,bj) = 0. _d 0        call ecco_zero(localobs,Nr,zeroRL,myThid)
276          enddo        call ecco_zero(localdif,Nr,zeroRL,myThid)
277          call ecco_zero(difmask,Nr,zeroRL,myThid)
278    
279          call ecco_zero(localdifmsk,Nr,zeroRL,myThid)
280          call ecco_zero(localdifsum,Nr,zeroRL,myThid)
281          call ecco_zero(localdifnum,Nr,zeroRL,myThid)
282    
283          domean=.FALSE.
284          doanom=.FALSE.
285          dovarwei=.FALSE.
286          dosmooth=.FALSE.
287          k2smooth=1
288          doclim=.FALSE.
289          nrecclim=nrecloc
290          do k2 = 1, NGENPPROC
291              if (preproc(k2).EQ.'mean') domean=.TRUE.
292              if (preproc(k2).EQ.'anom') doanom=.TRUE.
293              if (preproc(k2).EQ.'variaweight') dovarwei=.TRUE.
294              if (posproc(k2).EQ.'smooth') then
295                dosmooth=.TRUE.
296                k2smooth=k2
297              endif
298              if (preproc(k2).EQ.'clim') then
299                doclim=.TRUE.
300                nrecclim=preproc_i(k2)
301              endif
302        enddo        enddo
       call ecco_zero(localobs,nnzobs,zeroRL,myThid)  
       call ecco_zero(localdif,nnzobs,zeroRL,myThid)  
       call ecco_zero(difmask,nnzobs,zeroRL,myThid)  
303    
304  c--   Assign mask  c--   Assign mask
       do bj = jtlo,jthi  
         do bi = itlo,ithi  
           do k = 1,Nr  
             do j = 1-oly,sny+oly  
               do i = 1-olx,snx+olx  
305        if ( ylocmask .EQ. 'C' .OR. ylocmask .EQ. 'c' ) then        if ( ylocmask .EQ. 'C' .OR. ylocmask .EQ. 'c' ) then
306           localmask(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)          call ecco_cprsrl(maskC,nr,localmask,nr,myThid)
307        elseif ( ylocmask .EQ. 'S' .OR. ylocmask .EQ. 's' ) then        elseif ( ylocmask .EQ. 'S' .OR. ylocmask .EQ. 's' ) then
308           localmask(i,j,k,bi,bj) = maskS(i,j,k,bi,bj)          call ecco_cprsrl(maskS,nr,localmask,nr,myThid)
309        elseif ( ylocmask .EQ. 'W' .OR. ylocmask .EQ. 'w' ) then        elseif ( ylocmask .EQ. 'W' .OR. ylocmask .EQ. 'w' ) then
310           localmask(i,j,k,bi,bj) = maskW(i,j,k,bi,bj)          call ecco_cprsrl(maskW,nr,localmask,nr,myThid)
311        else        else
312           STOP 'cost_generic: wrong ylocmask'           STOP 'cost_generic: wrong ylocmask'
313        endif        endif
               enddo  
             enddo  
           enddo  
         enddo  
       enddo  
   
 c--   First, read tiled data.  
       doglobalread = .false.  
       ladinit      = .false.  
314    
315  c-- reset nrecloc, if needed, according to preproc  c-- set nrecloop to nrecloc
316        if ( preproc .EQ. 'climmon ') nrecloc=MIN(nrecloc,12)        nrecloop=nrecloc
317    
318        if ( .NOT. ( localobsfile.EQ.' ' ) ) then  c-- reset nrecloop, if needed, according to preproc
319          if ( doclim ) nrecloop=MIN(nrecloop,nrecclim)
320    
321  c--   loop over obsfile records  c--   loop over obsfile records
322        do irec = 1, nrecloc        do irec = 1, nrecloop
323    
324    c--     load weights
325            exst=.FALSE.
326            jrec=1
327            if( dovarwei ) jrec = irec
328            call cost_gencal(localbarfile, localerrfile,
329         &     jrec, localstartdate, localperiod, fname1,
330         &     fname3, localrec, obsrec, exst, mythid )
331            call ecco_zero(localweight,nnzobs,zeroRL,myThid)
332            if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) )
333         &  call ecco_readwei(fname3,localweight,localrec,nnzobs,mythid)
334    
335  c--     determine records and file names  c--     determine records and file names
336          exst=.FALSE.          exst=.FALSE.
# Line 174  c--     determine records and file names Line 339  c--     determine records and file names
339       &     fname2, localrec, obsrec, exst, mythid )       &     fname2, localrec, obsrec, exst, mythid )
340    
341  c--     load model average and observed average  c--     load model average and observed average
342          call cost_genread( fname1, localbar, irec, doglobalread,          call ecco_zero(localbar,nnzbar,zeroRL,myThid)
343       &                     ladinit, eccoiter, nnzbar, preproc,          call cost_genread( fname1, localbar, localtmp, irec, nnzbar,
344       &                     mythid, xx_localbar_mean_dummy )       &       nrecloc, preproc, preproc_c, preproc_i, preproc_r,
345         &       dummy, mythid )
346    
347          call ecco_zero(localobs,nnzobs,spzeroloc,myThid)          call ecco_zero(localobs,nnzobs,spzeroloc,myThid)
348          if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) )          if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) )
# Line 184  c--     load model average and observed Line 350  c--     load model average and observed
350       &         localobs, localrec, mythid )       &         localobs, localrec, mythid )
351    
352  c--     Compute masked model-data difference  c--     Compute masked model-data difference
353          call cost_gendif( localbar, nnzbar, localobs, nnzobs,          call ecco_diffmsk( localbar, nnzbar, localobs, nnzobs,
354       &     localmask, spminloc, spmaxloc, spzeroloc,       &     localmask, spminloc, spmaxloc, spzeroloc,
355       &     posproc, scalefile, localdif, difmask, myThid )       &     localdif, difmask, myThid )
356    
357  c--     Compute normalized model-obs cost function          if ( doanom ) call ecco_subtract( localdifmeanIn,
358          do bj = jtlo,jthi       &     nnzobs, localdif, nnzobs, myThid )
           do bi = itlo,ithi  
             localcost    = 0. _d 0  
             do k = 1,nnzobs  
              do j = jmin,jmax  
               do i = imin,imax  
                 localwww  = localweight(i,j,k,bi,bj)  
      &                    * difmask(i,j,k,bi,bj)  
                 junk      = localdif(i,j,k,bi,bj)  
                 localcost = localcost + junk*junk*localwww  
                 if ( localwww .ne. 0. )  
      &               num_local(bi,bj) = num_local(bi,bj) + 1. _d 0  
               enddo  
              enddo  
             enddo  
             objf_local(bi,bj) = objf_local(bi,bj) + localcost  
           enddo  
         enddo  
359    
360            if ( domean.OR.doanom )
361         &    call ecco_addmask(localdif,difmask, nnzobs,localdifsum,
362         &    localdifnum, nnzobs,myThid)
363    
364            if (addVariaCost) then
365    
366    #ifdef ALLOW_SMOOTH
367            if ( useSMOOTH.AND.dosmooth.AND.
368         &     (nnzbar.EQ.1).AND.(nnzobs.EQ.1) )
369         &  call smooth_hetero2d(localdif,maskc,
370         &     posproc_c(k2smooth),posproc_i(k2smooth),mythid)
371    #endif
372    
373    c--     Compute normalized model-obs cost function
374            call ecco_addcost(
375         I                   localdif, localweight, difmask, nnzobs,
376         I                   objf_local, num_local,
377         I                   myThid
378         &                   )
379    
380  c--     output model-data difference to disk  c--     output model-data difference to disk
381          if ( outlev.GT.0 ) then          if ( outlev.GT.0 ) then
# Line 220  c--     output model-data difference to Line 387  c--     output model-data difference to
387       &    WRITE_REC_XYZ_RL( fname3, localdif,irec, eccoiter, mythid )       &    WRITE_REC_XYZ_RL( fname3, localdif,irec, eccoiter, mythid )
388          endif          endif
389    
390            endif
391    
392        enddo        enddo
393  c--   End of loop over obsfile records.  c--   End of loop over obsfile records.
394    
395          call ecco_zero(localdifmeanOut,Nr,zeroRL,myThid)
396          call ecco_cp(localdifsum,nnzobs,localdifmeanOut,nnzobs,myThid)
397          call ecco_divfield(localdifmeanOut,nnzobs,localdifnum,myThid)
398          call ecco_cp(localdifnum,nnzobs,localdifmsk,nnzobs,myThid)
399          call ecco_divfield(localdifmsk,nnzobs,localdifnum,myThid)
400    
401          if ( domean ) then
402    c--     Compute normalized model-obs cost function
403            call ecco_addcost(
404         I       localdifmeanOut, localweight, localdifmsk, nnzobs,
405         I       objf_local, num_local, myThid)
406    
407    c--     output model-data difference to disk
408            if ( outlev.GT.0 ) then
409              il=ilnblnk(outname)
410              write(fname3(1:128),'(2a)') 'misfit_', outname(1:il)
411              if ( nnzobs.EQ.1 ) CALL
412         &    WRITE_REC_XY_RL(fname3,localdifmeanOut,1,eccoiter,mythid)
413              if ( nnzobs.EQ.nr ) CALL
414         &    WRITE_REC_XYZ_RL(fname3,localdifmeanOut,1,eccoiter,mythid)
415            endif
416        endif        endif
417    
418    
419    
420  #endif /* ALLOW_ECCO */  #endif /* ALLOW_ECCO */
421    
422          return
423        end        end
424    
425    
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.36
  ViewVC Help
Powered by ViewVC 1.1.22