/[MITgcm]/MITgcm/pkg/cost/cost_ctds.F
ViewVC logotype

Diff of /MITgcm/pkg/cost/cost_ctds.F

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

revision 1.1 by heimbach, Tue Feb 5 20:23:57 2002 UTC revision 1.1.2.2 by heimbach, Fri Feb 15 22:27:47 2002 UTC
# Line 0  Line 1 
1    C $Header$
2    
3    #include "COST_CPPOPTIONS.h"
4    
5    
6          subroutine cost_CTDS(
7         I                     myiter,
8         I                     mytime,
9         I                     mythid
10         &                   )
11    
12    c     ==================================================================
13    c     SUBROUTINE cost_CTDS
14    c     ==================================================================
15    c
16    c     o Evaluate cost function contribution of CTD temperature data.
17    c
18    c     started:  Elisabeth Remy eremy@ucsd.edu 30-Aug-2000
19    c
20    c
21    c     ==================================================================
22    c     SUBROUTINE cost_CTDS
23    c     ==================================================================
24    
25          implicit none
26    
27    c     == global variables ==
28    
29    #include "EEPARAMS.h"
30    #include "SIZE.h"
31    #include "GRID.h"
32    #include "DYNVARS.h"
33    
34    #include "cal.h"
35    #include "cost.h"
36    #include "ctrl.h"
37    #include "ctrl_dummy.h"
38    #include "optim.h"
39    
40    c     == routine arguments ==
41    
42          integer myiter
43          _RL     mytime
44          integer mythid
45    
46    c     == local variables ==
47    
48          integer bi,bj
49          integer i,j,k
50          integer itlo,ithi
51          integer jtlo,jthi
52          integer jmin,jmax
53          integer imin,imax
54          integer nrec
55          integer irec
56          integer ilu
57          
58          _RL fctile_ctds
59          _RL fcthread_ctds
60          _RL www (1-olx:snx+olx,1-oly:sny+oly)
61          _RL wtmp (1-olx:snx+olx,1-oly:sny+oly)
62          _RL tmpobs (1-olx:snx+olx,1-oly:sny+oly)
63          _RL tmpbar (1-olx:snx+olx,1-oly:sny+oly)
64          _RL cmask (1-olx:snx+olx,1-oly:sny+oly)  
65          _RL spval
66          
67          character*(80) fnamesalt
68    
69          logical doglobalread
70          logical ladinit
71    
72          character*(MAX_LEN_MBUF) msgbuf
73    
74    c     == external functions ==
75    
76          integer  ilnblnk
77          external ilnblnk
78    
79    c     == end of interface ==
80    
81          jtlo = mybylo(mythid)
82          jthi = mybyhi(mythid)
83          itlo = mybxlo(mythid)
84          ithi = mybxhi(mythid)
85          jmin = 1
86          jmax = sny
87          imin = 1
88          imax = snx
89          
90          spval = -9990.
91    
92    c--   Read state record from global file.
93          doglobalread = .false.
94          ladinit      = .false.
95          
96    c      call cal_TimePassed( topexstartdate, middate, difftime, mythid )
97    c      call cal_ToSeconds( difftime, diffsecs, mythid )
98    
99    #ifdef ALLOW_CTDS_COST_CONTRIBUTION
100    
101    #ifdef ECCO_VERBOSE
102          write(msgbuf,'(a)') ' '
103          call print_message( msgbuf, standardmessageunit,
104         &                    SQUEEZE_RIGHT , mythid)
105          write(msgbuf,'(a,i8.8)')
106         &  ' cost_CTDS: number of records to process =', nmonsrec
107          call print_message( msgbuf, standardmessageunit,
108         &                    SQUEEZE_RIGHT , mythid)
109          write(msgbuf,'(a)') ' '
110          call print_message( msgbuf, standardmessageunit,
111         &                    SQUEEZE_RIGHT , mythid)
112    #endif
113    
114          if (optimcycle .ge. 0) then
115            ilu=ilnblnk( sbarfile )
116            write(fnamesalt(1:80),'(2a,i10.10)') sbarfile(1:ilu),'.',
117         &                                   optimcycle
118          else
119            print*
120            print*,' cost_CTDS: optimcycle has a wrong value.'
121            print*,'                 optimcycle = ',optimcycle
122            print*
123            stop   '  ... stopped in cost_CTDS.'
124          endif
125          
126          fcthread_ctds = 0. _d 0
127    
128    c--   Loop over records.
129          do irec = 1,nmonsrec
130    
131    c--     Read time averages and the monthly mean data.
132            call active_read_xy( fnamesalt, sbar, irec,
133         &          doglobalread, ladinit,
134         &          optimcycle, mythid
135         &   , xx_sbar_mean_dummy )
136      
137    ce     &                       myiter, mytime )
138    
139    c--     Determine the record to be read on observations.
140    
141    c        levoff = mod(modelstartdate(1)/100,100)
142    c        levmon = (irec-1) + levoff
143    c        levmon = mod(levmon-1,12)+1
144                  
145    ccc        call cost_ReadCtd( mythid )
146            call mdsreadfield( ctdsfile, 64, 'RL', nr, ctdsobs, irec,
147         &                     mythid)
148    
149    c--     Loop over this thread's tiles.
150            do bj = jtlo,jthi
151              do bi = itlo,ithi
152    c--         Loop over the model layers
153    
154                fctile_ctds = 0. _d 0
155    
156                do k = 1,nr
157    
158    c--           Determine the weights to be used.
159                  do j = jmin,jmax
160                    do i = imin,imax
161                      cmask(i,j) = 1. _d 0
162                      if (ctdsobs(i,j,k,bi,bj) .eq. 0.) then
163                        cmask(i,j) = 0. _d 0
164                      endif
165                    
166                      if (ctdsobs(i,j,k,bi,bj) .lt. spval) then
167                        cmask(i,j) = 0. _d 0
168                      endif
169                      
170                    enddo
171                  enddo
172                  
173                  do j = jmin,jmax
174                    do i = imin,imax
175                      www(i,j)    = cosphi(i,j,bi,bj)*cmask(i,j)                
176                      tmpobs(i,j) = ctdsobs(i,j,k,bi,bj)
177                      tmpbar(i,j) = sbar(i,j,k,bi,bj)
178                      wtmp(i,j) = wsalt(k,bi,bj)
179                    enddo
180                  enddo
181    
182                  do j = jmin,jmax
183                    do i = imin,imax
184    c--               The array ctdsobs contains CTD salinity.
185                      if(cmask(i,j).eq.1.0) then
186                        print *,'cost_ctds',i,j,tmpbar(i,j),tmpobs(i,j)
187                      endif
188                      fctile_ctds = fctile_ctds +
189         &                             (wtmp(i,j)*www(i,j))*
190         &                             (tmpbar(i,j)-tmpobs(i,j))*
191         &                             (tmpbar(i,j)-tmpobs(i,j))
192    
193                    enddo
194                  enddo
195                enddo
196    c--         End of loop over layers.
197    
198                fcthread_ctds    = fcthread_ctds    + fctile_ctds
199                objf_ctds(bi,bj) = objf_ctds(bi,bj) + fctile_ctds  
200    
201    #ifdef ECCO_VERBOSE
202                write(msgbuf,'(a)') ' '
203                call print_message( msgbuf, standardmessageunit,
204         &                          SQUEEZE_RIGHT , mythid)
205                write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
206         &        ' COST_CTDS: irec,bi,bj            =  ',irec,bi,bj
207                call print_message( msgbuf, standardmessageunit,
208         &                          SQUEEZE_RIGHT , mythid)
209                write(msgbuf,'(a,d22.15)')
210         &        ' COST_CTDS: cost function         = ', fctile_ctds
211                call print_message( msgbuf, standardmessageunit,
212         &                      SQUEEZE_RIGHT , mythid)
213                write(msgbuf,'(a)') ' '
214                call print_message( msgbuf, standardmessageunit,
215         &                          SQUEEZE_RIGHT , mythid)
216    #endif
217    
218             enddo
219            enddo
220    
221    #ifdef ECCO_VERBOSE
222    c--     Print cost function for all tiles.
223            _GLOBAL_SUM_R8( fcthread_ctds , myThid )
224            write(msgbuf,'(a)') ' '
225            call print_message( msgbuf, standardmessageunit,
226         &                      SQUEEZE_RIGHT , mythid)
227            write(msgbuf,'(a,i8.8)')
228         &    ' cost_CTDS: irec =  ',irec
229            call print_message( msgbuf, standardmessageunit,
230         &                      SQUEEZE_RIGHT , mythid)
231            write(msgbuf,'(a,a,d22.15)')
232         &    ' global cost function value',
233         &    ' ( CTD sal. )  = ',fcthread_ctds
234            call print_message( msgbuf, standardmessageunit,
235         &                      SQUEEZE_RIGHT , mythid)
236            write(msgbuf,'(a)') ' '
237            call print_message( msgbuf, standardmessageunit,
238         &                      SQUEEZE_RIGHT , mythid)
239    #endif
240    
241          enddo
242    c--   End of second loop over records.
243    
244    #else
245    c--   Do not enter the calculation of the CTD temperature contribution
246    c--   to the final cost function.
247          
248          fctile_ctds   = 0. _d 0
249          fcthread_ctds = 0. _d 0
250    
251    crg
252          nrec = 1
253    crg
254    
255          _BEGIN_MASTER( mythid )
256            write(msgbuf,'(a)') ' '
257            call print_message( msgbuf, standardmessageunit,
258         &                      SQUEEZE_RIGHT , mythid)
259            write(msgbuf,'(a,a)')
260         &    ' cost_CTDS: no contribution of CTD temperature ',
261         &                    ' to cost function.'
262            call print_message( msgbuf, standardmessageunit,
263         &                      SQUEEZE_RIGHT , mythid)
264            write(msgbuf,'(a,a,i9.8)')
265         &    ' cost_CDTS: number of records that would have',
266         &                      ' been processed: ',nrec
267            call print_message( msgbuf, standardmessageunit,
268         &                      SQUEEZE_RIGHT , mythid)
269            write(msgbuf,'(a)') ' '
270            call print_message( msgbuf, standardmessageunit,
271         &                      SQUEEZE_RIGHT , mythid)
272          _END_MASTER( mythid )
273    #endif
274    
275          return
276          end
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.1.2.2
  ViewVC Help
Powered by ViewVC 1.1.22