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

Diff of /MITgcm/pkg/ecco/cost_obcse.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.3 by jmc, Tue Oct 9 00:02:50 2007 UTC
# Line 1  Line 1 
1    C $Header$
2    C $Name$
3    
4  #include "COST_CPPOPTIONS.h"  #include "COST_CPPOPTIONS.h"
5  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
# Line 7  Line 9 
9        subroutine cost_obcse(        subroutine cost_obcse(
10       I                       myiter,       I                       myiter,
11       I                       mytime,       I                       mytime,
      I                       startrec,  
      I                       endrec,  
12       I                       mythid       I                       mythid
13       &                     )       &                     )
14    
# Line 18  c     ================================== Line 18  c     ==================================
18  c  c
19  c     o cost function contribution obc  c     o cost function contribution obc
20  c  c
 c     o G. Gebbie, gebbie@mit.edu, 18-Mar-2003  
21  c     ==================================================================  c     ==================================================================
22  c     SUBROUTINE cost_obcse  c     SUBROUTINE cost_obcse
23  c     ==================================================================  c     ==================================================================
# Line 47  c     == routine arguments == Line 46  c     == routine arguments ==
46        integer myiter        integer myiter
47        _RL     mytime        _RL     mytime
48        integer mythid        integer mythid
 cgg(    
       integer startrec  
       integer endrec  
 cgg)  
49    
50  c     == local variables ==  c     == local variables ==
51    
# Line 64  c     == local variables == Line 59  c     == local variables ==
59        integer il        integer il
60        integer iobcs        integer iobcs
61        integer ip1        integer ip1
       integer nrec  
       integer ilfld  
       integer igg  
62    
63        _RL fctile        _RL fctile
64        _RL fcthread        _RL fcthread
65        _RL dummy        _RL dummy
       _RL gg  
       _RL tmpx  
       _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 ubottom  
66    
67        character*(80) fnamefld        character*(80) fnametheta
68          character*(80) fnamesalt
69          character*(80) fnameuvel
70          character*(80) fnamevvel
71    
72        logical doglobalread        logical doglobalread
73        logical ladinit        logical ladinit
# Line 107  c--   Read tiled data. Line 96  c--   Read tiled data.
96        doglobalread = .false.        doglobalread = .false.
97        ladinit      = .false.        ladinit      = .false.
98    
 c     Number of records to be used.  
       nrec = endrec-startrec+1  
   
99  #ifdef ALLOW_OBCSE_COST_CONTRIBUTION  #ifdef ALLOW_OBCSE_COST_CONTRIBUTION
100    
101        ip1 = 0        ip1 = 0
102        fcthread = 0. _d 0        fcthread = 0. _d 0
103    
104  #ifdef ECCO_VERBOSE  c--   Loop over records.
105        _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_obcse: 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  
106    
107        if (optimcycle .ge. 0) then  c--     temperature
108          ilfld=ilnblnk( xx_obcse_file )          iobcs = 1
109          write(fnamefld(1:80),'(2a,i10.10)')  c--     Read time averages and the monthly mean data.
110       &       xx_obcse_file(1:ilfld), '.', optimcycle          il = ilnblnk( tbarfile )
111        endif          write(fnametheta(1:80),'(2a,i10.10)')
112         &       tbarfile(1:il),'.',optimcycle
113            call active_read_xyz( fnametheta, tbar, irec,
114         &                        doglobalread, ladinit,
115         &                        optimcycle, mythid,
116         &                        xx_tbar_mean_dummy )
117    
118  c--   Loop over records.          call mdsreadfieldyz( OBEtFile, ReadBinaryPrec, 'RS',
119        do irec = 1,nrec       &                       nr, OBEt, irec, mythid)
120    
121          do bj = jtlo,jthi          do bj = jtlo,jthi
122            do bi = itlo,ithi            do bi = itlo,ithi
123              do j = jmin,jmax  c--         Loop over the model layers
124                barofield(j,bi,bj) = 0. _d 0              fctile = 0. _d 0
125                do k = 1,nr
126    c--           Compute model data misfit and cost function term for
127    c             the temperature field.
128                   do j = jmin,jmax
129                      i = OB_Ie(J,bi,bj)
130                      if (maskS(i+ip1,j,k,bi,bj) .ne. 0.) then
131                         fctile = fctile +
132         &                    wobcse(k,iobcs)*cosphi(i,j,bi,bj)*
133         &                    (tbar(i,j,k,bi,bj) - OBEt(j,k,bi,bj))*
134         &                    (tbar(i,j,k,bi,bj) - OBEt(j,k,bi,bj))
135                      endif
136                   enddo
137              enddo              enddo
138    c--         End of loop over layers.
139                fcthread         = fcthread           + fctile
140                objf_obcse(bi,bj) = objf_obcse(bi,bj) + fctile
141            enddo            enddo
142          enddo          enddo
                 
         call active_read_yz( fnamefld, tmpfield, irec, doglobalread,  
      &                       ladinit, optimcycle, mythid  
      &        , xx_obcse_dummy )  
   
 cgg    Need to solve for iobcs would have been.  
           gg    = (irec-1)/nobcs  
           igg   = int(gg)  
           iobcs = irec - igg*nobcs  
