/[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.14 by gforget, Tue Feb 28 00:51:07 2012 UTC revision 1.26 by gforget, Tue Oct 7 14:45:53 2014 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "COST_CPPOPTIONS.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, localstartdate, localperiod,       &     nrecloc, nrecobs,
14         &     localstartdate, localperiod,
15       &     ylocmask, localweight,       &     ylocmask, localweight,
16       &     spminloc, spmaxloc, spzeroloc,       &     spminloc, spmaxloc, spzeroloc,
17         &     preproc, preproc_c, preproc_i, preproc_r,
18         &     posproc, posproc_c, posproc_i, posproc_r,
19         &     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 35  c     == global variables == Line 37  c     == global variables ==
37  #ifdef ALLOW_CAL  #ifdef ALLOW_CAL
38  # include "cal.h"  # include "cal.h"
39  #endif  #endif
40  #ifdef ALLOW_COST  #ifdef ALLOW_ECCO
41  # include "ecco_cost.h"  # include "ecco.h"
42  # include "optim.h"  #endif
43  # ifdef ALLOW_SEAICE  #ifdef ALLOW_SEAICE
44  #  include "SEAICE_COST.h"  # include "SEAICE_COST.h"
 # endif  
45  #endif  #endif
46    
47  c     == routine arguments ==  c     == routine arguments ==
48    
       integer nnzbar  
       integer nnzobs  
       integer nrecloc  
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
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*(MAX_LEN_FNAM) preproc(NGENPPROC)
76          character*(MAX_LEN_FNAM) preproc_c(NGENPPROC)
77          character*(MAX_LEN_FNAM) posproc(NGENPPROC)
78          character*(MAX_LEN_FNAM) posproc_c(NGENPPROC)
79          character*(MAX_LEN_FNAM) outname
80    
81    #ifdef ALLOW_ECCO
82    
 #ifdef ALLOW_COST  
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        integer localrec, obsrec
       integer obsrec  
   
       logical doglobalread  
       logical ladinit  
   
       _RL spval  
       parameter (spval = -9999. )  
       _RL localwww  
       _RL localcost  
       _RL junk  
91    
92        _RL localmask  (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)        _RL localmask  (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
93        _RL localobs   (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)        _RL localobs   (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)
94        _RL cmask (1-olx:snx+olx,1-oly:sny+oly,nnzobs)        _RL localdif   (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)
95          _RL difmask    (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)
96    
97        character*(128) fname1, fname2        character*(128) fname1, fname2, fname3
       character*(MAX_LEN_MBUF) msgbuf  
98    
 cnew(  
       _RL daytime  
       _RL diffsecs  
       integer dayiter  
       integer daydate(4)  
       integer difftime(4)  
       integer middate(4)  
       integer yday, ymod  
       integer md, dd, sd, ld, wd  
       integer mody, modm  
       integer beginmodel, beginlocal  
99        logical exst        logical exst
 cnew)  
100    
101  c     == external functions ==  c     == external functions ==
102    
103        integer  ilnblnk        integer  ilnblnk
104        external ilnblnk        external ilnblnk
105    
106  c     == end of interface ==  CEOP
107    
108        jtlo = mybylo(mythid)        jtlo = mybylo(mythid)
109        jthi = mybyhi(mythid)        jthi = mybyhi(mythid)
110        itlo = mybxlo(mythid)        itlo = mybxlo(mythid)
111        ithi = mybxhi(mythid)        ithi = mybxhi(mythid)
       jmin = 1  
       jmax = sny  
       imin = 1  
       imax = snx  
112    
113  c--   Initialise local variables.  c--   Initialise local variables.
114    
       localwww = 0. _d 0  
   
115        do bj = jtlo,jthi        do bj = jtlo,jthi
116          do bi = itlo,ithi          do bi = itlo,ithi
117            objf_local(bi,bj) = 0. _d 0            objf_local(bi,bj) = 0. _d 0
118            num_local(bi,bj) = 0. _d 0            num_local(bi,bj) = 0. _d 0
           do k = 1,nnzobs  
             do j = jmin,jmax  
               do i = imin,imax  
                 localobs(i,j,k,bi,bj) = 0. _d 0  
               enddo  
             enddo  
           enddo  
119          enddo          enddo
120        enddo        enddo
121          call ecco_zero(localobs,nnzobs,zeroRL,myThid)
122          call ecco_zero(localdif,nnzobs,zeroRL,myThid)
123          call ecco_zero(difmask,nnzobs,zeroRL,myThid)
124    
125  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  
126        if ( ylocmask .EQ. 'C' .OR. ylocmask .EQ. 'c' ) then        if ( ylocmask .EQ. 'C' .OR. ylocmask .EQ. 'c' ) then
127           localmask(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)          call ecco_cprsrl(maskC,nr,localmask,nr,myThid)
128        elseif ( ylocmask .EQ. 'S' .OR. ylocmask .EQ. 's' ) then        elseif ( ylocmask .EQ. 'S' .OR. ylocmask .EQ. 's' ) then
129           localmask(i,j,k,bi,bj) = maskS(i,j,k,bi,bj)          call ecco_cprsrl(maskS,nr,localmask,nr,myThid)
130        elseif ( ylocmask .EQ. 'W' .OR. ylocmask .EQ. 'w' ) then        elseif ( ylocmask .EQ. 'W' .OR. ylocmask .EQ. 'w' ) then
131           localmask(i,j,k,bi,bj) = maskW(i,j,k,bi,bj)          call ecco_cprsrl(maskW,nr,localmask,nr,myThid)
132        else        else
133           STOP 'cost_generic: wrong ylocmask'           STOP 'cost_generic: wrong ylocmask'
134        endif        endif
               enddo  
             enddo  
           enddo  
         enddo  
       enddo  
135    
136  c--   First, read tiled data.  c-- reset nrecloc, if needed, according to preproc
137        doglobalread = .false.        if ( preproc(1) .EQ. 'climmon') nrecloc=MIN(nrecloc,12)
       ladinit      = .false.  
   
       write(fname1(1:128),'(80a)') ' '  
       il=ilnblnk( localbarfile )  
       write(fname1(1:128),'(2a,i10.10)')  
      &     localbarfile(1:il),'.',optimcycle  
138    
 cph      if ( .NOT. ( mult_local.EQ.0. .OR. localobsfile.EQ.' ' ) ) then  
139        if ( .NOT. ( localobsfile.EQ.' ' ) ) then        if ( .NOT. ( localobsfile.EQ.' ' ) ) then
140    
141  c--   Loop over records for the second time.  c--   loop over obsfile records
142        do irec = 1, nrecloc        do irec = 1, nrecloc
143    
144          if ( nnzbar .EQ. 1 ) then  c--     determine records and file names
145             call active_read_xy( fname1, localbar, irec, doglobalread,          exst=.FALSE.
146       &                      ladinit, optimcycle, mythid,          call cost_gencal(localbarfile, localobsfile,
147       &                      xx_localbar_mean_dummy )       &     irec, localstartdate, localperiod, fname1,
148          else       &     fname2, localrec, obsrec, exst, mythid )
149             call active_read_xyz( fname1, localbar, irec, doglobalread,  
150       &                       ladinit, optimcycle, mythid,  c--     load model average and observed average
151       &                       xx_localbar_mean_dummy )          call cost_genread( fname1, localbar, irec, nnzbar,
152          endif       &       preproc, dummy, mythid )
153    
154  cnew(          call ecco_zero(localobs,nnzobs,spzeroloc,myThid)
155          if ( localperiod .EQ. 86400. ) then          if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) )
156  c-- assume daily fields       &  call mdsreadfield( fname2, cost_iprec, cost_yftype, nnzobs,
            obsrec = irec  
            daytime = FLOAT(secondsperday*(irec-1)) + modelstart  
            dayiter = hoursperday*(irec-1) + modeliter0  
            call cal_getdate( dayiter, daytime, daydate, mythid )  
            call cal_convdate( daydate,yday,md,dd,sd,ld,wd,mythid )  
            ymod = localstartdate(1)/10000  
            if ( ymod .GE. yday ) then  
               call cal_FullDate( localstartdate(1), 0, middate, mythid)  
            else  
               middate(1) = yday*10000+100+1  
               call cal_FullDate( middate(1), 0, middate, mythid)  
            endif  
            call cal_TimePassed( middate, daydate, difftime, mythid )  
            call cal_ToSeconds( difftime, diffsecs, mythid )  
 c           localrec = floor(diffsecs/localperiod) + 1  
            localrec = int(diffsecs/localperiod) + 1  
         else  
 c-- assume monthly fields  
            beginlocal = localstartdate(1)/10000  
            beginmodel = modelstartdate(1)/10000  
            obsrec =  
      &           ( beginmodel - beginlocal )*nmonthyear  
      &         + ( mod(modelstartdate(1)/100,100)  
      &            -mod(localstartdate(1)/100,100) )  
      &         + irec  
            mody   = modelstartdate(1)/10000  
            modm   = modelstartdate(1)/100 - mody*100  
            yday   = mody + INT((modm-1+irec-1)/12)  
            localrec = 1 + MOD(modm-1+irec-1,12)  
         endif  
   
         il=ilnblnk(localobsfile)  
         write(fname2(1:128),'(2a,i4)')  
      &       localobsfile(1:il), '_', yday  
         inquire( file=fname2, exist=exst )  
         if ( (.NOT. exst).AND.( localperiod .NE. 86400. ) ) then  
            write(fname2(1:128),'(a)') localobsfile(1:il)  
            inquire( file=fname2, exist=exst )  
 #ifndef COST_GENERIC_ASSUME_CYCLIC  
 c assume we have one big file, one year after the other  
            localrec = obsrec  
 c otherwise assume climatology, used for each year  
 #endif  
         endif  
   
         if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) ) then  
           call mdsreadfield( fname2, cost_iprec, cost_yftype, nnzobs,  
157       &         localobs, localrec, mythid )       &         localobs, localrec, mythid )
         else  
           do bj = jtlo,jthi  
             do bi = itlo,ithi  
               do k = 1,nnzobs  
                 do j = jmin,jmax  
                   do i = imin,imax  
                      localobs(i,j,k,bi,bj) = spval  
                   enddo  
                 enddo  
               enddo  
             enddo  
           enddo  
         endif  
 cnew)  
   
         do bj = jtlo,jthi  
           do bi = itlo,ithi  
   
             localcost    = 0. _d 0  
   
 c--           Determine the mask on weights  
             do k = 1,nnzobs  
              do j = jmin,jmax  
               do i = imin,imax  
                cmask(i,j,k) = cosphi(i,j,bi,bj)*localmask(i,j,k,bi,bj)  
                 if ( localobs(i,j,k,bi,bj) .lt. spminloc .or.  
      &               localobs(i,j,k,bi,bj) .gt. spmaxloc .or.  
      &               localobs(i,j,k,bi,bj) .eq. spzeroloc ) then  
                    cmask(i,j,k) = 0. _d 0  
                 endif  
               enddo  
              enddo  
             enddo  
 c--  
             do k = 1,nnzobs  
              do j = jmin,jmax  
               do i = imin,imax  
                 localwww  = localweight(i,j,k,bi,bj)*cmask(i,j,k)  
                 junk      = ( localbar(i,j,k,bi,bj) -  
      &                        localobs(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  
158    
159              objf_local(bi,bj) = objf_local(bi,bj) + localcost  c--     Compute masked model-data difference
160            call ecco_diffmsk( localbar, nnzbar, localobs, nnzobs,
161            enddo       &     localmask, spminloc, spmaxloc, spzeroloc,
162          enddo       &     localdif, difmask, myThid )
163    
164    c--     Compute normalized model-obs cost function
165            call ecco_addcost(
166         I                   localdif, localweight, difmask, nnzobs,
167         I                   objf_local, num_local,
168         I                   myThid
169         &                   )
170    
171    c--     output model-data difference to disk
172            if ( outlev.GT.0 ) then
173              il=ilnblnk(outname)
174              write(fname3(1:128),'(2a)') 'misfit_', outname(1:il)
175              if ( nnzobs.EQ.1 ) CALL
176         &    WRITE_REC_XY_RL( fname3, localdif,irec, eccoiter, mythid )
177              if ( nnzobs.EQ.nr ) CALL
178         &    WRITE_REC_XYZ_RL( fname3, localdif,irec, eccoiter, mythid )
179            endif
180    
181        enddo        enddo
182  c--   End of second loop over records.  c--   End of loop over obsfile records.
183    
 c--   End of mult_local or localobsfile  
184        endif        endif
185    
186  #endif /* ifdef ALLOW_COST */  #endif /* ALLOW_ECCO */
187    
188        end        end
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.22