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

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

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

revision 1.4 by jmc, Tue Apr 28 18:13:28 2009 UTC revision 1.5 by mlosch, Wed Jan 19 13:46:19 2011 UTC
# Line 6  C $Name$ Line 6  C $Name$
6  # include "OBCS_OPTIONS.h"  # include "OBCS_OPTIONS.h"
7  #endif  #endif
8    
9    CBOP
10    C     !ROUTINE: COST_OBCSS
11    C     !INTERFACE:
12        subroutine cost_obcss(        subroutine cost_obcss(
13       I                       myiter,       I                       myiter,
14       I                       mytime,       I                       mytime,
15         I                       startrec,
16         I                       endrec,
17       I                       mythid       I                       mythid
18       &                     )       &                     )
19    
20    C     !DESCRIPTION: \bv
21  c     ==================================================================  c     ==================================================================
22  c     SUBROUTINE cost_obcss  c     SUBROUTINE cost_obcss
23  c     ==================================================================  c     ==================================================================
# Line 21  c Line 27  c
27  c     ==================================================================  c     ==================================================================
28  c     SUBROUTINE cost_obcss  c     SUBROUTINE cost_obcss
29  c     ==================================================================  c     ==================================================================
30    C     \ev
31    
32    C     !USES:
33    
34        implicit none        implicit none
35    
# Line 41  c     == global variables == Line 50  c     == global variables ==
50  #include "ctrl_dummy.h"  #include "ctrl_dummy.h"
51  #include "optim.h"  #include "optim.h"
52    
53    C     !INPUT/OUTPUT PARAMETERS:
54  c     == routine arguments ==  c     == routine arguments ==
55    
56        integer myiter        integer myiter
57        _RL     mytime        _RL     mytime
58        integer mythid        integer mythid
59          integer startrec
60          integer endrec
61    
62    C     !LOCAL VARIABLES:
63  c     == local variables ==  c     == local variables ==
64    
65        integer bi,bj        integer bi,bj
# Line 59  c     == local variables == Line 72  c     == local variables ==
72        integer il        integer il
73        integer iobcs        integer iobcs
74        integer jp1        integer jp1
75          integer nrec
76          integer ilfld
77          integer igg
78    
79        _RL fctile        _RL fctile
80        _RL fcthread        _RL fcthread
81        _RL dummy        _RL dummy
82          _RL gg
83          _RL tmpx
84          _RL tmpfield (1-olx:snx+olx,nr,nsx,nsy)
85          _RL maskxz   (1-olx:snx+olx,nr,nsx,nsy)
86    
87        character*(80) fnametheta        character*(80) fnamefld
       character*(80) fnamesalt  
       character*(80) fnameuvel  
       character*(80) fnamevvel  
