/[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.31 by gforget, Fri Jun 5 02:22:48 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, localbar, dummy,
12       &     nnzobs, localobsfile, mult_local,       &     nnzobs, localobsfile, mult_local,
13       &     nrecloc, nrecobs,       &     nrecloc, nrecobs,
14       &     localstartdate, localperiod,       &     localstartdate, localperiod,
15       &     ylocmask, localweight,       &     ylocmask, localweight,
16       &     spminloc, spmaxloc, spzeroloc,       &     spminloc, spmaxloc, spzeroloc,
17       &     preproc, posproc, scalefile,       &     preproc, preproc_c, preproc_i, preproc_r,
18         &     posproc, posproc_c, posproc_i, posproc_r,
19       &     outlev, outname,       &     outlev, outname,
20       &     objf_local, num_local,       &     objf_local, num_local,
21       &     myiter, mytime, mythid )       &     myiter, mytime, mythid )
22    
23  c     ==================================================================  C     !DESCRIPTION: \bv
24  c     SUBROUTINE cost_generic  C     Generic routine for evaluating time-dependent
 c     ==================================================================  
 c  
 c     o Generic routine for evaluating time-dependent  
25  c       cost function contribution  c       cost function contribution
26  c  C     \ev
 c     ==================================================================  
 c     SUBROUTINE cost_generic  
 c     ==================================================================  
27    
28    C     !USES:
29        implicit none        implicit none
30    
31  c     == global variables ==  c     == global variables ==
# Line 47  c     == global variables == Line 46  c     == global variables ==
46    
47  c     == routine arguments ==  c     == routine arguments ==
48    
       integer nnzbar, nnzobs  
       integer nrecloc, nrecobs  
49        integer myiter        integer myiter
50        integer mythid        integer mythid
51          integer nnzbar, nnzobs
52          integer nrecloc, nrecobs
53        integer localstartdate(4)        integer localstartdate(4)
54        integer outlev        integer outlev
55          integer preproc_i(NGENPPROC)
56          integer posproc_i(NGENPPROC)
57    
58          _RL objf_local(nsx,nsy)
59          _RL num_local(nsx,nsy)
60        _RL localbar   (1-olx:snx+olx,1-oly:sny+oly,nnzbar,nsx,nsy)        _RL localbar   (1-olx:snx+olx,1-oly:sny+oly,nnzbar,nsx,nsy)
61        _RL localweight(1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)        _RL localweight(1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)
62        _RL xx_localbar_mean_dummy        _RL dummy
63        _RL mult_local        _RL mult_local
64        _RL mytime        _RL mytime
65        _RL localperiod        _RL localperiod
66        _RL spminloc        _RL spminloc
67        _RL spmaxloc        _RL spmaxloc
68        _RL spzeroloc        _RL spzeroloc
69        _RL objf_local(nsx,nsy)        _RL preproc_r(NGENPPROC)
70        _RL num_local(nsx,nsy)        _RL posproc_r(NGENPPROC)
71    
72        character*(1) ylocmask        character*(1) ylocmask
73        character*(MAX_LEN_FNAM) localbarfile        character*(MAX_LEN_FNAM) localbarfile
74        character*(MAX_LEN_FNAM) localobsfile        character*(MAX_LEN_FNAM) localobsfile
75        character*(16) preproc        character*(MAX_LEN_FNAM) preproc(NGENPPROC)
76        character*(16) posproc        character*(MAX_LEN_FNAM) preproc_c(NGENPPROC)
77        character*(MAX_LEN_FNAM) scalefile        character*(MAX_LEN_FNAM) posproc(NGENPPROC)
78          character*(MAX_LEN_FNAM) posproc_c(NGENPPROC)
79        character*(MAX_LEN_FNAM) outname        character*(MAX_LEN_FNAM) outname
80    
81  #ifdef ALLOW_ECCO  #ifdef ALLOW_ECCO
# Line 79  c     == routine arguments == Line 83  c     == routine arguments ==
83  c     == local variables ==  c     == local variables ==
84    
85        integer bi,bj        integer bi,bj
       integer i,j,k  
86        integer itlo,ithi        integer itlo,ithi
87        integer jtlo,jthi        integer jtlo,jthi
       integer jmin,jmax  
       integer imin,imax  
