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

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

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

revision 1.1 by heimbach, Thu Nov 6 22:10:07 2003 UTC revision 1.2 by heimbach, Mon Oct 11 16:38:53 2004 UTC
# Line 7  Line 7 
7        subroutine cost_obcsw(        subroutine cost_obcsw(
8       I                       myiter,       I                       myiter,
9       I                       mytime,       I                       mytime,
      I                       startrec,  
      I                       endrec,  
10       I                       mythid       I                       mythid
11       &                     )       &                     )
12    
# Line 46  c     == routine arguments == Line 44  c     == routine arguments ==
44        integer myiter        integer myiter
45        _RL     mytime        _RL     mytime
46        integer mythid        integer mythid
       integer startrec  
       integer endrec  
47    
48  c     == local variables ==  c     == local variables ==
49    
# Line 61  c     == local variables == Line 57  c     == local variables ==
57        integer il        integer il
58        integer iobcs        integer iobcs
59        integer ip1        integer ip1
       integer nrec  
       integer ilfld  
       integer igg  
60    
61        _RL fctile        _RL fctile
62        _RL fcthread        _RL fcthread
63        _RL dummy        _RL dummy
       _RL gg  
       _RL tmpx  
 cgg(  
       _RL tmpfield (1-oly:sny+oly,nr,nsx,nsy)  
       _RL barofield (1-oly:sny+oly,nsx,nsy)  
       _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)  
       _RL utop  
64    
65        character*(80) fnamefld        character*(80) fnametheta
66          character*(80) fnamesalt
67          character*(80) fnameuvel
68          character*(80) fnamevvel
69    
70        logical doglobalread        logical doglobalread
71        logical ladinit        logical ladinit
# Line 105  c--   Read tiled data. Line 94  c--   Read tiled data.
94        doglobalread = .false.        doglobalread = .false.
95        ladinit      = .false.        ladinit      = .false.
96    
 c     Number of records to be used.  
       nrec = endrec-startrec+1  
   
97  #ifdef ALLOW_OBCSW_COST_CONTRIBUTION  #ifdef ALLOW_OBCSW_COST_CONTRIBUTION
98    
99        ip1 = 1        ip1 = 1
100        fcthread = 0. _d 0        fcthread = 0. _d 0
101    
102  #ifdef ECCO_VERBOSE  c--   Loop over records.
103        _BEGIN_MASTER( mythid )        do irec = 1,nmonsrec
       write(msgbuf,'(a)') ' '  
       call print_message( msgbuf, standardmessageunit,  
      &                    SQUEEZE_RIGHT , mythid)  
       write(msgbuf,'(a)') ' '  
       call print_message( msgbuf, standardmessageunit,  
      &                    SQUEEZE_RIGHT , mythid)  
       write(msgbuf,'(a,i9.8)')  
      &  ' cost_obcsw: number of records to process: ',nrec  
       call print_message( msgbuf, standardmessageunit,  
      &                    SQUEEZE_RIGHT , mythid)  
       write(msgbuf,'(a)') ' '  
       call print_message( msgbuf, standardmessageunit,  
      &                    SQUEEZE_RIGHT , mythid)  
       _END_MASTER( mythid )  
 #endif  