88    
89        logical doglobalread        logical doglobalread
90        logical ladinit        logical ladinit
# Line 82  c     == external functions == Line 99  c     == external functions ==
99        external ilnblnk        external ilnblnk
100    
101  c     == end of interface ==  c     == end of interface ==
102    CEOP
103    
104        jtlo = mybylo(mythid)        jtlo = mybylo(mythid)
105        jthi = mybyhi(mythid)        jthi = mybyhi(mythid)
# Line 96  c--   Read tiled data. Line 114  c--   Read tiled data.
114        doglobalread = .false.        doglobalread = .false.
115        ladinit      = .false.        ladinit      = .false.
116    
117    c     Number of records to be used.
118          nrec = endrec-startrec+1
119    
120  #ifdef ALLOW_OBCSS_COST_CONTRIBUTION  #ifdef ALLOW_OBCSS_COST_CONTRIBUTION
121    
122        jp1 = 1        jp1 = 1
123        fcthread = 0. _d 0        fcthread = 0. _d 0
124    
125  c--   Loop over records.  #ifdef ECCO_VERBOSE
126        do irec = 1,nmonsrec        _BEGIN_MASTER( mythid )
127          write(msgbuf,'(a)') ' '
128  c--     temperature        call print_message( msgbuf, standardmessageunit,
129          iobcs = 1       &                    SQUEEZE_RIGHT , mythid)
130  c--     Read time averages and the monthly mean data.        write(msgbuf,'(a)') ' '
131          il = ilnblnk( tbarfile )        call print_message( msgbuf, standardmessageunit,
132          write(fnametheta(1:80),'(2a,i10.10)')       &                    SQUEEZE_RIGHT , mythid)
133       &       tbarfile(1:il),'.',optimcycle        write(msgbuf,'(a,i9.8)')
134          call active_read_xyz( fnametheta, tbar, irec,       &  ' cost_obcss: number of records to process: ',nrec
135       &                        doglobalread, ladinit,        call print_message( msgbuf, standardmessageunit,
136       &                        optimcycle, mythid,       &                    SQUEEZE_RIGHT , mythid)
137       &                        xx_tbar_mean_dummy )        write(msgbuf,'(a)') ' '
138          call print_message( msgbuf, standardmessageunit,
139         &                    SQUEEZE_RIGHT , mythid)
140          _END_MASTER( mythid )
141    #endif
142    
143          call mdsreadfieldxz( OBStFile, readBinaryPrec, 'RS',        if (optimcycle .ge. 0) then
144       &                       nr, OBSt, irec, mythid)          ilfld=ilnblnk( xx_obcss_file )
145            write(fnamefld(1:80),'(2a,i10.10)')
146         &       xx_obcss_file(1:ilfld), '.', optimcycle
147          endif
148    
149          do bj = jtlo,jthi  c--   Loop over records.
150            do bi = itlo,ithi        do irec = 1,nrec
 c--         Loop over the model layers  
             fctile = 0. _d 0  
             do k = 1,nr  
 c--           Compute model data misfit and cost function term for  
 c             the temperature field.  
                do i = imin,imax  
                   j = OB_Js(I,bi,bj)  
                   if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then  
                      fctile = fctile +  
      &                    wobcss(k,iobcs)*cosphi(i,j,bi,bj)*  
      &                    (tbar(i,j,k,bi,bj) - OBSt(i,k,bi,bj))*  
      &                    (tbar(i,j,k,bi,bj) - OBSt(i,k,bi,bj))  
                   endif  
                enddo  
             enddo  
 c--         End of loop over layers.  
             fcthread         = fcthread           + fctile  
             objf_obcss(bi,bj) = objf_obcss(bi,bj) + fctile  
           enddo  
         enddo  
151    
152  c--     salt          call active_read_xz( fnamefld, tmpfield, irec, doglobalread,
153          iobcs = 2       &                       ladinit, optimcycle, mythid
154  c--     Read time averages and the monthly mean data.       &        , xx_obcss_dummy )
         il = ilnblnk( sbarfile )  
         write(fnamesalt(1:80),'(2a,i10.10)')  
      &       sbarfile(1:il),'.',optimcycle  
         call active_read_xyz( fnamesalt, sbar, irec,  
      &                        doglobalread, ladinit,  
      &                        optimcycle, mythid,  
      &                        xx_sbar_mean_dummy )  
155    
156          call mdsreadfieldxz( OBSsFile, readBinaryPrec, 'RS',  cgg    Need to solve for iobcs would have been.
157       &                       nr, OBSs, irec, mythid)            gg    = (irec-1)/nobcs
158              igg   = int(gg)
159              iobcs = irec - igg*nobcs
160    cgg          print*,'S IREC, IOBCS',irec, iobcs
161    
162          do bj = jtlo,jthi            call active_read_xz( 'maskobcss', maskxz,
163            do bi = itlo,ithi       &                         iobcs,
164  c--         Loop over the model layers       &                         doglobalread, ladinit, 0,
165              fctile = 0. _d 0       &                         mythid, dummy )
             do k = 1,nr  
 c--           Compute model data misfit and cost function term for  
 c             the temperature field.  
                do i = imin,imax  
                   j = OB_Js(I,bi,bj)  
                   if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then  
                      fctile = fctile +  
      &                    wobcss(k,iobcs)*cosphi(i,j,bi,bj)*  
      &                    (sbar(i,j,k,bi,bj) - OBSs(i,k,bi,bj))*  
      &                    (sbar(i,j,k,bi,bj) - OBSs(i,k,bi,bj))  
                   endif  
                enddo  
             enddo  
 c--         End of loop over layers.  
             fcthread         = fcthread           + fctile  
             objf_obcss(bi,bj) = objf_obcss(bi,bj) + fctile  
           enddo  
         enddo  