143    
144  c--     Loop over this thread's tiles.  c--     salt
145            iobcs = 2
146    c--     Read time averages and the monthly mean data.
147            il = ilnblnk( sbarfile )
148            write(fnamesalt(1:80),'(2a,i10.10)')
149         &       sbarfile(1:il),'.',optimcycle
150            call active_read_xyz( fnamesalt, sbar, irec,
151         &                        doglobalread, ladinit,
152         &                        optimcycle, mythid,
153         &                        xx_sbar_mean_dummy )
154    
155            call mdsreadfieldyz( OBEsFile, readBinaryPrec, 'RS',
156         &                       nr, OBEs, irec, mythid)
157    
158          do bj = jtlo,jthi          do bj = jtlo,jthi
159            do bi = itlo,ithi            do bi = itlo,ithi
160    c--         Loop over the model layers
161            call active_read_yz( 'maskobcse', maskyz,              fctile = 0. _d 0
162       &                         iobcs,              do k = 1,nr
163       &                         doglobalread, ladinit, 0,  c--           Compute model data misfit and cost function term for
164       &                         mythid, dummy )  c             the temperature field.
165                   do j = jmin,jmax
166  cgg     Need to change xx field to baroclinic vel.                    i = OB_Ie(J,bi,bj)
167  cgg     Level Nr contains barotropic vel, remove it.                    if (maskS(i+ip1,j,k,bi,bj) .ne. 0.) then
168  cgg     The deepest wet level velocity should be                       fctile = fctile +
169  cgg     computed in order to ensure zero integrated flux.       &                    wobcse(k,iobcs)*cosphi(i,j,bi,bj)*
170          if (iobcs .eq. 3 .or. iobcs .eq. 4) then       &                    (sbar(i,j,k,bi,bj) - OBEs(j,k,bi,bj))*
171  c--     Loop over this thread's tiles.       &                    (sbar(i,j,k,bi,bj) - OBEs(j,k,bi,bj))
172            do bj = jtlo,jthi                    endif
173              do bi = itlo,ithi                 enddo
               do j = jmin,jmax  
   
 cgg         The barotropic velocity is stored in level 1.  
 cgg         Penalize for deviation in it.  
                 barofield(j,bi,bj) = tmpfield(j,1,bi,bj)  
      &                             * maskyz(j,1,bi,bj)  
   
                 tmpfield(j,1,bi,bj) = 0.d0  
                 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  
