/[MITgcm]/MITgcm/pkg/ctrl/ctrl_init_wet.F
ViewVC logotype

Contents of /MITgcm/pkg/ctrl/ctrl_init_wet.F

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


Revision 1.3 - (show annotations) (download)
Tue Nov 16 05:42:12 2004 UTC (19 years, 6 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint57b_post, checkpoint57c_pre, checkpoint56b_post, checkpoint56c_post, checkpoint57a_post, checkpoint57a_pre, checkpoint57, checkpoint56, checkpoint57c_post, checkpoint56a_post
Changes since 1.2: +2 -1 lines
More on dsvd vs. MITgcm interfacing
o handling of g_, ad, via admtlm_vector (mds...vector)
o use ctrl_pack/unpack for admtlm_vector I/O
o use optimcycle for dsvd iteration
o make sure norm is w.r.t. derived quantities

1 C
2 C $Header: /u/gcmpack/MITgcm/pkg/ctrl/ctrl_init_wet.F,v 1.2 2003/11/06 22:05:08 heimbach Exp $
3 C $Name: $
4
5 #include "CTRL_CPPOPTIONS.h"
6
7 subroutine ctrl_init_wet( mythid )
8
9 c ==================================================================
10 c SUBROUTINE ctrl_init_wet
11 c ==================================================================
12
13 implicit none
14
15 c == global variables ==
16
17 #include "EEPARAMS.h"
18 #include "SIZE.h"
19 #include "PARAMS.h"
20 #include "GRID.h"
21 #include "ctrl.h"
22
23 #ifdef ALLOW_OBCS_CONTROL
24 # include "OBCS.h"
25 #endif
26
27 c == routine arguments ==
28
29 integer mythid
30
31 c == local variables ==
32
33 integer bi,bj
34 integer i,j,k
35 integer itlo,ithi
36 integer jtlo,jthi
37 integer jmin,jmax
38 integer imin,imax
39 integer ntmp
40 integer iobcs
41 integer nwetc3d
42
43 _RL dummy
44
45 character*(80) ymaskobcs
46
47 c-- Set loop ranges.
48 jtlo = mybylo(mythid)
49 jthi = mybyhi(mythid)
50 itlo = mybxlo(mythid)
51 ithi = mybxhi(mythid)
52 jmin = 1
53 jmax = sny
54 imin = 1
55 imax = snx
56
57 c-- Determine the number of wet points in each tile:
58 c-- maskc, masks, and maskw.
59
60 c-- Initialise the counters.
61 do bj = jtlo,jthi
62 do bi = itlo,ithi
63 do k = 1,nr
64 nwetctile(bi,bj,k) = 0
65 nwetstile(bi,bj,k) = 0
66 nwetwtile(bi,bj,k) = 0
67 nwetvtile(bi,bj,k) = 0
68 enddo
69 enddo
70 enddo
71
72 #ifdef ALLOW_OBCS_CONTROL
73 c-- Initialise obcs counters.
74 do bj = jtlo,jthi
75 do bi = itlo,ithi
76 do k = 1,nr
77 do iobcs = 1,nobcs
78 #ifdef ALLOW_OBCSN_CONTROL
79 nwetobcsn(bi,bj,k,iobcs) = 0
80 #endif
81 #ifdef ALLOW_OBCSS_CONTROL
82 nwetobcss(bi,bj,k,iobcs) = 0
83 #endif
84 #ifdef ALLOW_OBCSW_CONTROL
85 nwetobcsw(bi,bj,k,iobcs) = 0
86 #endif
87 #ifdef ALLOW_OBCSE_CONTROL
88 nwetobcse(bi,bj,k,iobcs) = 0
89 #endif
90 enddo
91 enddo
92 enddo
93 enddo
94 #endif
95
96 c-- Count wet points on each tile.
97 do bj = jtlo,jthi
98 do bi = itlo,ithi
99 do k = 1,nr
100 do j = jmin,jmax
101 do i = imin,imax
102 c-- Center mask.
103 if (hFacC(i,j,k,bi,bj) .ne. 0.) then
104 nwetctile(bi,bj,k) = nwetctile(bi,bj,k) + 1
105 endif
106 c-- South mask.
107 if (maskS(i,j,k,bi,bj) .eq. 1.) then
108 nwetstile(bi,bj,k) = nwetstile(bi,bj,k) + 1
109 endif
110 c-- West mask.
111 if (maskW(i,j,k,bi,bj) .eq. 1.) then
112 nwetwtile(bi,bj,k) = nwetwtile(bi,bj,k) + 1
113 endif
114 #if (defined (ALLOW_EFLUXP0_CONTROL))
115 c-- Vertical mask.
116 if (hFacV(i,j,k,bi,bj) .ne. 0.) then
117 nwetvtile(bi,bj,k) = nwetvtile(bi,bj,k) + 1
118 endif
119 #endif
120 enddo
121 enddo
122 enddo
123 enddo
124 enddo
125
126 #ifdef ALLOW_OBCSN_CONTROL
127 c-- Count wet points at Northern boundary.
128 c-- mask conventions are adopted from obcs_apply_ts, obcs_apply_uv
129 ymaskobcs = 'maskobcsn'
130 call ctrl_mask_set_xz( 0, OB_Jn, nwetobcsn, ymaskobcs, mythid )
131 #endif
132
133 #ifdef ALLOW_OBCSS_CONTROL
134 c-- Count wet points at Southern boundary.
135 c-- mask conventions are adopted from obcs_apply_ts, obcs_apply_uv
136 ymaskobcs = 'maskobcss'
137 call ctrl_mask_set_xz( 1, OB_Js, nwetobcss, ymaskobcs, mythid )
138 #endif
139
140 #ifdef ALLOW_OBCSW_CONTROL
141 c-- Count wet points at Western boundary.
142 c-- mask conventions are adopted from obcs_apply_ts, obcs_apply_uv
143 ymaskobcs = 'maskobcsw'
144 call ctrl_mask_set_yz( 1, OB_Iw, nwetobcsw, ymaskobcs, mythid )
145 #endif
146
147 #ifdef ALLOW_OBCSE_CONTROL
148 c-- Count wet points at Eastern boundary.
149 c-- mask conventions are adopted from obcs_apply_ts, obcs_apply_uv
150 ymaskobcs = 'maskobcse'
151 call ctrl_mask_set_yz( 0, OB_Ie, nwetobcse, ymaskobcs, mythid )
152 #endif
153
154 _BEGIN_MASTER( mythid )
155 c-- Determine the total number of control variables.
156 nvartype = 0
157 nvarlength = 0
158 do i = 1,maxcvars
159 c
160 if ( ncvarindex(i) .ne. -1 ) then
161 nvartype = nvartype + 1
162 do bj = jtlo,jthi
163 do bi = itlo,ithi
164 do k = 1,ncvarnrmax(i)
165 if ( ncvargrd(i) .eq. 'c' ) then
166 nvarlength = nvarlength +
167 & ncvarrecs(i)*nwetctile(bi,bj,k)
168 else if ( ncvargrd(i) .eq. 's' ) then
169 nvarlength = nvarlength +
170 & ncvarrecs(i)*nwetstile(bi,bj,k)
171 else if ( ncvargrd(i) .eq. 'w' ) then
172 nvarlength = nvarlength +
173 & ncvarrecs(i)*nwetwtile(bi,bj,k)
174 else if ( ncvargrd(i) .eq. 'v' ) then
175 nvarlength = nvarlength +
176 & ncvarrecs(i)*nwetvtile(bi,bj,k)
177 else if ( ncvargrd(i) .eq. 'm' ) then
178 #ifdef ALLOW_OBCS_CONTROL
179 do iobcs = 1, nobcs
180 cgg This overcounts the number of o.b. control points by a factor of "nobcs".
181 cgg As an ad-hoc solution I've divided by nobcs everywhere.
182 if ( i .eq. 11 ) then
183 #ifdef ALLOW_OBCSN_CONTROL
184 nvarlength = nvarlength +
185 & (ncvarrecs(i)/nobcs)
186 & *nwetobcsn(bi,bj,k,iobcs)
187 #endif
188 else if ( i .eq. 12 ) then
189 #ifdef ALLOW_OBCSS_CONTROL
190 nvarlength = nvarlength +
191 & (ncvarrecs(i)/nobcs)
192 & *nwetobcss(bi,bj,k,iobcs)
193 #endif
194 else if ( i .eq. 13 ) then
195 #ifdef ALLOW_OBCSW_CONTROL
196 nvarlength = nvarlength +
197 & (ncvarrecs(i)/nobcs)
198 & *nwetobcsw(bi,bj,k,iobcs)
199 #endif
200 else if ( i .eq. 14 ) then
201 #ifdef ALLOW_OBCSE_CONTROL
202 nvarlength = nvarlength +
203 & (ncvarrecs(i)/nobcs)
204 & *nwetobcse(bi,bj,k,iobcs)
205 #endif
206 end if
207 enddo
208 #endif
209 else
210 print*,'ctrl_init: invalid grid location'
211 print*,' control variable = ',ncvarindex(i)
212 print*,' grid location = ',ncvargrd(i)
213 stop ' ... stopped in ctrl_init'
214 endif
215 enddo
216 enddo
217 enddo
218 endif
219 enddo
220
221 cph(
222 print *, 'ph-wet 1: nvarlength = ', nvarlength
223 print *, 'ph-wet 2: surface wet C = ', nwetctile(1,1,1)
224 print *, 'ph-wet 3: surface wet W = ', nwetwtile(1,1,1)
225 print *, 'ph-wet 4: surface wet S = ', nwetstile(1,1,1)
226 print *, 'ph-wet 4a:surface wet V = ', nwetvtile(1,1,1)
227 nwetc3d = 0
228 do k = 1, Nr
229 nwetc3d = nwetc3d + nwetctile(1,1,k)
230 end do
231 print *, 'ph-wet 5: 3D wet points = ', nwetc3d
232 do i = 1, maxcvars
233 print *, 'ph-wet 6: no recs for i = ', i, ncvarrecs(i)
234 end do
235 print *, 'ph-wet 7: ',
236 & 2*nwetc3d +
237 & ncvarrecs(3)*nwetctile(1,1,1) +
238 & ncvarrecs(4)*nwetctile(1,1,1) +
239 & ncvarrecs(5)*nwetwtile(1,1,1) +
240 & ncvarrecs(6)*nwetstile(1,1,1)
241 print *, 'ph-wet 8: ',
242 & 2*nwetc3d +
243 & ncvarrecs(7)*nwetctile(1,1,1) +
244 & ncvarrecs(8)*nwetctile(1,1,1) +
245 & ncvarrecs(9)*nwetwtile(1,1,1) +
246 & ncvarrecs(10)*nwetstile(1,1,1)
247 #ifdef ALLOW_OBCSN_CONTROL
248 print *, 'ph-wet 9: surface wet obcsn = '
249 & , nwetobcsn(1,1,1,1), nwetobcsn(1,1,1,2)
250 & , nwetobcsn(1,1,1,3), nwetobcsn(1,1,1,4)
251 #endif
252 #ifdef ALLOW_OBCSS_CONTROL
253 print *, 'ph-wet 10: surface wet obcss = '
254 & , nwetobcss(1,1,1,1), nwetobcss(1,1,1,2)
255 & , nwetobcss(1,1,1,3), nwetobcss(1,1,1,4)
256 #endif
257 #ifdef ALLOW_OBCSW_CONTROL
258 print *, 'ph-wet 11: surface wet obcsw = '
259 & , nwetobcsw(1,1,1,1), nwetobcsw(1,1,1,2)
260 & , nwetobcsw(1,1,1,3), nwetobcsw(1,1,1,4)
261 #endif
262 #ifdef ALLOW_OBCSE_CONTROL
263 print *, 'ph-wet 12: surface wet obcse = '
264 & , nwetobcse(1,1,1,1), nwetobcse(1,1,1,2)
265 & , nwetobcse(1,1,1,3), nwetobcse(1,1,1,4)
266 #endif
267 cph)
268
269 CALL GLOBAL_SUM_INT( nvarlength, myThid )
270
271 print *, 'ph-wet 13: global nvarlength vor k=', k, nvarlength
272
273 c
274 c Summation of wet point counters
275 c
276 do k = 1, nr
277
278 ntmp=0
279 do bj=1,nSy
280 do bi=1,nSx
281 ntmp=ntmp+nWetcTile(bi,bj,k)
282 enddo
283 enddo
284 CALL GLOBAL_SUM_INT( ntmp, myThid )
285 nWetcGlobal(k)=ntmp
286 print *, 'ph-wet 14a: global nWet... k=', k, ntmp
287
288 ntmp=0
289 do bj=1,nSy
290 do bi=1,nSx
291 ntmp=ntmp+nWetsTile(bi,bj,k)
292 enddo
293 enddo
294 CALL GLOBAL_SUM_INT( ntmp, myThid )
295 nWetsGlobal(k)=ntmp
296 print *, 'ph-wet 14b: global nWet... k=', k, ntmp
297
298 ntmp=0
299 do bj=1,nSy
300 do bi=1,nSx
301 ntmp=ntmp+nWetwTile(bi,bj,k)
302 enddo
303 enddo
304 CALL GLOBAL_SUM_INT( ntmp, myThid )
305 nWetwGlobal(k)=ntmp
306 print *, 'ph-wet 14c: global nWet... k=', k, ntmp
307
308 ntmp=0
309 do bj=1,nSy
310 do bi=1,nSx
311 ntmp=ntmp+nWetvTile(bi,bj,k)
312 enddo
313 enddo
314 CALL GLOBAL_SUM_INT( ntmp, myThid )
315 nWetvGlobal(k)=ntmp
316 print *, 'ph-wet 14d: global nWet... k=', k, ntmp
317
318 #ifdef ALLOW_OBCSN_CONTROL
319 do iobcs = 1, nobcs
320 ntmp=0
321 do bj=1,nSy
322 do bi=1,nSx
323 ntmp=ntmp+nwetobcsn(bi,bj,k,iobcs)
324 enddo
325 enddo
326 CALL GLOBAL_SUM_INT( ntmp, myThid )
327 nwetobcsnglo(k,iobcs)=ntmp
328 print *, 'ph-wet 15a: global nWet... k=', k, iobcs, ntmp
329 enddo
330 #endif
331 #ifdef ALLOW_OBCSS_CONTROL
332 do iobcs = 1, nobcs
333 ntmp=0
334 do bj=1,nSy
335 do bi=1,nSx
336 ntmp=ntmp+nwetobcss(bi,bj,k,iobcs)
337 enddo
338 enddo
339 CALL GLOBAL_SUM_INT( ntmp, myThid )
340 nwetobcssglo(k,iobcs)=ntmp
341 print *, 'ph-wet 15b: global nWet... k=', k, iobcs, ntmp
342 enddo
343 #endif
344 #ifdef ALLOW_OBCSW_CONTROL
345 do iobcs = 1, nobcs
346 ntmp=0
347 do bj=1,nSy
348 do bi=1,nSx
349 ntmp=ntmp+nwetobcsw(bi,bj,k,iobcs)
350 enddo
351 enddo
352 CALL GLOBAL_SUM_INT( ntmp, myThid )
353 nwetobcswglo(k,iobcs)=ntmp
354 print *, 'ph-wet 15c: global nWet... k=', k, iobcs, ntmp
355 enddo
356 #endif
357 #ifdef ALLOW_OBCSE_CONTROL
358 do iobcs = 1, nobcs
359 ntmp=0
360 do bj=1,nSy
361 do bi=1,nSx
362 ntmp=ntmp+nwetobcse(bi,bj,k,iobcs)
363 enddo
364 enddo
365 CALL GLOBAL_SUM_INT( ntmp, myThid )
366 nwetobcseglo(k,iobcs)=ntmp
367 print *, 'ph-wet 15d: global nWet... k=', k, iobcs, ntmp
368 enddo
369 #endif
370
371 enddo
372
373 print*, 'ctrl_init: no. of control variables: ', nvartype
374 print*, 'ctrl_init: control vector length: ', nvarlength
375 _END_MASTER( mythid )
376
377 c Set unit weight to 1
378 c
379 do bj=1,nSy
380 do bi=1,nSx
381 do k=1, nr
382 wunit(k,bi,bj) = 1. _d 0
383 enddo
384 enddo
385 enddo
386
387 c write masks and weights to files to be read by a master process
388 c
389 call active_write_xyz( 'hFacC', hFacC, 1, 0, mythid, dummy)
390 call active_write_xyz( 'maskC', maskC, 1, 0, mythid, dummy)
391 call active_write_xyz( 'maskW', maskW, 1, 0, mythid, dummy)
392 call active_write_xyz( 'maskS', maskS, 1, 0, mythid, dummy)
393 #if (defined (ALLOW_EFLUXP0_CONTROL))
394 call active_write_xyz( 'hFacV', hFacV, 1, 0, mythid, dummy)
395 #endif
396
397 return
398 end

  ViewVC Help
Powered by ViewVC 1.1.22