88        integer irec        integer irec
89        integer  il        integer  il
90        integer localrec, obsrec        integer localrec, obsrec
91          integer nrecloop
92    
93        logical doglobalread        _RL localtmp   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
94        logical ladinit        _RL localmask  (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
95          _RL localobs   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
96        _RL localwww        _RL localdif   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
97        _RL localcost        _RL difmask    (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
98        _RL junk        _RL localbarmean   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
99          _RL localbarnum   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
100        _RL localmask  (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)        _RL localobsmean   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
101        _RL localobs   (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)        _RL localobsnum   (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
       _RL localdif   (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)  
       _RL difmask  (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)  
102    
103        character*(128) fname1, fname2, fname3        character*(128) fname1, fname2, fname3
104          character*200 msgbuf
105    
106        logical exst        logical exst
107    
# Line 109  c     == external functions == Line 110  c     == external functions ==
110        integer  ilnblnk        integer  ilnblnk
111        external ilnblnk        external ilnblnk
112    
113  c     == end of interface ==  CEOP
114    
115        jtlo = mybylo(mythid)        jtlo = mybylo(mythid)
116        jthi = mybyhi(mythid)        jthi = mybyhi(mythid)
117        itlo = mybxlo(mythid)        itlo = mybxlo(mythid)
118        ithi = mybxhi(mythid)        ithi = mybxhi(mythid)
       jmin = 1  
       jmax = sny  
       imin = 1  
       imax = snx  
119    
120  c--   Initialise local variables.  c--   Initialise local variables.
121    
       localwww = 0. _d 0  
   
122        do bj = jtlo,jthi        do bj = jtlo,jthi
123          do bi = itlo,ithi          do bi = itlo,ithi
124            objf_local(bi,bj) = 0. _d 0            objf_local(bi,bj) = 0. _d 0
125            num_local(bi,bj) = 0. _d 0            num_local(bi,bj) = 0. _d 0
126          enddo          enddo
127        enddo        enddo
128        call ecco_zero(localobs,nnzobs,zeroRL,myThid)        call ecco_zero(localtmp,Nr,zeroRL,myThid)
129        call ecco_zero(localdif,nnzobs,zeroRL,myThid)        call ecco_zero(localobs,Nr,zeroRL,myThid)
130        call ecco_zero(difmask,nnzobs,zeroRL,myThid)        call ecco_zero(localdif,Nr,zeroRL,myThid)
131          call ecco_zero(difmask,Nr,zeroRL,myThid)
132    
133          call ecco_zero(localbarmean,Nr,zeroRL,myThid)
134          call ecco_zero(localbarnum,Nr,zeroRL,myThid)
135          call ecco_zero(localobsmean,Nr,zeroRL,myThid)
136          call ecco_zero(localobsnum,Nr,zeroRL,myThid)
137    
138  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  
139        if ( ylocmask .EQ. 'C' .OR. ylocmask .EQ. 'c' ) then        if ( ylocmask .EQ. 'C' .OR. ylocmask .EQ. 'c' ) then
140           localmask(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)          call ecco_cprsrl(maskC,nr,localmask,nr,myThid)
141        elseif ( ylocmask .EQ. 'S' .OR. ylocmask .EQ. 's' ) then        elseif ( ylocmask .EQ. 'S' .OR. ylocmask .EQ. 's' ) then
142           localmask(i,j,k,bi,bj) = maskS(i,j,k,bi,bj)          call ecco_cprsrl(maskS,nr,localmask,nr,myThid)
143        elseif ( ylocmask .EQ. 'W' .OR. ylocmask .EQ. 'w' ) then        elseif ( ylocmask .EQ. 'W' .OR. ylocmask .EQ. 'w' ) then
144           localmask(i,j,k,bi,bj) = maskW(i,j,k,bi,bj)          call ecco_cprsrl(maskW,nr,localmask,nr,myThid)
145        else        else
146           STOP 'cost_generic: wrong ylocmask'           STOP 'cost_generic: wrong ylocmask'
147        endif        endif
               enddo  
             enddo  
           enddo  
         enddo  
       enddo  
   
 c--   First, read tiled data.  
       doglobalread = .false.  
       ladinit      = .false.  
148    
149  c-- reset nrecloc, if needed, according to preproc  c-- set nrecloop to nrecloc
150        if ( preproc .EQ. 'climmon ') nrecloc=MIN(nrecloc,12)        nrecloop=nrecloc
   
       if ( .NOT. ( localobsfile.EQ.' ' ) ) then  
151    
152    C Extra time loop to compute time-mean fields and costs
153          if ( (.NOT. ( localobsfile.EQ.' ' ) )
154         &   .AND. ( preproc(1) .EQ. 'mean' .OR.
155         &           preproc(1) .EQ. 'anom' .OR.
156         &     (preproc(1) .EQ. 'clim' .AND. preproc(2) .EQ. 'anom'))) then
157  c--   loop over obsfile records  c--   loop over obsfile records
158        do irec = 1, nrecloc        do irec = 1, nrecloop
159    
160  c--     determine records and file names  c--     determine records and file names
161          exst=.FALSE.          exst=.FALSE.
# Line 174  c--     determine records and file names Line 164  c--     determine records and file names
164       &     fname2, localrec, obsrec, exst, mythid )       &     fname2, localrec, obsrec, exst, mythid )
165    
166  c--     load model average and observed average  c--     load model average and observed average
167          call cost_genread( fname1, localbar, irec, doglobalread,          call ecco_zero(localbar,nnzbar,zeroRL,myThid)
168       &                     ladinit, eccoiter, nnzbar, preproc,          call cost_genread( fname1, localbar, localtmp, irec, nnzbar,
169       &                     mythid, xx_localbar_mean_dummy )       &       nrecloc, preproc(2), preproc_c(2),
170         &       preproc_i(2), preproc_r(2),
171         &       dummy, mythid )
172    
173          call ecco_zero(localobs,nnzobs,spzeroloc,myThid)          call ecco_zero(localobs,Nr,spzeroloc,myThid)
174          if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) )          if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) )
175       &  call mdsreadfield( fname2, cost_iprec, cost_yftype, nnzobs,       &  call mdsreadfield( fname2, cost_iprec, cost_yftype, nnzobs,
176       &         localobs, localrec, mythid )       &         localobs, localrec, mythid )
177    
178            call ecco_zero(localdif,nnzbar,spzeroloc,myThid)
179    C mask the model field is observation if out-of-bound or missing.
180            call ecco_obsmsk( localbar, nnzbar, localobs, nnzobs,
181         &     localmask, spminloc, spmaxloc, spzeroloc,
182         &     localdif, difmask, myThid )
183    C Accumulate the model fields to compute time-mean at end of the loop.
184            call ecco_addmask(localdif,difmask, nnzbar,localbarmean,
185         &    localbarnum, nnzbar,myThid)
186    
187            call ecco_zero(localdif,nnzobs,spzeroloc,myThid)
188    C Accumulate the observation fields to compute time-mean at end of the loop.
189            call ecco_obsmsk( localobs, nnzobs, localobs, nnzobs,
190         &     localmask, spminloc, spmaxloc, spzeroloc,
191         &     localdif, difmask, myThid )
192            call ecco_addmask(localdif,difmask, nnzobs,localobsmean,
193         &    localobsnum, nnzobs,myThid)
194          enddo
195    C Compute the time-mean fields
196          call ecco_divfield(localbarmean,nnzbar,localbarnum,myThid)
197          call ecco_divfield(localobsmean,nnzobs,localobsnum,myThid)
198    
199    C Here only computes time-mean costs.
200          if ( preproc(1) .EQ. 'mean') then
201  c--     Compute masked model-data difference  c--     Compute masked model-data difference
202          call cost_gendif( localbar, nnzbar, localobs, nnzobs,  
203            call ecco_diffmsk( localbarmean, nnzbar, localobsmean, nnzobs,
204       &     localmask, spminloc, spmaxloc, spzeroloc,       &     localmask, spminloc, spmaxloc, spzeroloc,
205       &     posproc, scalefile, localdif, difmask, myThid )       &     localdif, difmask, myThid )
206    
207    #ifdef ALLOW_SMOOTH
208          if ( (useSMOOTH).AND.(posproc(1).EQ.'smooth').AND.
209         &     (nnzbar.EQ.1).AND.(nnzobs.EQ.1) )
210         &  call smooth_hetero2d(localdif,maskc,
211         &     posproc_c(1),posproc_i(1),mythid)
212    #endif
213    
214  c--     Compute normalized model-obs cost function  c--     Compute normalized model-obs cost function
215          do bj = jtlo,jthi          call ecco_addcost(
216            do bi = itlo,ithi       I                   localdif, localweight, difmask, nnzobs,
217              localcost    = 0. _d 0       I                   objf_local, num_local,
218              do k = 1,nnzobs       I                   myThid
219               do j = jmin,jmax       &                   )
220                do i = imin,imax  
221                  localwww  = localweight(i,j,k,bi,bj)  c--     output model-data difference to disk
222       &                    * difmask(i,j,k,bi,bj)          if ( outlev.GT.0 ) then
223                  junk      = localdif(i,j,k,bi,bj)            il=ilnblnk(outname)
224                  localcost = localcost + junk*junk*localwww            write(fname3(1:128),'(2a)') 'misfit_', outname(1:il)
225                  if ( localwww .ne. 0. )            if ( nnzobs.EQ.1 ) CALL
226       &               num_local(bi,bj) = num_local(bi,bj) + 1. _d 0       &    WRITE_REC_XY_RL( fname3, localdif,1, eccoiter, mythid )
227                enddo            if ( nnzobs.EQ.nr ) CALL
228               enddo       &    WRITE_REC_XYZ_RL( fname3, localdif,1, eccoiter, mythid )
229              enddo          endif
230              objf_local(bi,bj) = objf_local(bi,bj) + localcost         endif
231            enddo        endif
232          enddo  c
233    C The main time loop to compute costs
234    c-- reset nrecloop, if needed, according to preproc
235          if ( preproc(1) .EQ. 'clim') nrecloop=MIN(nrecloop,preproc_i(1))
236    
237          if ( .NOT. ( localobsfile.EQ.' ' ) ) then
238    
239          if ( preproc(1) .NE. 'mean') then
240    c--   loop over obsfile records
241          do irec = 1, nrecloop
242    
243    c--     determine records and file names
244            exst=.FALSE.
245            call cost_gencal(localbarfile, localobsfile,
246         &     irec, localstartdate, localperiod, fname1,
247         &     fname2, localrec, obsrec, exst, mythid )
248    
249    c--     load model average and observed average
250            call ecco_zero(localbar,nnzbar,zeroRL,myThid)
251            call cost_genread( fname1, localbar, localtmp, irec, nnzbar,
252         &       nrecloc, preproc, preproc_c, preproc_i, preproc_r,
253         &       dummy, mythid )
254    
255            call ecco_zero(localobs,Nr,spzeroloc,myThid)
256            if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) )
257         &  call mdsreadfield( fname2, cost_iprec, cost_yftype, nnzobs,
258         &         localobs, localrec, mythid )
259          if ( preproc(1) .EQ. 'anom' .OR. preproc(2) .EQ. 'anom') then
260    c--     Compute masked time-anomaly model-data difference
261            call ecco_diffanommsk( localbar, localbarmean, nnzbar,
262         &     localobs, localobsmean, nnzobs,
263         &     localmask, spminloc, spmaxloc, spzeroloc,
264         &     localdif, difmask, myThid )
265          else
266    c--     Compute masked (total) model-data difference
267            call ecco_diffmsk( localbar, nnzbar, localobs, nnzobs,
268         &     localmask, spminloc, spmaxloc, spzeroloc,
269         &     localdif, difmask, myThid )
270          endif
271    
272    #ifdef ALLOW_SMOOTH
273          if ( (useSMOOTH).AND.(posproc(1).EQ.'smooth').AND.
274         &     (nnzbar.EQ.1).AND.(nnzobs.EQ.1) )
275         &  call smooth_hetero2d(localdif,maskc,
276         &     posproc_c(1),posproc_i(1),mythid)
277    #endif
278    
279    c--     Compute normalized model-obs cost function
280            call ecco_addcost(
281         I                   localdif, localweight, difmask, nnzobs,
282         I                   objf_local, num_local,
283         I                   myThid
284         &                   )
285    
286  c--     output model-data difference to disk  c--     output model-data difference to disk
287          if ( outlev.GT.0 ) then          if ( outlev.GT.0 ) then
# Line 223  c--     output model-data difference to Line 296  c--     output model-data difference to
296        enddo        enddo
297  c--   End of loop over obsfile records.  c--   End of loop over obsfile records.
298    
299          endif /* if preproc(1) .NE. 'mean' */
300        endif        endif
301    
302  #endif /* ALLOW_ECCO */  #endif /* ALLOW_ECCO */
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.31
  ViewVC Help
Powered by ViewVC 1.1.22