104    
105        if (optimcycle .ge. 0) then  c--     temperature
106          ilfld=ilnblnk( xx_obcsw_file )          iobcs = 1
107          write(fnamefld(1:80),'(2a,i10.10)')  c--     Read time averages and the monthly mean data.
108       &       xx_obcsw_file(1:ilfld), '.', optimcycle          il = ilnblnk( tbarfile )
109        endif          write(fnametheta(1:80),'(2a,i10.10)')
110         &       tbarfile(1:il),'.',optimcycle
111            call active_read_xyz( fnametheta, tbar, irec,
112         &                        doglobalread, ladinit,
113         &                        optimcycle, mythid,
114         &                        xx_tbar_mean_dummy )
115    
116  c--   Loop over records.          call mdsreadfieldyz( OBWtFile, readBinaryPrec, 'RS',
117        do irec = 1,nrec       &                       nr, OBWt, irec, mythid)
118    
119          do bj = jtlo,jthi          do bj = jtlo,jthi
120            do bi = itlo,ithi            do bi = itlo,ithi
121              do j = jmin,jmax  c--         Loop over the model layers
122                barofield(j,bi,bj) = 0. _d 0              fctile = 0. _d 0
123                do k = 1,nr
124    c--           Compute model data misfit and cost function term for
125    c             the temperature field.
126                   do j = jmin,jmax
127                      i = OB_Iw(J,bi,bj)
128                      if (maskS(i+ip1,j,k,bi,bj) .ne. 0.) then
129                         fctile = fctile +
130         &                    wobcsw(k,iobcs)*cosphi(i,j,bi,bj)*
131         &                    (tbar(i,j,k,bi,bj) - OBWt(j,k,bi,bj))*
132         &                    (tbar(i,j,k,bi,bj) - OBWt(j,k,bi,bj))
133                      endif
134                   enddo
135              enddo              enddo
136    c--         End of loop over layers.
137                fcthread         = fcthread           + fctile
138                objf_obcsw(bi,bj) = objf_obcsw(bi,bj) + fctile
139            enddo            enddo
140          enddo          enddo
141    
142          call active_read_yz( fnamefld, tmpfield, irec, doglobalread,  c--     salt
143       &                       ladinit, optimcycle, mythid          iobcs = 2
144       &        , xx_obcsw_dummy )  c--     Read time averages and the monthly mean data.
145            il = ilnblnk( sbarfile )
146  cgg    Need to solve for iobcs would have been.          write(fnamesalt(1:80),'(2a,i10.10)')
147            gg    = (irec-1)/nobcs       &       sbarfile(1:il),'.',optimcycle
148            igg   = int(gg)          call active_read_xyz( fnamesalt, sbar, irec,
149            iobcs = irec - igg*nobcs       &                        doglobalread, ladinit,
150         &                        optimcycle, mythid,
151            call active_read_yz( 'maskobcsw', maskyz,       &                        xx_sbar_mean_dummy )
152       &                         iobcs,  
153       &                         doglobalread, ladinit, 0,          call mdsreadfieldyz( OBWsFile, readBinaryPrec, 'RS',
154       &                         mythid, dummy )       &                       nr, OBWs, irec, mythid)
155    
156  cgg     Need to change xx field to baroclinic vel.          do bj = jtlo,jthi
157  cgg     Level Nr contains barotropic vel, remove it.            do bi = itlo,ithi
158  cgg     The deepest wet level velocity should be  c--         Loop over the model layers
159  cgg     computed in order to ensure zero integrated flux.              fctile = 0. _d 0
160          if (iobcs .eq. 3 .or. iobcs .eq. 4) then              do k = 1,nr
161  c--     Loop over this thread's tiles.  c--           Compute model data misfit and cost function term for
162            do bj = jtlo,jthi  c             the temperature field.
163              do bi = itlo,ithi                 do j = jmin,jmax
164                do j = jmin,jmax                    i = OB_Iw(J,bi,bj)
165                      if (maskS(i+ip1,j,k,bi,bj) .ne. 0.) then
166  cgg         The barotropic velocity is stored in level 1.                       fctile = fctile +
167  cgg         Penalize for deviation in it.       &                    wobcsw(k,iobcs)*cosphi(i,j,bi,bj)*
168                  barofield(j,bi,bj) = tmpfield(j,1,bi,bj)       &                    (sbar(i,j,k,bi,bj) - OBWs(j,k,bi,bj))*
169       &                             * maskyz(j,1,bi,bj)       &                    (sbar(i,j,k,bi,bj) - OBWs(j,k,bi,bj))
170                      endif
171                  tmpfield(j,1,bi,bj) = 0.d0                 enddo
                 utop = 0.d0  
   
                 do k = 1,Nr  
 cgg    If cells are not full, this should be modified with hFac.  
 cgg      
 cgg    The xx field (tmpfldxz) does not contain the velocity at the  
 cgg    lowest wet level. This velocity is not independent; it must  
 cgg    exactly balance the volume flux, since we are dealing with  
 cgg    the baroclinic velocity structure..  
                   utop = tmpfield(j,k,bi,bj)*  
      &               maskyz(j,k,bi,bj) * delZ(k) + utop  
                 enddo  
                 tmpfield(j,1,bi,bj) = -utop / delZ(1)  
               enddo  
