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

Annotation of /MITgcm/pkg/ecco/cost_xbt.F

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


Revision 1.1 - (hide annotations) (download)
Thu Nov 6 22:10:08 2003 UTC (20 years, 8 months ago) by heimbach
Branch: MAIN
CVS Tags: branch-netcdf, checkpoint52b_pre, checkpoint52, checkpoint52a_post, checkpoint52b_post, checkpoint52c_post, ecco_c52_e35, checkpoint52a_pre, checkpoint51u_post
o merging from ecco-branch
o pkg/ecco now containes ecco-specific part of cost function
o top level routines the_main_loop, forward_step
  supersede those in model/src/
  previous input data.cost now in data.ecco
  (new namelist ecco_cost_nml)

1 heimbach 1.1 C $Header: /u/gcmpack/MITgcm/pkg/cost/Attic/cost_xbt.F,v 1.1.2.3 2002/04/04 10:58:59 heimbach Exp $
2    
3     #include "COST_CPPOPTIONS.h"
4    
5    
6     subroutine cost_XBT(
7     I myiter,
8     I mytime,
9     I mythid
10     & )
11    
12     c ==================================================================
13     c SUBROUTINE cost_XBT
14     c ==================================================================
15     c
16     c o Evaluate cost function contribution of XBT 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_XBT
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 "ecco_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_xbt
59     _RL fcthread_xbt
60     _RL www (1-olx:snx+olx,1-oly:sny+oly)
61     _RL tmpobs (1-olx:snx+olx,1-oly:sny+oly)
62     _RL cmask (1-olx:snx+olx,1-oly:sny+oly)
63     _RL spval
64     _RL ztop,rl_35,rl_0
65    
66     character*(80) fnametheta
67    
68     logical doglobalread
69     logical ladinit
70    
71     character*(MAX_LEN_MBUF) msgbuf
72    
73     c == external functions ==
74    
75     integer ilnblnk
76     external ilnblnk
77     _RL SW_PTMP
78     external SW_PTMP
79    
80     c == end of interface ==
81    
82     jtlo = mybylo(mythid)
83     jthi = mybyhi(mythid)
84     itlo = mybxlo(mythid)
85     ithi = mybxhi(mythid)
86     jmin = 1
87     jmax = sny
88     imin = 1
89     imax = snx
90    
91     spval = -2.
92     ztop = -.981*1.027
93     rl_35 = 35.0
94     rl_0 = 0.0
95    
96     c-- Read state record from global file.
97     doglobalread = .false.
98     ladinit = .false.
99    
100     #ifdef ALLOW_XBT_COST_CONTRIBUTION
101    
102     #ifdef ECCO_VERBOSE
103     write(msgbuf,'(a)') ' '
104     call print_message( msgbuf, standardmessageunit,
105     & SQUEEZE_RIGHT , mythid)
106     write(msgbuf,'(a,i8.8)')
107     & ' cost_XBT: number of records to process =', nmonsrec
108     call print_message( msgbuf, standardmessageunit,
109     & SQUEEZE_RIGHT , mythid)
110     write(msgbuf,'(a)') ' '
111     call print_message( msgbuf, standardmessageunit,
112     & SQUEEZE_RIGHT , mythid)
113     #endif
114    
115     if (optimcycle .ge. 0) then
116     ilu=ilnblnk( tbarfile )
117     write(fnametheta(1:80),'(2a,i10.10)') tbarfile(1:ilu),'.',
118     & optimcycle
119     else
120     print*
121     print*,' cost_XBT: optimcycle has a wrong value.'
122     print*,' optimcycle = ',optimcycle
123     print*
124     stop ' ... stopped in cost_XBT.'
125     endif
126    
127     fcthread_xbt = 0. _d 0
128    
129     c-- Loop over records.
130     do irec = 1,nmonsrec
131    
132     c-- Read time averages and the monthly mean data.
133     call active_read_xy( fnametheta, tbar, irec,
134     & doglobalread, ladinit,
135     & optimcycle, mythid
136     & , xx_tbar_mean_dummy )
137    
138     c-- Determine the record to be read on observations.
139    
140     call mdsreadfield( xbtfile, cost_iprec, 'RL', nr, xbtobs,
141     & irec, mythid)
142    
143     c-- Loop over this thread's tiles.
144     do bj = jtlo,jthi
145     do bi = itlo,ithi
146     c-- Loop over the model layers
147    
148     fctile_xbt = 0. _d 0
149    
150     do k = 1,nr
151    
152     c-- Determine the weights to be used.
153     do j = jmin,jmax
154     do i = imin,imax
155     cmask(i,j) = 1. _d 0
156     if (xbtobs(i,j,k,bi,bj) .eq. 0.) then
157     cmask(i,j) = 0. _d 0
158     endif
159    
160     if (xbtobs(i,j,k,bi,bj) .lt. spval) then
161     cmask(i,j) = 0. _d 0
162     endif
163    
164     cph(
165     cph print *, 'WARNING: SPECIFIC SETUP FOR ECCO'
166     cph below statement could be replaced by following
167     cph to make it independnet of Nr:
168     cph
169     cph if ( rC(K) .GT. -1000. ) then
170     cph)
171     c-- set cmask=0 in areas shallower than 1000m
172     if (_hFacC(i,j,13,bi,bj) .eq. 0.) then
173     cmask(i,j) = 0. _d 0
174     endif
175    
176     if (_hFacC(i,j,k,bi,bj) .eq. 0.) then
177     cmask(i,j) = 0. _d 0
178     endif
179    
180     enddo
181     enddo
182    
183     do j = jmin,jmax
184     do i = imin,imax
185     www(i,j) = cosphi(i,j,bi,bj)*cmask(i,j)
186     tmpobs(i,j) = SW_PTMP(rl_35,
187     $ xbtobs(i,j,k,bi,bj),ztop*rc(k),rl_0)
188     c-- The array ctdtobs contains CTD temperature.
189    
190     fctile_xbt = fctile_xbt +
191     & (wtheta2(i,j,k,bi,bj)*www(i,j))*
192     & (tbar(i,j,k,bi,bj)-tmpobs(i,j))*
193     & (tbar(i,j,k,bi,bj)-tmpobs(i,j))
194     enddo
195     enddo
196     c-- End of loop over layers.
197     enddo
198    
199     fcthread_xbt = fcthread_xbt + fctile_xbt
200     objf_xbt(bi,bj) = objf_xbt(bi,bj) + fctile_xbt
201    
202     #ifdef ECCO_VERBOSE
203     write(msgbuf,'(a)') ' '
204     call print_message( msgbuf, standardmessageunit,
205     & SQUEEZE_RIGHT , mythid)
206     write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
207     & ' COST_XBT: irec,bi,bj = ',irec,bi,bj
208     call print_message( msgbuf, standardmessageunit,
209     & SQUEEZE_RIGHT , mythid)
210     write(msgbuf,'(a,d22.15)')
211     & ' COST_XBT: cost function = ', fctile_xbt
212     call print_message( msgbuf, standardmessageunit,
213     & SQUEEZE_RIGHT , mythid)
214     write(msgbuf,'(a)') ' '
215     call print_message( msgbuf, standardmessageunit,
216     & SQUEEZE_RIGHT , mythid)
217     #endif
218    
219     enddo
220     enddo
221    
222     #ifdef ECCO_VERBOSE
223     c-- Print cost function for all tiles.
224     _GLOBAL_SUM_R8( fcthread_xbt , myThid )
225     write(msgbuf,'(a)') ' '
226     call print_message( msgbuf, standardmessageunit,
227     & SQUEEZE_RIGHT , mythid)
228     write(msgbuf,'(a,i8.8)')
229     & ' cost_XBT: irec = ',irec
230     call print_message( msgbuf, standardmessageunit,
231     & SQUEEZE_RIGHT , mythid)
232     write(msgbuf,'(a,a,d22.15)')
233     & ' global cost function value',
234     & ' ( XBT temp. ) = ',fcthread_xbt
235     call print_message( msgbuf, standardmessageunit,
236     & SQUEEZE_RIGHT , mythid)
237     write(msgbuf,'(a)') ' '
238     call print_message( msgbuf, standardmessageunit,
239     & SQUEEZE_RIGHT , mythid)
240     #endif
241    
242     enddo
243     c-- End of second loop over records.
244    
245     #else
246     c-- Do not enter the calculation of the XBT temperature contribution
247     c-- to the final cost function.
248    
249     fctile_xbt = 0. _d 0
250     fcthread_xbt = 0. _d 0
251    
252     crg
253     nrec = 1
254     crg
255    
256     _BEGIN_MASTER( mythid )
257     write(msgbuf,'(a)') ' '
258     call print_message( msgbuf, standardmessageunit,
259     & SQUEEZE_RIGHT , mythid)
260     write(msgbuf,'(a,a)')
261     & ' cost_XBT: no contribution of XBT temperature ',
262     & ' to cost function.'
263     call print_message( msgbuf, standardmessageunit,
264     & SQUEEZE_RIGHT , mythid)
265     write(msgbuf,'(a,a,i9.8)')
266     & ' cost_XBT: number of records that would have',
267     & ' been processed: ',nrec
268     call print_message( msgbuf, standardmessageunit,
269     & SQUEEZE_RIGHT , mythid)
270     write(msgbuf,'(a)') ' '
271     call print_message( msgbuf, standardmessageunit,
272     & SQUEEZE_RIGHT , mythid)
273     _END_MASTER( mythid )
274     #endif
275    
276     return
277     end

  ViewVC Help
Powered by ViewVC 1.1.22