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

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

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

revision 1.2 by heimbach, Mon Oct 11 16:38:53 2004 UTC revision 1.6 by heimbach, Sat May 3 00:08:42 2008 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "COST_CPPOPTIONS.h"  #include "COST_CPPOPTIONS.h"
5    
# Line 59  c     == local variables == Line 60  c     == local variables ==
60        integer irec        integer irec
61        integer ilw        integer ilw
62        integer nf, nl        integer nf, nl
63          integer minrec
64    
65        _RL fctilew        _RL fctilew
66        _RL fcthread_wdrift        _RL fcthread_wdrift
# Line 97  c--   Read tiled data. Line 99  c--   Read tiled data.
99          write(fnamew(1:80),'(2a,i10.10)')          write(fnamew(1:80),'(2a,i10.10)')
100       &    wbarfile(1:ilw),'.',optimcycle       &    wbarfile(1:ilw),'.',optimcycle
101        endif        endif
         
       if (nmonsrec.lt.12) then  
         print*  
         print*,' cost_Driftw: assimilation period smaller'  
         print*,'             than 1 year, no drift computed.'  
         print*  
         stop   '  ... stopped in cost_Driftw.'  
       endif  
102    
103        fcthread_wdrift = 0. _d 0        fcthread_wdrift = 0. _d 0
104          
105        do bj = jtlo,jthi        do bj = jtlo,jthi
106          do bi = itlo,ithi          do bi = itlo,ithi
107            do k = 1,nr                do k = 1,nr
108              do j = jmin,jmax              do j = jmin,jmax
109                do i = imin,imax                do i = imin,imax
110                  wfmean(i,j,k,bi,bj) = 0.0                  wfmean(i,j,k,bi,bj) = 0.0
111                  wlmean(i,j,k,bi,bj) = 0.0                  wlmean(i,j,k,bi,bj) = 0.0
112                enddo                enddo
113              enddo                              enddo
114            enddo            enddo
115          enddo          enddo
116        enddo        enddo
117          
118        nf = 0        nf = 0
119        nl = 0        nl = 0
120    c--   Number of full years
121          nfmin = MAX(INT(FLOAT(nmonsrec)/12.),1)
122    c--   Prevent code from crashing if integrated for less than a year
123          minrec = MIN(nmonsrec,12)
124    
125  c--   Loop over records.  c--   Loop over records.
126        do irec = 1,12        do irec = 1,minrec
127    
128  c--     Read time averages and the monthly mean data.  c--     Read time averages and the monthly mean data.
129          call active_read_xyz( fnamew, wbar, irec,          call active_read_xyz( fnamew, wbar, irec,
130       &                        doglobalread, ladinit,       &                        doglobalread, ladinit,
131       &                        optimcycle, mythid,       &                        optimcycle, mythid,
132       &                        xx_wbar_mean_dummy )       &                        xx_wbar_mean_dummy )
133        
134          nf = nf + 1          nf = nf + 1
135          do bj = jtlo,jthi          do bj = jtlo,jthi
136            do bi = itlo,ithi            do bi = itlo,ithi
137              do k = 1,nr                  do k = 1,nr
138                do j = jmin,jmax                do j = jmin,jmax
139                  do i = imin,imax                  do i = imin,imax
140                    wfmean(i,j,k,bi,bj) = wfmean(i,j,k,bi,bj) +                    wfmean(i,j,k,bi,bj) = wfmean(i,j,k,bi,bj) +
141       &                  wbar(i,j,k,bi,bj)       &                  wbar(i,j,k,bi,bj)
142                  enddo                  enddo
143                enddo                                enddo
144              enddo              enddo
145            enddo            enddo
146          enddo          enddo
147            
148        enddo        enddo
           
       do irec = nmonsrec-12+1, nmonsrec  
