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

Contents of /MITgcm/pkg/ecco/cost_obcsw.F

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


Revision 1.1 - (show annotations) (download)
Thu Nov 6 22:10:07 2003 UTC (20 years, 7 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
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