172              enddo              enddo
173    c--         End of loop over layers.
174                fcthread         = fcthread           + fctile
175                objf_obcsw(bi,bj) = objf_obcsw(bi,bj) + fctile
176            enddo            enddo
177            enddo
178    
179          endif  c--     uvel
180            iobcs = 3
181    c--     Read time averages and the monthly mean data.
182            il = ilnblnk( ubarfile )
183            write(fnameuvel(1:80),'(2a,i10.10)')
184         &       ubarfile(1:il),'.',optimcycle
185            call active_read_xyz( fnameuvel, ubar, irec,
186         &                        doglobalread, ladinit,
187         &                        optimcycle, mythid,
188         &                        dummy )
189    
190  c--     Loop over this thread's tiles.          call mdsreadfieldyz( OBWuFile, readBinaryPrec, 'RS',
191         &                       nr, OBWu, irec, mythid)
192          do bj = jtlo,jthi          do bj = jtlo,jthi
193            do bi = itlo,ithi            do bi = itlo,ithi
194    c--         Loop over the model layers
 c--         Determine the weights to be used.  
195              fctile = 0. _d 0              fctile = 0. _d 0
196                do k = 1,nr
197              if (iobcs .eq. 3 .or. iobcs .eq. 4) then  c--           Compute model data misfit and cost function term for
198                do j = jmin,jmax  c             the temperature field.
199                  i = OB_Iw(j,bi,bj)                 do j = jmin,jmax
200                  tmpx = barofield(j,bi,bj)                    i = OB_Iw(J,bi,bj)
201                    fctile = fctile                    if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
202       &               + wbaro*cosphi(i,j,bi,bj)                       fctile = fctile +
203       &                 *tmpx*tmpx       &                    wobcsw(k,iobcs)*cosphi(i,j,bi,bj)*
204       &                 *maskyz(j,1,bi,bj)       &                    (ubar(i,j,k,bi,bj) - OBWu(j,k,bi,bj))*
205  cgg                print*,'fctile barotropic W',fctile       &                    (ubar(i,j,k,bi,bj) - OBWu(j,k,bi,bj))
206                enddo                    endif
207              endif                 enddo
   
             do k = 1, Nr  
               do j = jmin,jmax  
                 i = OB_Iw(j,bi,bj)  
 cgg                if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then  
                   tmpx = tmpfield(j,k,bi,bj)  
                   fctile = fctile  
      &                 + wobcswLev(j,k,bi,bj,iobcs)*cosphi(i,j,bi,bj)  
      &                 *tmpx*tmpx  
      &                 *maskyz(j,k,bi,bj)  
 cgg                endif  
               enddo  
208              enddo              enddo
209    c--         End of loop over layers.
210                fcthread         = fcthread           + fctile
211                objf_obcsw(bi,bj) = objf_obcsw(bi,bj) + fctile
212              enddo
213            enddo
214    
215    c--     vvel
216            iobcs = 4
217    c--     Read time averages and the monthly mean data.
218            il = ilnblnk( vbarfile )
219            write(fnamevvel(1:80),'(2a,i10.10)')
220         &       vbarfile(1:il),'.',optimcycle
221            call active_read_xyz( fnamevvel, vbar, irec,
222         &                        doglobalread, ladinit,
223         &                        optimcycle, mythid,
224         &                        dummy )
225    
226            call mdsreadfieldyz( OBWvFile, readBinaryPrec, 'RS',
227         &                       nr, OBWv, irec, mythid)
228    
229            do bj = jtlo,jthi
230              do bi = itlo,ithi
231    c--         Loop over the model layers
232                fctile = 0. _d 0
233                do k = 1,nr
234    c--           Compute model data misfit and cost function term for
235    c             the temperature field.
236                   do j = jmin,jmax
237                      i = OB_Iw(J,bi,bj)
238                      if (maskS(i,j,k,bi,bj) .ne. 0.) then
239                         fctile = fctile +
240         &                    wobcsw(k,iobcs)*cosphi(i,j,bi,bj)*
241         &                    (vbar(i,j,k,bi,bj) - OBWv(j,k,bi,bj))*
242         &                    (vbar(i,j,k,bi,bj) - OBWv(j,k,bi,bj))
243                      endif
244                   enddo
245                enddo
246    c--         End of loop over layers.
247                fcthread         = fcthread           + fctile
248              objf_obcsw(bi,bj) = objf_obcsw(bi,bj) + fctile              objf_obcsw(bi,bj) = objf_obcsw(bi,bj) + fctile
             fcthread         = fcthread + fctile  
249            enddo            enddo
250          enddo          enddo
251    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22