149    
150          do irec = nmonsrec-minrec+1, nmonsrec
151    
152  c--     Read time averages and the monthly mean data.  c--     Read time averages and the monthly mean data.
153          call active_read_xyz( fnamew, wbar, irec,          call active_read_xyz( fnamew, wbar, irec,
154       &                        doglobalread, ladinit,       &                        doglobalread, ladinit,
155       &                        optimcycle, mythid,       &                        optimcycle, mythid,
156       &                        xx_wbar_mean_dummy )       &                        xx_wbar_mean_dummy )
157        
158          nl = nl + 1            nl = nl + 1
159        
160          do bj = jtlo,jthi          do bj = jtlo,jthi
161            do bi = itlo,ithi            do bi = itlo,ithi
162              do k = 1,nr                  do k = 1,nr
163                do j = jmin,jmax                do j = jmin,jmax
164                  do i = imin,imax                  do i = imin,imax
165                    wlmean(i,j,k,bi,bj) = wlmean(i,j,k,bi,bj) +                    wlmean(i,j,k,bi,bj) = wlmean(i,j,k,bi,bj) +
166       &                  wbar(i,j,k,bi,bj)       &                  wbar(i,j,k,bi,bj)
167                  enddo                  enddo
168                enddo                                enddo
169              enddo              enddo
170            enddo            enddo
171          enddo          enddo
172          
173        enddo        enddo
174    
175            
176        do bj = jtlo,jthi        do bj = jtlo,jthi
177          do bi = itlo,ithi          do bi = itlo,ithi
178    
179  c--       Loop over the model layers  c--       Loop over the model layers
180            fctilew = 0. _d 0            fctilew = 0. _d 0
181              
182            do k = 1,nr            do k = 1,nr
183    
184  c--         Compute model misfit and cost function term for  c--         Compute model misfit and cost function term for
185  c           the vertical velovity field. The error is 1e-4 m/s.  c           the vertical velovity field. The error is 1e-4 m/s.
186              do j = jmin,jmax              do j = jmin,jmax
187                do i = imin,imax               do i = imin,imax
188                    if (_hFacC(i,j,k,bi,bj) .ne. 0.) then                if (_hFacC(i,j,k,bi,bj) .ne. 0.) then
189                       fctilew = fctilew +                       fctilew = fctilew +
190       &                  (2.5e11*cosphi(i,j,bi,bj)*       &                  (2.5e11*cosphi(i,j,bi,bj)*
191       &          (wlmean(i,j,k,bi,bj)/nl - wfmean(i,j,k,bi,bj)/nf)*       &          (wlmean(i,j,k,bi,bj)/nl - wfmean(i,j,k,bi,bj)/nf)*
192       &          (wlmean(i,j,k,bi,bj)/nl - wfmean(i,j,k,bi,bj)/nf))       &          (wlmean(i,j,k,bi,bj)/nl - wfmean(i,j,k,bi,bj)/nf))
193                    endif                       if ( cosphi(i,j,bi,bj) .ne. 0. )
194                enddo       &                 num_wdrift(bi,bj) = num_wdrift(bi,bj) + 1. _d 0
195                  endif
196                 enddo
197              enddo              enddo
198    
199            enddo            enddo
# Line 213  c--       Print cost function for each t Line 213  c--       Print cost function for each t
213       &                          SQUEEZE_RIGHT , mythid)       &                          SQUEEZE_RIGHT , mythid)
214            write(msgbuf,'(a,d22.15)')            write(msgbuf,'(a,d22.15)')
215       &        '     cost function (wvel) = ',       &        '     cost function (wvel) = ',
216       &        fctilet+fctiles       &        fcthread_wdrift
217            call print_message( msgbuf, standardmessageunit,            call print_message( msgbuf, standardmessageunit,
218       &                      SQUEEZE_RIGHT , mythid)       &                      SQUEEZE_RIGHT , mythid)
219            write(msgbuf,'(a)') ' '            write(msgbuf,'(a)') ' '
# Line 226  c--       Print cost function for each t Line 226  c--       Print cost function for each t
226    
227  #ifdef ECCO_VERBOSE  #ifdef ECCO_VERBOSE
228  c--     Print cost function for all tiles.  c--     Print cost function for all tiles.
229          _GLOBAL_SUM_R8( fcthread_tdrift , myThid )          _GLOBAL_SUM_R8( fcthread_wdrift , myThid )
         _GLOBAL_SUM_R8( fcthread_tdrifs , myThid )  
230          write(msgbuf,'(a)') ' '          write(msgbuf,'(a)') ' '
231          call print_message( msgbuf, standardmessageunit,          call print_message( msgbuf, standardmessageunit,
232       &                      SQUEEZE_RIGHT , mythid)       &                      SQUEEZE_RIGHT , mythid)
233          write(msgbuf,'(a,i8.8)')          write(msgbuf,'(a,i8.8)')
234       &    ' cost_Drift: irec = ',irec       &    ' cost_Driftw: irec = ',irec
235          call print_message( msgbuf, standardmessageunit,          call print_message( msgbuf, standardmessageunit,
236       &                      SQUEEZE_RIGHT , mythid)       &                      SQUEEZE_RIGHT , mythid)
237          write(msgbuf,'(a,a,d22.15)')          write(msgbuf,'(a,a,d22.15)')
238       &    '  cost function value',       &    '  cost function value',
239       &    ' (wvel) = ',fcthread_tdrifw       &    ' (wvel) = ',fcthread_wdrift
240          call print_message( msgbuf, standardmessageunit,          call print_message( msgbuf, standardmessageunit,
241       &                      SQUEEZE_RIGHT , mythid)       &                      SQUEEZE_RIGHT , mythid)
242          write(msgbuf,'(a)') ' '          write(msgbuf,'(a)') ' '
# Line 245  c--     Print cost function for all tile Line 244  c--     Print cost function for all tile
244       &                      SQUEEZE_RIGHT , mythid)       &                      SQUEEZE_RIGHT , mythid)
245  #endif  #endif
246    
 #else  
 c--   Do not enter the calculation of the drift contribution to  
 c--   the final cost function.  
   
       _BEGIN_MASTER( mythid )  
         write(msgbuf,'(a)') ' '  
         call print_message( msgbuf, standardmessageunit,  
      &                      SQUEEZE_RIGHT , mythid)  
         write(msgbuf,'(a,a)')  
      &    ' cost_Driftw: no contribution of drift between the first ',  
      &                 'and the last year to cost function.'  
         call print_message( msgbuf, standardmessageunit,  
      &                      SQUEEZE_RIGHT , mythid)  
         write(msgbuf,'(a,a,i9.8)')  
      &    ' cost_Driftw: number of records that would have',  
      &                ' been processed: ',nmonsrec  
         call print_message( msgbuf, standardmessageunit,  
      &                      SQUEEZE_RIGHT , mythid)  
         write(msgbuf,'(a)') ' '  
         call print_message( msgbuf, standardmessageunit,  
      &                      SQUEEZE_RIGHT , mythid)  
       _END_MASTER( mythid )  
247  #endif  #endif
248    
249        return        return
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22