174              enddo              enddo
175    c--         End of loop over layers.
176                fcthread         = fcthread           + fctile
177                objf_obcse(bi,bj) = objf_obcse(bi,bj) + fctile
178            enddo            enddo
179          endif          enddo
180    
181  c--         Determine the weights to be used.  c--     uvel
182              fctile = 0. _d 0          iobcs = 3
183    c--     Read time averages and the monthly mean data.
184            il = ilnblnk( ubarfile )
185            write(fnameuvel(1:80),'(2a,i10.10)')
186         &       ubarfile(1:il),'.',optimcycle
187            call active_read_xyz( fnameuvel, ubar, irec,
188         &                        doglobalread, ladinit,
189         &                        optimcycle, mythid,
190         &                        dummy )
191    
192              if (iobcs .eq. 3 .or. iobcs .eq. 4) then          call mdsreadfieldyz( OBEuFile, readBinaryPrec, 'RS',
193                do j = jmin,jmax       &                       nr, OBEu, irec, mythid)
194                  i = OB_Ie(j,bi,bj)          do bj = jtlo,jthi
195                  tmpx = barofield(j,bi,bj)            do bi = itlo,ithi
196                    fctile = fctile  c--         Loop over the model layers
197       &               + wbaro*cosphi(i,j,bi,bj)              fctile = 0. _d 0
198       &                 *tmpx*tmpx              do k = 1,nr
199       &                 *maskyz(j,1,bi,bj)  c--           Compute model data misfit and cost function term for
200                  print*,'fctile barotropic E',fctile  c             the temperature field.
201                enddo                 do j = jmin,jmax
202              endif                    i = OB_Ie(J,bi,bj)
203                      if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
204              do k = 1, Nr                       fctile = fctile +
205                do j = jmin,jmax       &                    wobcse(k,iobcs)*cosphi(i,j,bi,bj)*
206                  i = OB_Ie(j,bi,bj)       &                    (ubar(i,j,k,bi,bj) - OBEu(j,k,bi,bj))*
207  cgg                if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then       &                    (ubar(i,j,k,bi,bj) - OBEu(j,k,bi,bj))
208                    tmpx = tmpfield(j,k,bi,bj)                    endif
209                    fctile = fctile                 enddo
      &                 + wobcseLev(j,k,bi,bj,iobcs)*cosphi(i,j,bi,bj)  
      &                 *tmpx*tmpx  
      &                 *maskyz(j,k,bi,bj)  
 cgg                endif  
               enddo  
210              enddo              enddo
211    c--         End of loop over layers.
212                fcthread         = fcthread           + fctile
213                objf_obcse(bi,bj) = objf_obcse(bi,bj) + fctile
214              enddo
215            enddo
216    
217    c--     vvel
218            iobcs = 4
219    c--     Read time averages and the monthly mean data.
220            il = ilnblnk( vbarfile )
221            write(fnamevvel(1:80),'(2a,i10.10)')
222         &       vbarfile(1:il),'.',optimcycle
223            call active_read_xyz( fnamevvel, vbar, irec,
224         &                        doglobalread, ladinit,
225         &                        optimcycle, mythid,
226         &                        dummy )
227    
228            call mdsreadfieldyz( OBEvFile, readBinaryPrec, 'RS',
229         &                       nr, OBEv, irec, mythid)
230    
231            do bj = jtlo,jthi
232              do bi = itlo,ithi
233    c--         Loop over the model layers
234                fctile = 0. _d 0
235                do k = 1,nr
236    c--           Compute model data misfit and cost function term for
237    c             the temperature field.
238                   do j = jmin,jmax
239                      i = OB_Ie(J,bi,bj)
240                      if (maskS(i,j,k,bi,bj) .ne. 0.) then
241                         fctile = fctile +
242         &                    wobcse(k,iobcs)*cosphi(i,j,bi,bj)*
243         &                    (vbar(i,j,k,bi,bj) - OBEv(j,k,bi,bj))*
244         &                    (vbar(i,j,k,bi,bj) - OBEv(j,k,bi,bj))
245                      endif
246                   enddo
247                enddo
248    c--         End of loop over layers.
249                fcthread         = fcthread           + fctile
250              objf_obcse(bi,bj) = objf_obcse(bi,bj) + fctile              objf_obcse(bi,bj) = objf_obcse(bi,bj) + fctile
             fcthread         = fcthread + fctile  
251            enddo            enddo
252          enddo          enddo
253    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22