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

Annotation of /MITgcm/pkg/ecco/cost_obcsw.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:07 2003 UTC (20 years, 6 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint52l_pre, checkpoint52e_pre, hrcube4, checkpoint52n_post, checkpoint52j_post, checkpoint53d_post, checkpoint54a_pre, checkpoint55c_post, checkpoint54e_post, checkpoint52e_post, checkpoint54a_post, checkpoint53c_post, checkpoint55d_pre, hrcube_1, branch-netcdf, checkpoint52d_pre, checkpoint52l_post, checkpoint52k_post, checkpoint52b_pre, checkpoint54b_post, checkpoint53b_pre, checkpoint55b_post, checkpoint54d_post, checkpoint52m_post, checkpoint55, checkpoint53a_post, checkpoint54, checkpoint54f_post, checkpoint53b_post, checkpoint52a_pre, checkpoint53, checkpoint52, checkpoint52d_post, checkpoint52a_post, checkpoint52b_post, checkpoint53g_post, checkpoint52f_post, checkpoint52c_post, ecco_c52_e35, hrcube5, checkpoint52i_post, checkpoint52j_pre, checkpoint53f_post, checkpoint55a_post, checkpoint53d_pre, checkpoint54c_post, checkpoint52i_pre, checkpoint51u_post, checkpoint52h_pre, checkpoint52f_pre, hrcube_2, hrcube_3
Branch point for: netcdf-sm0
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
2     #include "COST_CPPOPTIONS.h"
3     #ifdef ALLOW_OBCS
4     # include "OBCS_OPTIONS.h"
5     #endif
6    
7     subroutine cost_obcsw(
8     I myiter,
9     I mytime,
10     I startrec,
11     I endrec,
12     I mythid
13     & )
14    
15     c ==================================================================
16     c SUBROUTINE cost_obcsw
17     c ==================================================================
18     c
19     c o cost function contribution obc
20     c
21     c ==================================================================
22     c SUBROUTINE cost_obcsw
23     c ==================================================================
24    
25     implicit none
26    
27     c == global variables ==
28    
29     #include "EEPARAMS.h"
30     #include "SIZE.h"
31     #include "PARAMS.h"
32     #include "GRID.h"
33     #include "DYNVARS.h"
34     #ifdef ALLOW_OBCS
35     # include "OBCS.h"
36     #endif
37    
38     #include "cal.h"
39     #include "ecco_cost.h"
40     #include "ctrl.h"
41     #include "ctrl_dummy.h"
42     #include "optim.h"
43    
44     c == routine arguments ==
45    
46     integer myiter
47     _RL mytime
48     integer mythid
49     integer startrec
50     integer endrec
51    
52     c == local variables ==
53    
54     integer bi,bj
55     integer i,j,k
56     integer itlo,ithi
57     integer jtlo,jthi
58     integer jmin,jmax
59     integer imin,imax
60     integer irec
61     integer il
62     integer iobcs
63     integer ip1
64     integer nrec
65     integer ilfld
66     integer igg
67    
68     _RL fctile
69     _RL fcthread
70     _RL dummy
71     _RL gg
72     _RL tmpx
73     cgg(
74     _RL tmpfield (1-oly:sny+oly,nr,nsx,nsy)
75     _RL barofield (1-oly:sny+oly,nsx,nsy)
76     _RL maskyz (1-oly:sny+oly,nr,nsx,nsy)
77     _RL utop
78    
79     character*(80) fnamefld
80    
81     logical doglobalread
82     logical ladinit
83    
84     #ifdef ECCO_VERBOSE
85     character*(MAX_LEN_MBUF) msgbuf
86     #endif
87    
88     c == external functions ==
89    
90     integer ilnblnk
91     external ilnblnk
92    
93     c == end of interface ==
94    
95     jtlo = mybylo(mythid)
96     jthi = mybyhi(mythid)
97     itlo = mybxlo(mythid)
98     ithi = mybxhi(mythid)
99     jmin = 1
100     jmax = sny
101     imin = 1
102     imax = snx
103    
104     c-- Read tiled data.
105     doglobalread = .false.
106     ladinit = .false.
107    
108     c Number of records to be used.
109     nrec = endrec-startrec+1
110    
111     #ifdef ALLOW_OBCSW_COST_CONTRIBUTION
112    
113     ip1 = 1
114     fcthread = 0. _d 0
115    
116     #ifdef ECCO_VERBOSE
117     _BEGIN_MASTER( mythid )
118     write(msgbuf,'(a)') ' '
119     call print_message( msgbuf, standardmessageunit,
120     & SQUEEZE_RIGHT , mythid)
121     write(msgbuf,'(a)') ' '
122     call print_message( msgbuf, standardmessageunit,
123     & SQUEEZE_RIGHT , mythid)
124     write(msgbuf,'(a,i9.8)')
125     & ' cost_obcsw: number of records to process: ',nrec
126     call print_message( msgbuf, standardmessageunit,
127     & SQUEEZE_RIGHT , mythid)
128     write(msgbuf,'(a)') ' '
129     call print_message( msgbuf, standardmessageunit,
130     & SQUEEZE_RIGHT , mythid)
131     _END_MASTER( mythid )
132     #endif
133    
134     if (optimcycle .ge. 0) then
135     ilfld=ilnblnk( xx_obcsw_file )
136     write(fnamefld(1:80),'(2a,i10.10)')
137     & xx_obcsw_file(1:ilfld), '.', optimcycle
138     endif
139    
140     c-- Loop over records.
141     do irec = 1,nrec
142    
143     do bj = jtlo,jthi
144     do bi = itlo,ithi
145     do j = jmin,jmax
146     barofield(j,bi,bj) = 0. _d 0
147     enddo
148     enddo
149     enddo
150    
151     call active_read_yz( fnamefld, tmpfield, irec, doglobalread,
152     & ladinit, optimcycle, mythid
153     & , xx_obcsw_dummy )
154    
155     cgg Need to solve for iobcs would have been.
156     gg = (irec-1)/nobcs
157     igg = int(gg)
158     iobcs = irec - igg*nobcs
159    
160     call active_read_yz( 'maskobcsw', maskyz,
161     & iobcs,
162     & doglobalread, ladinit, 0,
163     & mythid, dummy )
164    
165     cgg Need to change xx field to baroclinic vel.
166     cgg Level Nr contains barotropic vel, remove it.
167     cgg The deepest wet level velocity should be
168     cgg computed in order to ensure zero integrated flux.
169     if (iobcs .eq. 3 .or. iobcs .eq. 4) then
170     c-- Loop over this thread's tiles.
171     do bj = jtlo,jthi
172     do bi = itlo,ithi
173     do j = jmin,jmax
174    
175     cgg The barotropic velocity is stored in level 1.
176     cgg Penalize for deviation in it.
177     barofield(j,bi,bj) = tmpfield(j,1,bi,bj)
178     & * maskyz(j,1,bi,bj)
179    
180     tmpfield(j,1,bi,bj) = 0.d0
181     utop = 0.d0
182    
183     do k = 1,Nr
184     cgg If cells are not full, this should be modified with hFac.
185     cgg
186     cgg The xx field (tmpfldxz) does not contain the velocity at the
187     cgg lowest wet level. This velocity is not independent; it must
188     cgg exactly balance the volume flux, since we are dealing with
189     cgg the baroclinic velocity structure..
190     utop = tmpfield(j,k,bi,bj)*
191     & maskyz(j,k,bi,bj) * delZ(k) + utop
192     enddo
193     tmpfield(j,1,bi,bj) = -utop / delZ(1)
194     enddo
195     enddo
196     enddo
197    
198     endif
199    
200     c-- Loop over this thread's tiles.
201     do bj = jtlo,jthi
202     do bi = itlo,ithi
203    
204     c-- Determine the weights to be used.
205     fctile = 0. _d 0
206    
207     if (iobcs .eq. 3 .or. iobcs .eq. 4) then
208     do j = jmin,jmax
209     i = OB_Iw(j,bi,bj)
210     tmpx = barofield(j,bi,bj)
211     fctile = fctile
212     & + wbaro*cosphi(i,j,bi,bj)
213     & *tmpx*tmpx
214     & *maskyz(j,1,bi,bj)
215     cgg print*,'fctile barotropic W',fctile
216     enddo
217     endif
218    
219     do k = 1, Nr
220     do j = jmin,jmax
221     i = OB_Iw(j,bi,bj)
222     cgg if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
223     tmpx = tmpfield(j,k,bi,bj)
224     fctile = fctile
225     & + wobcswLev(j,k,bi,bj,iobcs)*cosphi(i,j,bi,bj)
226     & *tmpx*tmpx
227     & *maskyz(j,k,bi,bj)
228     cgg endif
229     enddo
230     enddo
231    
232     objf_obcsw(bi,bj) = objf_obcsw(bi,bj) + fctile
233     fcthread = fcthread + fctile
234     enddo
235     enddo
236    
237     #ifdef ECCO_VERBOSE
238     c-- Print cost function for all tiles.
239     _GLOBAL_SUM_R8( fcthread , myThid )
240     write(msgbuf,'(a)') ' '
241     call print_message( msgbuf, standardmessageunit,
242     & SQUEEZE_RIGHT , mythid)
243     write(msgbuf,'(a,i8.8)')
244     & ' cost_obcsw: irec = ',irec
245     call print_message( msgbuf, standardmessageunit,
246     & SQUEEZE_RIGHT , mythid)
247     write(msgbuf,'(a,a,d22.15)')
248     & ' global cost function value',
249     & ' (obcsw) = ',fcthread
250     call print_message( msgbuf, standardmessageunit,
251     & SQUEEZE_RIGHT , mythid)
252     write(msgbuf,'(a)') ' '
253     call print_message( msgbuf, standardmessageunit,
254     & SQUEEZE_RIGHT , mythid)
255     #endif
256    
257     enddo
258     c-- End of loop over records.
259    
260     #endif
261    
262     return
263     end
264    
265    
266    
267    
268    
269    
270    

  ViewVC Help
Powered by ViewVC 1.1.22