166    
 c--     uvel  
         iobcs = 3  
 c--     Read time averages and the monthly mean data.  
         il = ilnblnk( ubarfile )  
         write(fnameuvel(1:80),'(2a,i10.10)')  
      &       ubarfile(1:il),'.',optimcycle  
         call active_read_xyz( fnameuvel, ubar, irec,  
      &                        doglobalread, ladinit,  
      &                        optimcycle, mythid,  
      &                        dummy )  
167    
168          call mdsreadfieldxz( OBSuFile, readBinaryPrec, 'RS',  c--     Loop over this thread's tiles.
      &                       nr, OBSu, irec, mythid)  
169          do bj = jtlo,jthi          do bj = jtlo,jthi
170            do bi = itlo,ithi            do bi = itlo,ithi
 c--         Loop over the model layers  
             fctile = 0. _d 0  
             do k = 1,nr  
 c--           Compute model data misfit and cost function term for  
 c             the temperature field.  
                do i = imin,imax  
                   j = OB_Js(I,bi,bj)  
                   if (maskW(i,j,k,bi,bj) .ne. 0.) then  
                      fctile = fctile +  
      &                    wobcss(k,iobcs)*cosphi(i,j,bi,bj)*  
      &                    (ubar(i,j,k,bi,bj) - OBSu(i,k,bi,bj))*  
      &                    (ubar(i,j,k,bi,bj) - OBSu(i,k,bi,bj))  
                   endif  
                enddo  
             enddo  
 c--         End of loop over layers.  
             fcthread         = fcthread           + fctile  
             objf_obcss(bi,bj) = objf_obcss(bi,bj) + fctile  
           enddo  
         enddo  
   
 c--     vvel  
         iobcs = 4  
 c--     Read time averages and the monthly mean data.  
         il = ilnblnk( vbarfile )  
         write(fnamevvel(1:80),'(2a,i10.10)')  
      &       vbarfile(1:il),'.',optimcycle  
         call active_read_xyz( fnamevvel, vbar, irec,  
      &                        doglobalread, ladinit,  
      &                        optimcycle, mythid,  
      &                        dummy )  
   
         call mdsreadfieldxz( OBSvFile, readBinaryPrec, 'RS',  
      &                       nr, OBSv, irec, mythid)  
171    
172          do bj = jtlo,jthi  c--         Determine the weights to be used.
           do bi = itlo,ithi  
 c--         Loop over the model layers  
173              fctile = 0. _d 0              fctile = 0. _d 0
174              do k = 1,nr  
175  c--           Compute model data misfit and cost function term for              do k = 1, Nr
176  c             the temperature field.                do i = imin,imax
177                 do i = imin,imax                   j = OB_Js(I,bi,bj)
178                    j = OB_Js(I,bi,bj)  cgg                 if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
179                    if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then                    tmpx = tmpfield(i,k,bi,bj)
180                       fctile = fctile +  CMM                  fctile = fctile + wobcss2(i,k,bi,bj,iobcs)
181       &                    wobcss(k,iobcs)*cosphi(i,j,bi,bj)*                    fctile = fctile + wobcss(k,iobcs)
182       &                    (vbar(i,j,k,bi,bj) - OBSv(i,k,bi,bj))*       &                        *tmpx*tmpx*maskxz(i,k,bi,bj)
183       &                    (vbar(i,j,k,bi,bj) - OBSv(i,k,bi,bj))  cgg                endif
184                    endif  CMM                  if (wobcss2(i,k,bi,bj,iobcs)*maskxz(i,k,bi,bj).ne.0.)
185                 enddo                    if (wobcss(k,iobcs)*maskxz(i,k,bi,bj).ne.0.)
186         &                    num_obcss(bi,bj) = num_obcss(bi,bj) + 1. _d 0
187    cgg                print*,'S fctile',fctile
188                  enddo
189              enddo              enddo
190  c--         End of loop over layers.  
             fcthread         = fcthread           + fctile  
191              objf_obcss(bi,bj) = objf_obcss(bi,bj) + fctile              objf_obcss(bi,bj) = objf_obcss(bi,bj) + fctile
192                fcthread         = fcthread + fctile
193            enddo            enddo
194          enddo          enddo
195    
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22