/[MITgcm]/MITgcm_contrib/verification_other/global1x1_tot/code/ecco_cost_weights.F
ViewVC logotype

Contents of /MITgcm_contrib/verification_other/global1x1_tot/code/ecco_cost_weights.F

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

Revision 1.1 - (show annotations) (download)
Sat Feb 4 02:45:01 2012 UTC (13 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64k, checkpoint64x, checkpoint64z, checkpoint64o, checkpoint64p, checkpoint64r, checkpoint64w, checkpoint64v, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint64m, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65u, checkpoint65j, checkpoint67a, checkpoint67b, checkpoint65i, checkpoint67d, checkpoint65m, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65e, checkpoint64i, checkpoint64h, checkpoint65h, checkpoint65, checkpoint64j, checkpoint65n, HEAD 
move un-tested set-up from MITgcm/verification/global1x1_tot to Contrib/verification_other
1 C $Header: /u/gcmpack/MITgcm/verification/global1x1_tot/code/ecco_cost_weights.F,v 1.5 2007/10/09 02:29:10 jmc Exp $
2 C $Name: $
3
4 #include "COST_CPPOPTIONS.h"
5
6
7 subroutine ecco_cost_weights( mythid )
8
9 c ==================================================================
10 c SUBROUTINE ecco_cost_weights
11 c ==================================================================
12 c
13 c o Read the weights used for the cost function evaluation.
14 c
15 c started: Christian Eckert eckert@mit.edu 30-Jun-1999
16 c
17 c changed: Christian Eckert eckert@mit.edu 25-Feb-2000
18 c
19 c - Restructured the code in order to create a package
20 c for the MITgcmUV.
21 c
22 c Christian Eckert eckert@mit.edu 02-May-2000
23 c
24 c - corrected typo in mdsreadfield( sflux_errfile );
25 c wp --> wsflux. Spotted by Patrick Heimbach.
26 c
27 c ==================================================================
28 c SUBROUTINE ecco_cost_weights
29 c ==================================================================
30
31 implicit none
32
33 c == global variables ==
34
35 #include "EEPARAMS.h"
36 #include "SIZE.h"
37 #include "PARAMS.h"
38 #include "GRID.h"
39
40 #include "ctrl.h"
41 #include "ecco_cost.h"
42
43 c == routine arguments ==
44
45 integer mythid
46
47 c == local variables ==
48
49 integer bi,bj
50 integer i,j,k
51 integer itlo,ithi
52 integer jtlo,jthi
53 integer jmin,jmax
54 integer imin,imax
55 integer gwunit
56 integer irec,nnz
57 integer ilo,ihi
58 integer iobcs
59
60 _RL factor
61 _RL wti(nr)
62 _RL wsi(nr)
63 _RL wui(nr)
64 _RL wvi(nr)
65 _RL whflux0m
66 _RL wsflux0m
67 _RL wtau0m
68 _RL ratio
69 _RL dummy
70
71 c == external ==
72
73 integer ifnblnk
74 external ifnblnk
75 integer ilnblnk
76 external ilnblnk
77
78 c == end of interface ==
79
80 jtlo = mybylo(mythid)
81 jthi = mybyhi(mythid)
82 itlo = mybxlo(mythid)
83 ithi = mybxhi(mythid)
84 jmin = 1-oly
85 jmax = sny+oly
86 imin = 1-olx
87 imax = snx+olx
88
89 c-- Initialize background weights
90 wtau0 = 0.
91 whflux0 = 0.
92 wsflux0 = 0.
93 whflux0m = 0
94 wsflux0m = 0.
95 watemp0 = 0.
96 waqh0 = 0.
97 wwind0 = 0.
98
99 c-- Initialize variance (weight) fields.
100 do k = 1,nr
101 wti(k) = 0. _d 0
102 wsi(k) = 0. _d 0
103 wui(k) = 0. _d 0
104 wvi(k) = 0. _d 0
105 enddo
106 do bj = jtlo,jthi
107 do bi = itlo,ithi
108 do j = jmin,jmax
109 do i = imin,imax
110 whflux (i,j,bi,bj) = 0. _d 0
111 whfluxm (i,j,bi,bj) = 0. _d 0
112 wsflux (i,j,bi,bj) = 0. _d 0
113 wsfluxm (i,j,bi,bj) = 0. _d 0
114 wtauu (i,j,bi,bj) = 0. _d 0
115 wtauum (i,j,bi,bj) = 0. _d 0
116 wtauv (i,j,bi,bj) = 0. _d 0
117 wtauvm (i,j,bi,bj) = 0. _d 0
118 watemp (i,j,bi,bj) = 0. _d 0
119 waqh (i,j,bi,bj) = 0. _d 0
120 wuwind (i,j,bi,bj) = 0. _d 0
121 wvwind (i,j,bi,bj) = 0. _d 0
122 wsst (i,j,bi,bj) = 0. _d 0
123 wsss (i,j,bi,bj) = 0. _d 0
124 wtp (i,j,bi,bj) = 0. _d 0
125 wers (i,j,bi,bj) = 0. _d 0
126 wgfo (i,j,bi,bj) = 0. _d 0
127 wp (i,j,bi,bj) = 0. _d 0
128 wudrift (i,j,bi,bj) = 0. _d 0
129 wvdrift (i,j,bi,bj) = 0. _d 0
130 cph(
131 whflux2 (i,j,bi,bj) = 0. _d 0
132 wsflux2 (i,j,bi,bj) = 0. _d 0
133 wtauu2 (i,j,bi,bj) = 0. _d 0
134 wtauv2 (i,j,bi,bj) = 0. _d 0
135 cph)
136 enddo
137 enddo
138 enddo
139 enddo
140 do bj = jtlo,jthi
141 do bi = itlo,ithi
142 do k = 1,Nr
143 wtheta (k,bi,bj) = 0. _d 0
144 wsalt (k,bi,bj) = 0. _d 0
145 wctdt (k,bi,bj) = 0. _d 0
146 wctds (k,bi,bj) = 0. _d 0
147 do j = jmin,jmax
148 do i = imin,imax
149 wtheta2 (i,j,k,bi,bj) = 0. _d 0
150 wsalt2 (i,j,k,bi,bj) = 0. _d 0
151 wthetaLev (i,j,k,bi,bj) = 0. _d 0
152 wsaltLev (i,j,k,bi,bj) = 0. _d 0
153 enddo
154 enddo
155 enddo
156 enddo
157 enddo
158
159 #if (defined (ALLOW_OBCS_COST_CONTRIBUTION) || \
160 defined (ALLOW_OBCS_CONTROL))
161 do iobcs = 1,nobcs
162 do k = 1,Nr
163 #if (defined (ALLOW_OBCSN_CONTROL) || \
164 defined (ALLOW_OBCSN_COST_CONTRIBUTION))
165 wobcsn(k,iobcs) = 0. _d 0
166 #endif
167 #if (defined (ALLOW_OBCSS_CONTROL) || \
168 defined (ALLOW_OBCSS_COST_CONTRIBUTION))
169 wobcss(k,iobcs) = 0. _d 0
170 #endif
171 #if (defined (ALLOW_OBCSW_CONTROL) || \
172 defined (ALLOW_OBCSW_COST_CONTRIBUTION))
173 wobcsw(k,iobcs) = 0. _d 0
174 #endif
175 #if (defined (ALLOW_OBCSE_CONTROL) || \
176 defined (ALLOW_OBCSE_COST_CONTRIBUTION))
177 wobcse(k,iobcs) = 0. _d 0
178 #endif
179 enddo
180 enddo
181 #endif
182
183 c-- Build area weighting matrix used in the cost function
184 c-- contributions.
185
186 c-- Define frame.
187 do j = jmin,jmax
188 do i = imin,imax
189 c-- North/South and West/East edges set to zero.
190 if ( (j .lt. 1) .or. (j .gt. sny) .or.
191 & (i .lt. 1) .or. (i .gt. snx) ) then
192 frame(i,j) = 0. _d 0
193 else
194 frame(i,j) = 1. _d 0
195 endif
196 enddo
197 enddo
198
199 c-- First account for the grid used.
200 if (usingCartesianGrid) then
201 factor = 0. _d 0
202 else if (usingSphericalPolarGrid) then
203 factor = 1. _d 0
204 endif
205
206 do bj = jtlo,jthi
207 do bi = itlo,ithi
208 do j = jmin,jmax
209 do i = imin,imax
210 cds cosphi(i,j,bi,bj) = cos(yc(i,j,bi,bj)*deg2rad*factor)*
211 cds & frame(i,j)
212 cosphi(i,j,bi,bj) = frame(i,j)
213 enddo
214 enddo
215 enddo
216 enddo
217
218 c-- Read error information and set up weight matrices.
219 _BEGIN_MASTER(myThid)
220 ilo = ifnblnk(data_errfile)
221 ihi = ilnblnk(data_errfile)
222 CALL OPEN_COPY_DATA_FILE(
223 I data_errfile(ilo:ihi),
224 I 'ECCO_COST_WEIGHTS',
225 O gwunit,
226 I myThid )
227
228 read(gwunit,*)
229 #if (defined (ALLOW_HFLUX_COST_CONTRIBUTION) || defined (ALLOW_HFLUX_CONTROL))
230 & whflux0
231 #elif (defined (ALLOW_ATEMP_COST_CONTRIBUTION) || defined (ALLOW_ATEMP_CONTROL))
232 & watemp0
233 #endif
234 #if (defined (ALLOW_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SFLUX_CONTROL))
235 & , wsflux0
236 #elif (defined (ALLOW_AQH_COST_CONTRIBUTION) || defined (ALLOW_AQH_CONTROL))
237 & , waqh0
238 #endif
239 #if (defined (ALLOW_USTRESS_COST_CONTRIBUTION) || defined (ALLOW_USTRESS_CONTROL))
240 & , wtau0
241 #elif (defined (ALLOW_UWIND_COST_CONTRIBUTION) || defined (ALLOW_UWIND_CONTROL))
242 & , wwind0
243 #endif
244 & , ratio
245 #if (defined (ALLOW_OBCS_COST_CONTRIBUTION) || defined (ALLOW_OBCS_CONTROL))
246 & , wbaro
247 #endif
248
249 do k = 1,nr
250 read(gwunit,*) wti(k), wsi(k)
251 #if (defined (ALLOW_OBCS_COST_CONTRIBUTION) || \
252 defined (ALLOW_OBCS_CONTROL))
253 & , wvi(k)
254 #endif
255 end do
256 close(gwunit)
257
258 whflux0m = whflux0
259 wsflux0m = wsflux0
260 wtau0m = wtau0
261 _END_MASTER(myThid)
262
263 _BARRIER
264
265 cph jmin = 1
266 cph jmax = sny
267 cph imin = 1
268 cph imax = snx
269
270 do bj = jtlo,jthi
271 do bi = itlo,ithi
272
273 wsfluxmm(bi,bj) = 1.
274 whfluxmm(bi,bj) = 1.
275
276 c-- The "classic" state estimation tool wastes memory here;
277 c-- as long as there is not more information available there
278 c-- is no need to add the zonal and meridional directions.
279 do k = 1,nr
280 wtheta(k,bi,bj) = wti(k)
281 wsalt (k,bi,bj) = wsi(k)
282 wcurrent(k,bi,bj) = wvi(k)
283 enddo
284
285 do k = 1,nr
286 #ifdef ALLOW_OBCSN_COST_CONTRIBUTION
287 wobcsn(k,1) = wti(k)
288 wobcsn(k,2) = wsi(k)
289 wobcsn(k,3) = wti(k)*0.02
290 wobcsn(k,4) = wti(k)*0.02
291 #endif
292 #ifdef ALLOW_OBCSS_COST_CONTRIBUTION
293 wobcss(k,1) = wti(k)
294 wobcss(k,2) = wsi(k)
295 wobcss(k,3) = wti(k)*0.02
296 wobcss(k,4) = wti(k)*0.02
297 #endif
298 #ifdef ALLOW_OBCSW_COST_CONTRIBUTION
299 wobcsw(k,1) = wti(k)
300 wobcsw(k,2) = wsi(k)
301 wobcsw(k,3) = wti(k)*0.02
302 wobcsw(k,4) = wti(k)*0.02
303 #endif
304 #ifdef ALLOW_OBCSE_COST_CONTRIBUTION
305 wobcse(k,1) = wti(k)
306 wobcse(k,2) = wsi(k)
307 wobcse(k,3) = wti(k)*0.02
308 wobcse(k,4) = wti(k)*0.02
309 #endif
310 enddo
311 enddo
312 enddo
313
314 #ifdef ALLOW_SALT0_COST_CONTRIBUTION
315 if ( salt0errfile .NE. ' ' ) then
316 call mdsreadfield( salt0errfile, 32, 'RL', Nr,
317 & wsaltLev, 1, mythid)
318 do bj = jtlo,jthi
319 do bi = itlo,ithi
320 do k = 1,nr
321 do j = jmin,jmax
322 do i = imin,imax
323 c-- Test for missing values.
324 if ( wsaltLev(i,j,k,bi,bj).eq.0 ) then
325 wsaltLev(i,j,k,bi,bj) = 0. _d 0
326 else
327 wsaltLev(i,j,k,bi,bj)=frame(i,j)/
328 $ ( wsaltLev(i,j,k,bi,bj)*wsaltLev(i,j,k,bi,bj) )
329 endif
330 enddo
331 enddo
332 enddo
333 enddo
334 enddo
335 else
336 do bj = jtlo,jthi
337 do bi = itlo,ithi
338 do k = 1,nr
339 do j = jmin,jmax
340 do i = imin,imax
341 wsaltLev(i,j,k,bi,bj)=ratio/(wsalt(k,bi,bj)
342 $ *wsalt(k,bi,bj)) *frame(i,j)
343 enddo
344 enddo
345 enddo
346 enddo
347 enddo
348 endif
349 call active_write_xyz( 'wsaltLev', wsaltLev,
350 & 1, 0, mythid, dummy)
351 #endif
352
353 #ifdef ALLOW_THETA0_COST_CONTRIBUTION
354 if ( temp0errfile .NE. ' ' ) then
355 call mdsreadfield( temp0errfile, 32, 'RL', Nr,
356 & wthetaLev, 1, mythid)
357 do bj = jtlo,jthi
358 do bi = itlo,ithi
359 do k = 1,nr
360 do j = jmin,jmax
361 do i = imin,imax
362 c-- Test for missing values.
363 if ( wthetaLev(i,j,k,bi,bj).eq.0 ) then
364 wthetaLev(i,j,k,bi,bj) = 0. _d 0
365 else
366 wthetaLev(i,j,k,bi,bj)=frame(i,j)/
367 $ ( wthetaLev(i,j,k,bi,bj)*wthetaLev(i,j,k,bi,bj) )
368 endif
369 enddo
370 enddo
371 enddo
372 enddo
373 enddo
374 else
375 do bj = jtlo,jthi
376 do bi = itlo,ithi
377 do k = 1,nr
378 do j = jmin,jmax
379 do i = imin,imax
380 wthetaLev(i,j,k,bi,bj)=ratio/(wtheta(k,bi,bj)
381 $ *wtheta(k,bi,bj)) *frame(i,j)
382 enddo
383 enddo
384 enddo
385 enddo
386 enddo
387 endif
388 call active_write_xyz( 'wthetaLev', wthetaLev,
389 & 1, 0, mythid, dummy)
390 #endif
391
392 #if (defined (ALLOW_ARGO_SALT_COST_CONTRIBUTION) || \
393 defined (ALLOW_CTDS_COST_CONTRIBUTION)|| \
394 defined (ALLOW_CTDSCLIM_COST_CONTRIBUTION))
395
396 if ( salterrfile .NE. ' ' ) then
397 call mdsreadfield( salterrfile, 32, 'RL', Nr, wsalt2, 1,
398 & mythid)
399
400 do bj = jtlo,jthi
401 do bi = itlo,ithi
402 do k = 1,nr
403 do j = jmin,jmax
404 do i = imin,imax
405 c-- Test for missing values.
406 if (wsalt(k,bi,bj).eq.0 .and.
407 $ wsalt2(i,j,k,bi,bj).eq.0) then
408 wsalt2(i,j,k,bi,bj) = 0. _d 0
409 else
410 cph new weights by G. Forget dont need MAX
411 wsalt2(i,j,k,bi,bj)=frame(i,j)/(
412 $ wsalt2(i,j,k,bi,bj)*wsalt2(i,j,k,bi,bj) )
413 endif
414 enddo
415 enddo
416 enddo
417 enddo
418 enddo
419 else
420 do bj = jtlo,jthi
421 do bi = itlo,ithi
422 do k = 1,nr
423 do j = jmin,jmax
424 do i = imin,imax
425 wsalt2(i,j,k,bi,bj)=ratio/(wsalt(k,bi,bj)
426 $ *wsalt(k,bi,bj)) *frame(i,j)
427 enddo
428 enddo
429 enddo
430 enddo
431 enddo
432 endif
433 #endif
434
435 #if (defined (ALLOW_ARGO_THETA_COST_CONTRIBUTION) || \
436 defined (ALLOW_CTDT_COST_CONTRIBUTION) || \
437 defined (ALLOW_CTDTCLIM_COST_CONTRIBUTION) || \
438 defined (ALLOW_XBT_COST_CONTRIBUTION))
439
440 if ( temperrfile .NE. ' ' ) then
441 call mdsreadfield( temperrfile, 32, 'RL', Nr, wtheta2, 1,
442 & mythid)
443 do bj = jtlo,jthi
444 do bi = itlo,ithi
445 do k = 1,nr
446 do j = jmin,jmax
447 do i = imin,imax
448 c-- Test for missing values.
449 if (wtheta(k,bi,bj).eq.0 .and.
450 $ wtheta2(i,j,k,bi,bj).eq.0) then
451 wtheta2(i,j,k,bi,bj) = 0. _d 0
452 else
453 cph new weights by G. Forget dont need MAX
454 wtheta2(i,j,k,bi,bj)=frame(i,j)/(
455 $ wtheta2(i,j,k,bi,bj)*wtheta2(i,j,k,bi,bj) )
456 endif
457 enddo
458 enddo
459 enddo
460 enddo
461 enddo
462 else
463 do bj = jtlo,jthi
464 do bi = itlo,ithi
465 do k = 1,nr
466 do j = jmin,jmax
467 do i = imin,imax
468 if (wtheta(k,bi,bj).eq.0 ) then
469 wtheta2(i,j,k,bi,bj) = 0. _d 0
470 else
471 wtheta2(i,j,k,bi,bj) = ratio/(wtheta(k,bi,bj)
472 $ *wtheta(k,bi,bj))*frame(i,j)
473 endif
474 enddo
475 enddo
476 enddo
477 enddo
478 enddo
479 endif
480 #endif
481
482 k = 1
483 do bj = jtlo,jthi
484 do bi = itlo,ithi
485 do j = jmin,jmax
486 do i = imin,imax
487 if (_hFacC(i,j,k,bi,bj) .eq. 0.) then
488 wsst(i,j,bi,bj) = 0. _d 0
489 wsss(i,j,bi,bj) = 0. _d 0
490 else
491 cph wsst(i,j,bi,bj) = wtheta(k,bi,bj)*10.
492 cph wsss(i,j,bi,bj) = wsalt(k,bi,bj)*10.
493 cph factor 5. by D Stammer
494 wsst(i,j,bi,bj) = wtheta(k,bi,bj)*frame(i,j)
495 wsss(i,j,bi,bj) = wsalt(k,bi,bj)*frame(i,j)
496 endif
497 enddo
498 enddo
499 enddo
500 enddo
501 #ifdef ALLOW_EGM96_ERROR_DIAG
502 c-- Read egm-96 geoid covariance. Data in units of meters.
503 nnz = 1
504 irec = 1
505 call mdsreadfield( geoid_errfile, cost_iprec, cost_yftype, nnz,
506 & wp, irec, mythid )
507 c-- Set all tile edges to zero.
508 do bj = jtlo,jthi
509 do bi = itlo,ithi
510 do j = jmin,jmax
511 do i = imin,imax
512 wp(i,j,bi,bj) = wp(i,j,bi,bj)*frame(i,j)
513 cph-indonesian(
514 if ( xC(i,j,bi,bj) .GT. 115. .AND.
515 & xC(i,j,bi,bj) .LT. 130. .AND.
516 & yC(i,j,bi,bj) .GT. -10. .AND.
517 & yC(i,j,bi,bj) .LT. 10. ) then
518 cph wp(i,j,bi,bj) = wp(i,j,bi,bj)*10000.
519 wp(i,j,bi,bj) = 0.
520 endif
521 cph-indonesian)
522 cph-medit(
523 if ( ( xC(i,j,bi,bj) .GT. 355. .AND.
524 & xC(i,j,bi,bj) .LT. 360. .AND.
525 & yC(i,j,bi,bj) .GT. 30. .AND.
526 & yC(i,j,bi,bj) .LT. 48. )
527 & .OR.
528 & ( xC(i,j,bi,bj) .GT. 0. .AND.
529 & xC(i,j,bi,bj) .LT. 39. .AND.
530 & yC(i,j,bi,bj) .GT. 30. .AND.
531 & yC(i,j,bi,bj) .LT. 48. ) ) then
532 wp(i,j,bi,bj) = wp(i,j,bi,bj)*10.
533 endif
534 cph-medit)
535 enddo
536 enddo
537 enddo
538 enddo
539 #else
540 do bj = jtlo,jthi
541 do bi = itlo,ithi
542 do j = jmin,jmax
543 do i = imin,imax
544 wp(i,j,bi,bj) = frame(i,j)
545 enddo
546 enddo
547 enddo
548 enddo
549 #endif
550
551 #ifdef ALLOW_SSH_COST_CONTRIBUTION
552 c-- Read T/P SSH anomaly rms field. Data in units of centimeters.
553 nnz = 1
554 irec = 1
555 call mdsreadfield( ssh_errfile, cost_iprec, cost_yftype, nnz,
556 & wtp, irec, mythid )
557
558 do bj = jtlo,jthi
559 do bi = itlo,ithi
560 k = 1
561 do j = jmin,jmax
562 do i = imin,imax
563 c-- Unit conversion to meters. ERS error is set to
564 c-- T/P error + 5cm
565 if (_hFacC(i,j,k,bi,bj) .eq. 0.) then
566 wtp (i,j,bi,bj) = 0. _d 0
567 wers(i,j,bi,bj) = 0. _d 0
568 wgfo(i,j,bi,bj) = 0. _d 0
569 else
570 wtp (i,j,bi,bj) = ( wtp(i,j,bi,bj) * 0.01 * 0.5 )
571 & *frame(i,j)
572 wers(i,j,bi,bj) = ( wtp(i,j,bi,bj) + 0.05 )
573 & *frame(i,j)
574 wgfo(i,j,bi,bj) = ( wtp(i,j,bi,bj) + 0.05 )
575 & *frame(i,j)
576 endif
577 cph-indonesian(
578 if ( xC(i,j,bi,bj) .GT. 115. .AND.
579 & xC(i,j,bi,bj) .LT. 130. .AND.
580 & yC(i,j,bi,bj) .GT. -10. .AND.
581 & yC(i,j,bi,bj) .LT. 10. ) then
582 wtp(i,j,bi,bj) = 0.
583 wers(i,j,bi,bj) = 0.
584 wgfo(i,j,bi,bj) = 0.
585 endif
586 cph-indonesian)
587 cph-medit(
588 if ( ( xC(i,j,bi,bj) .GT. 355. .AND.
589 & xC(i,j,bi,bj) .LT. 360. .AND.
590 & yC(i,j,bi,bj) .GT. 30. .AND.
591 & yC(i,j,bi,bj) .LT. 48. )
592 & .OR.
593 & ( xC(i,j,bi,bj) .GT. 0. .AND.
594 & xC(i,j,bi,bj) .LT. 39. .AND.
595 & yC(i,j,bi,bj) .GT. 30. .AND.
596 & yC(i,j,bi,bj) .LT. 48. ) ) then
597 wtp(i,j,bi,bj) = wtp(i,j,bi,bj) *10.
598 wers(i,j,bi,bj) = wers(i,j,bi,bj)*10.
599 wgfo(i,j,bi,bj) = wgfo(i,j,bi,bj)*10.
600 endif
601 cph-medit)
602 enddo
603 enddo
604 enddo
605 enddo
606 #endif /* ALLOW_SSH_COST_CONTRIBUTION */
607
608 c-- Read zonal wind stress variance.
609 #if (defined (ALLOW_SCAT_COST_CONTRIBUTION))
610
611 nnz = 1
612 irec = 1
613 call mdsreadfield( scatx_errfile, cost_iprec, cost_yftype, nnz,
614 & wscatx, irec, mythid )
615 call mdsreadfield( scaty_errfile, cost_iprec, cost_yftype, nnz,
616 & wscaty, irec, mythid )
617
618 do bj = jtlo,jthi
619 do bi = itlo,ithi
620 k = 1
621 do j = jmin,jmax
622 do i = imin,imax
623 c-- Test for missing values.
624 if (wscatx(i,j,bi,bj) .lt. -9900.) then
625 wscatx(i,j,bi,bj) = 0. _d 0
626 endif
627 c-- Rescale dyn -> N/M^2
628 wscatx(i,j,bi,bj) = wscatx(i,j,bi,bj)
629 c-- Missing values over water should have larger errors
630 if ( wscatx(i,j,bi,bj).EQ.0. .AND.
631 & maskW(i,j,k,bi,bj).NE.0. )
632 & wscatx(i,j,bi,bj) = 4.*wtau0
633 c-- Cut off extreme values
634 if ( wscatx(i,j,bi,bj).GT.0.15 )
635 & wscatx(i,j,bi,bj) = 0.15
636 c-- Set mimimum background
637 wscatx(i,j,bi,bj) = max(wscatx(i,j,bi,bj),wtau0)
638 wscatx(i,j,bi,bj) = wscatx(i,j,bi,bj)*maskW(i,j,k,bi,bj)
639 & *frame(i,j)
640 c
641 if (wscaty(i,j,bi,bj) .lt. -9900.) then
642 wscaty(i,j,bi,bj) = 0. _d 0
643 endif
644 c-- Rescale dyn -> N/M^2
645 wscaty(i,j,bi,bj) = wscaty(i,j,bi,bj)
646 c-- Missing values over water should have larger errors
647 if ( wscaty(i,j,bi,bj).EQ.0. .AND.
648 & maskS(i,j,k,bi,bj).NE.0. )
649 & wscaty(i,j,bi,bj) = 4.*wtau0
650 c-- Cut off extreme values
651 if ( wscaty(i,j,bi,bj).GT.0.15 )
652 & wscaty(i,j,bi,bj) = 0.15
653 c-- Set mimimum background
654 wscaty(i,j,bi,bj) = max(wscaty(i,j,bi,bj),wtau0)
655 wscaty(i,j,bi,bj) = wscaty(i,j,bi,bj)*maskS(i,j,k,bi,bj)
656 & *frame(i,j)
657 enddo
658 enddo
659 enddo
660 enddo
661
662 #endif
663
664 c-- Read zonal wind stress variance.
665 #if (defined (ALLOW_STRESS_MEAN_COST_CONTRIBUTION))
666 nnz = 1
667 irec = 1
668 cph call mdsreadfield( tauum_errfile, cost_iprec, cost_yftype,
669 cph & nnz, wtauum, irec, mythid )
670 cph call mdsreadfield( tauvm_errfile, cost_iprec, cost_yftype,
671 cph & nnz, wtauvm, irec, mythid )
672
673 do bj = jtlo,jthi
674 do bi = itlo,ithi
675 k = 1
676 do j = jmin,jmax
677 do i = imin,imax
678 c-- Test for missing values.
679 if (wtauum(i,j,bi,bj) .lt. -9900.) then
680 wtauum(i,j,bi,bj) = 0. _d 0
681 endif
682 wtauum(i,j,bi,bj) = max(wtauum(i,j,bi,bj),wtau0m)
683 & *frame(i,j)
684 c-- Test for missing values.
685 if (wtauvm(i,j,bi,bj) .lt. -9900.) then
686 wtauvm(i,j,bi,bj) = 0. _d 0
687 endif
688 wtauvm(i,j,bi,bj) = max(wtauvm(i,j,bi,bj),wtau0m)
689 & *frame(i,j)
690 enddo
691 enddo
692 enddo
693 enddo
694 #endif
695
696 #if (defined (ALLOW_USTRESS_COST_CONTRIBUTION))
697 nnz = 1
698 ce irec = 2
699 ce( due to Patrick's processing:
700 irec = 1
701 ce)
702 if ( tauu_errfile .NE. ' ' ) then
703 call mdsreadfield( tauu_errfile, cost_iprec, cost_yftype, nnz,
704 & wtauu, irec, mythid )
705 endif
706
707 do bj = jtlo,jthi
708 do bi = itlo,ithi
709 k = 1
710 do j = jmin,jmax
711 do i = imin,imax
712 c-- Test for missing values.
713 if (wtauu(i,j,bi,bj) .lt. -9900.) then
714 wtauu(i,j,bi,bj) = 0. _d 0
715 endif
716 c-- Rescale dyn -> N/M^2
717 wtauu(i,j,bi,bj) = wtauu(i,j,bi,bj)*0.1
718 c-- Missing values over water should have larger errors
719 if ( wtauu(i,j,bi,bj).EQ.0. .AND.
720 & maskW(i,j,k,bi,bj).NE.0. )
721 & wtauu(i,j,bi,bj) = 4.*wtau0
722 c-- Cut off extreme values
723 if ( wtauu(i,j,bi,bj).GT.0.12 )
724 & wtauu(i,j,bi,bj) = 0.12
725 c-- Set mimimum background
726 wtauu(i,j,bi,bj) = max(wtauu(i,j,bi,bj),wtau0)
727 wtauu(i,j,bi,bj) = wtauu(i,j,bi,bj)*maskW(i,j,k,bi,bj)
728 & *frame(i,j)
729 cph(
730 cph wtauu2(i,j,bi,bj) = 2.*wtau0*maskW(i,j,k,bi,bj)*frame(i,j)
731 wtauu2(i,j,bi,bj) = wtauu(i,j,bi,bj)
732 cph)
733 enddo
734 enddo
735 enddo
736 enddo
737
738 #elif (defined (ALLOW_UWIND_COST_CONTRIBUTION))
739
740 nnz = 1
741 ce irec = 2
742 ce( due to Patrick's processing:
743 irec = 1
744 ce)
745 if ( uwind_errfile .NE. ' ' ) then
746 call mdsreadfield( uwind_errfile, cost_iprec, cost_yftype, nnz,
747 & wuwind, irec, mythid )
748 endif
749
750 do bj = jtlo,jthi
751 do bi = itlo,ithi
752 k = 1
753 do j = jmin,jmax
754 do i = imin,imax
755 c-- Test for missing values.
756 if (wuwind(i,j,bi,bj) .lt. -9900.) then
757 wuwind(i,j,bi,bj) = 0. _d 0
758 endif
759 wuwind(i,j,bi,bj) = wuwind(i,j,bi,bj)
760 wuwind(i,j,bi,bj) = max(wuwind(i,j,bi,bj),wwind0)
761 wuwind(i,j,bi,bj) = wuwind(i,j,bi,bj)*maskC(i,j,k,bi,bj)
762 & *frame(i,j)
763 enddo
764 enddo
765 enddo
766 enddo
767 #endif
768
769 c-- Read meridional wind stress variance.
770 #if (defined (ALLOW_VSTRESS_COST_CONTRIBUTION))
771 nnz = 1
772 ce irec = 2
773 ce( due to Patrick's processing:
774 irec = 1
775 ce)
776
777 if ( tauv_errfile .NE. ' ' ) then
778 call mdsreadfield( tauv_errfile, cost_iprec, cost_yftype, nnz,
779 & wtauv, irec, mythid )
780 endif
781
782 do bj = jtlo,jthi
783 do bi = itlo,ithi
784 do j = jmin,jmax
785 do i = imin,imax
786 c-- Test for missing values.
787 if (wtauv(i,j,bi,bj) .lt. -9900.) then
788 wtauv(i,j,bi,bj) = 0. _d 0
789 endif
790 c-- Rescape dyn -> dyn
791 wtauv(i,j,bi,bj) = wtauv(i,j,bi,bj)*0.1
792 c-- Missing values over water should have larger errors
793 if ( wtauv(i,j,bi,bj).EQ.0. .AND.
794 & maskS(i,j,k,bi,bj).NE.0. )
795 & wtauv(i,j,bi,bj) = 4.*wtau0
796 c-- Cut off extreme values
797 if ( wtauv(i,j,bi,bj).GT.0.12 )
798 & wtauv(i,j,bi,bj) = 0.12
799 c-- Set mimimum background
800 wtauv(i,j,bi,bj) = max(wtauv(i,j,bi,bj),wtau0)
801 wtauv(i,j,bi,bj) = wtauv(i,j,bi,bj)*maskS(i,j,k,bi,bj)
802 & *frame(i,j)
803 cph(
804 cph wtauv2(i,j,bi,bj) = 2.*wtau0*maskS(i,j,k,bi,bj)*frame(i,j)
805 wtauv2(i,j,bi,bj) = wtauv(i,j,bi,bj)
806 cph)
807 enddo
808 enddo
809 enddo
810 enddo
811
812 #elif (defined (ALLOW_VWIND_COST_CONTRIBUTION))
813
814 nnz = 1
815 ce irec = 2
816 ce( due to Patrick's processing:
817 irec = 1
818 ce)
819
820 if ( vwind_errfile .NE. ' ' ) then
821 call mdsreadfield( vwind_errfile, cost_iprec, cost_yftype, nnz,
822 & wvwind, irec, mythid )
823 endif
824
825 do bj = jtlo,jthi
826 do bi = itlo,ithi
827 do j = jmin,jmax
828 do i = imin,imax
829 c-- Test for missing values.
830 if (wvwind(i,j,bi,bj) .lt. -9900.) then
831 wvwind(i,j,bi,bj) = 0. _d 0
832 endif
833 wvwind(i,j,bi,bj) = wvwind(i,j,bi,bj)
834 wvwind(i,j,bi,bj) = max(wvwind(i,j,bi,bj),wwind0)
835 wvwind(i,j,bi,bj) = wvwind(i,j,bi,bj)*maskC(i,j,k,bi,bj)
836 & *frame(i,j)
837 enddo
838 enddo
839 enddo
840 enddo
841 #endif
842
843 #if (defined (ALLOW_HFLUX_COST_CONTRIBUTION))
844 c-- Read heat flux flux variance.
845 nnz = 1
846 c-- First record in data file: mean field.
847 c-- Second record in data file: rms field.
848 ce irec = 2
849 ce( due to Patrick's processing:
850 irec = 1
851 ce)
852 if ( hflux_errfile .NE. ' ' ) then
853 call mdsreadfield( hflux_errfile, cost_iprec, cost_yftype, nnz,
854 & whflux, irec, mythid )
855 endif
856
857 do bj = jtlo,jthi
858 do bi = itlo,ithi
859 do j = jmin,jmax
860 do i = imin,imax
861 c-- Test for missing values.
862 if (whflux(i,j,bi,bj) .lt. -9900.) then
863 whflux(i,j,bi,bj) = 0. _d 0
864 endif
865 c-- Data are in units of W/m**2.
866 whflux(i,j,bi,bj) = whflux(i,j,bi,bj)/3.
867 whflux(i,j,bi,bj) = max(whflux(i,j,bi,bj),whflux0)
868 & *frame(i,j)
869 whfluxm(i,j,bi,bj) = max(whfluxm(i,j,bi,bj),whflux0m)
870 & *frame(i,j)
871 cph(
872 cph whflux2(i,j,bi,bj) = 2.*whflux0*frame(i,j)
873 whflux2(i,j,bi,bj) = whflux(i,j,bi,bj)
874 cph)
875 enddo
876 enddo
877 enddo
878 enddo
879 #elif (defined (ALLOW_ATEMP_COST_CONTRIBUTION))
880 c-- Read atmos. temp. variance.
881 nnz = 1
882 c-- First record in data file: mean field.
883 c-- Second record in data file: rms field.
884 ce irec = 2
885 ce( due to Patrick's processing:
886 irec = 1
887 ce)
888 if ( atemp_errfile .NE. ' ' ) then
889 call mdsreadfield( atemp_errfile, cost_iprec, cost_yftype, nnz,
890 & watemp, irec, mythid )
891 endif
892
893 do bj = jtlo,jthi
894 do bi = itlo,ithi
895 do j = jmin,jmax
896 do i = imin,imax
897 c-- Test for missing values.
898 if (watemp(i,j,bi,bj) .lt. -9900.) then
899 watemp(i,j,bi,bj) = 0. _d 0
900 endif
901 c-- Data are in units of deg.
902 watemp(i,j,bi,bj) = watemp(i,j,bi,bj)
903 watemp(i,j,bi,bj) = max(watemp(i,j,bi,bj),watemp0)
904 & *frame(i,j)
905 enddo
906 enddo
907 enddo
908 enddo
909 #endif
910
911 #if (defined (ALLOW_SFLUX_COST_CONTRIBUTION))
912 c-- Read salt flux variance. Second read: data in units of m/s.
913 nnz = 1
914 c-- First record in data file: mean field.
915 c-- Second record in data file: rms field.
916 ce irec = 2
917 ce( due to Patrick's processing:
918 irec = 1
919 ce)
920 if ( sflux_errfile .NE. ' ' ) then
921 call mdsreadfield( sflux_errfile, cost_iprec, cost_yftype, nnz,
922 & wsflux, irec, mythid )
923 endif
924
925 do bj = jtlo,jthi
926 do bi = itlo,ithi
927 do j = jmin,jmax
928 do i = imin,imax
929 c-- Test for missing values.
930 if (wsflux(i,j,bi,bj) .lt. -9900.) then
931 wsflux(i,j,bi,bj) = 0. _d 0
932 endif
933 c-- Data are in units of m/s.
934 wsflux(i,j,bi,bj) = wsflux(i,j,bi,bj) / 3.
935 wsflux(i,j,bi,bj) = max(wsflux(i,j,bi,bj),wsflux0)
936 & *frame(i,j)
937 wsfluxm(i,j,bi,bj) = max(wsfluxm(i,j,bi,bj),wsflux0m)
938 & *frame(i,j)
939 cph(
940 cph wsflux2(i,j,bi,bj) = 2.*wsflux0*frame(i,j)
941 wsflux2(i,j,bi,bj) = wsflux(i,j,bi,bj)
942 cph)
943 enddo
944 enddo
945 enddo
946 enddo
947 #elif (defined (ALLOW_AQH_COST_CONTRIBUTION))
948 c-- Secific humid. variance. Second read: data in units of m/s.
949 nnz = 1
950 c-- First record in data file: mean field.
951 c-- Second record in data file: rms field.
952 ce irec = 2
953 ce( due to Patrick's processing:
954 irec = 1
955 ce)
956 if ( aqh_errfile .NE. ' ' ) then
957 call mdsreadfield( aqh_errfile, cost_iprec, cost_yftype, nnz,
958 & waqh, irec, mythid )
959 endif
960
961 do bj = jtlo,jthi
962 do bi = itlo,ithi
963 do j = jmin,jmax
964 do i = imin,imax
965 c-- Test for missing values.
966 if (waqh(i,j,bi,bj) .lt. -9900.) then
967 waqh(i,j,bi,bj) = 0. _d 0
968 endif
969 c-- Data are in units of
970 waqh(i,j,bi,bj) = waqh(i,j,bi,bj)
971 waqh(i,j,bi,bj) = max(waqh(i,j,bi,bj),waqh0)
972 & *frame(i,j)
973 enddo
974 enddo
975 enddo
976 enddo
977 #endif
978
979 c-- Units have to be checked!
980 do bj = jtlo,jthi
981 do bi = itlo,ithi
982 do j = jmin,jmax
983 do i = imin,imax
984 if (wtp(i,j,bi,bj) .ne. 0.) then
985 wtp (i,j,bi,bj) = 1./wtp(i,j,bi,bj)/wtp(i,j,bi,bj)
986 endif
987 if (wers(i,j,bi,bj) .ne. 0.) then
988 wers(i,j,bi,bj) = 1./wers(i,j,bi,bj)/wers(i,j,bi,bj)
989 endif
990 if (wgfo(i,j,bi,bj) .ne. 0.) then
991 wgfo(i,j,bi,bj) = 1./wgfo(i,j,bi,bj)/wgfo(i,j,bi,bj)
992 endif
993 cph(
994 cph sst, sss: reduce weights by factor 2
995 if (wsst(i,j,bi,bj) .ne. 0.) then
996 wsst(i,j,bi,bj) = 1./wsst(i,j,bi,bj)/wsst(i,j,bi,bj)/2.
997 endif
998 if (wsss(i,j,bi,bj) .ne. 0.) then
999 wsss(i,j,bi,bj) = 1./wsss(i,j,bi,bj)/wsss(i,j,bi,bj)/2.
1000 endif
1001 cph)
1002 if (wp(i,j,bi,bj) .ne. 0.) then
1003 wp(i,j,bi,bj) = 1./wp(i,j,bi,bj)/wp(i,j,bi,bj)
1004 endif
1005 if (wtauu(i,j,bi,bj) .ne. 0.) then
1006 wtauu(i,j,bi,bj) =
1007 & 1./wtauu(i,j,bi,bj)/wtauu(i,j,bi,bj)
1008 else
1009 wtauu(i,j,bi,bj) = 0.0 _d 0
1010 endif
1011 if (wtauum(i,j,bi,bj) .ne. 0.) then
1012 wtauum(i,j,bi,bj) =
1013 & 1./wtauum(i,j,bi,bj)/wtauum(i,j,bi,bj)
1014 else
1015 wtauum(i,j,bi,bj) = 0.0 _d 0
1016 endif
1017 if (wscatx(i,j,bi,bj) .ne. 0.) then
1018 wscatx(i,j,bi,bj) =
1019 & 1./wscatx(i,j,bi,bj)/wscatx(i,j,bi,bj)
1020 else
1021 wscatx(i,j,bi,bj) = 0.0 _d 0
1022 endif
1023 if (wtauv(i,j,bi,bj) .ne. 0.) then
1024 wtauv(i,j,bi,bj) =
1025 & 1./wtauv(i,j,bi,bj)/wtauv(i,j,bi,bj)
1026 else
1027 wtauv(i,j,bi,bj) = 0.0 _d 0
1028 endif
1029 if (wtauvm(i,j,bi,bj) .ne. 0.) then
1030 wtauvm(i,j,bi,bj) =
1031 & 1./wtauvm(i,j,bi,bj)/wtauvm(i,j,bi,bj)
1032 else
1033 wtauvm(i,j,bi,bj) = 0.0 _d 0
1034 endif
1035 if (wscaty(i,j,bi,bj) .ne. 0.) then
1036 wscaty(i,j,bi,bj) =
1037 & 1./wscaty(i,j,bi,bj)/wscaty(i,j,bi,bj)
1038 else
1039 wscaty(i,j,bi,bj) = 0.0 _d 0
1040 endif
1041 if (whflux(i,j,bi,bj) .ne. 0.) then
1042 whflux(i,j,bi,bj) =
1043 & 1./whflux(i,j,bi,bj)/whflux(i,j,bi,bj)
1044 else
1045 whflux(i,j,bi,bj) = 0.0 _d 0
1046 endif
1047 if (whfluxm(i,j,bi,bj) .ne. 0.) then
1048 whfluxm(i,j,bi,bj) =
1049 & 1./whfluxm(i,j,bi,bj)/whfluxm(i,j,bi,bj)
1050 else
1051 whfluxm(i,j,bi,bj) = 0.0 _d 0
1052 endif
1053 if (wsflux(i,j,bi,bj) .ne. 0.) then
1054 wsflux(i,j,bi,bj) =
1055 & 1./wsflux(i,j,bi,bj)/wsflux(i,j,bi,bj)
1056 else
1057 wsflux(i,j,bi,bj) = 0.0 _d 0
1058 endif
1059 if (wsfluxm(i,j,bi,bj) .ne. 0.) then
1060 wsfluxm(i,j,bi,bj) =
1061 & 1./wsfluxm(i,j,bi,bj)/wsfluxm(i,j,bi,bj)
1062 else
1063 wsfluxm(i,j,bi,bj) = 0.0 _d 0
1064 endif
1065 cph)
1066 if (wuwind(i,j,bi,bj) .ne. 0.) then
1067 wuwind(i,j,bi,bj) =
1068 & 1./wuwind(i,j,bi,bj)/wuwind(i,j,bi,bj)
1069 else
1070 wuwind(i,j,bi,bj) = 0.0 _d 0
1071 endif
1072 if (wvwind(i,j,bi,bj) .ne. 0.) then
1073 wvwind(i,j,bi,bj) =
1074 & 1./wvwind(i,j,bi,bj)/wvwind(i,j,bi,bj)
1075 else
1076 wvwind(i,j,bi,bj) = 0.0 _d 0
1077 endif
1078 if (watemp(i,j,bi,bj) .ne. 0.) then
1079 watemp(i,j,bi,bj) =
1080 & 1./watemp(i,j,bi,bj)/watemp(i,j,bi,bj)
1081 else
1082 watemp(i,j,bi,bj) = 0.0 _d 0
1083 endif
1084 if (waqh(i,j,bi,bj) .ne. 0.) then
1085 waqh(i,j,bi,bj) =
1086 & 1./waqh(i,j,bi,bj)/waqh(i,j,bi,bj)
1087 else
1088 waqh(i,j,bi,bj) = 0.0 _d 0
1089 endif
1090
1091 if (wsfluxmm(bi,bj).ne.0.)
1092 & wsfluxmm(bi,bj) = 1./wsfluxmm(bi,bj)*wsfluxmm(bi,bj)
1093 if (whfluxmm(bi,bj).ne.0.)
1094 & whfluxmm(bi,bj) = 1./whfluxmm(bi,bj)*whfluxmm(bi,bj)
1095 cph(
1096 if (whflux2(i,j,bi,bj) .ne. 0.) then
1097 whflux2(i,j,bi,bj) =
1098 & 1./whflux2(i,j,bi,bj)/whflux2(i,j,bi,bj)
1099 else
1100 whflux2(i,j,bi,bj) = 0.0 _d 0
1101 endif
1102 if (wsflux2(i,j,bi,bj) .ne. 0.) then
1103 wsflux2(i,j,bi,bj) =
1104 & 1./wsflux2(i,j,bi,bj)/wsflux2(i,j,bi,bj)
1105 else
1106 wsflux2(i,j,bi,bj) = 0.0 _d 0
1107 endif
1108 if (wtauu2(i,j,bi,bj) .ne. 0.) then
1109 wtauu2(i,j,bi,bj) =
1110 & 1./wtauu2(i,j,bi,bj)/wtauu2(i,j,bi,bj)
1111 else
1112 wtauu2(i,j,bi,bj) = 0.0 _d 0
1113 endif
1114 if (wtauv2(i,j,bi,bj) .ne. 0.) then
1115 wtauv2(i,j,bi,bj) =
1116 & 1./wtauv2(i,j,bi,bj)/wtauv2(i,j,bi,bj)
1117 else
1118 wtauv2(i,j,bi,bj) = 0.0 _d 0
1119 endif
1120 cph)
1121 enddo
1122 enddo
1123
1124 cph(
1125 cph reduce wtheta, wsalt by factor 2.
1126 do k = 1,nr
1127 if (wtheta(k,bi,bj) .ne. 0.) then
1128 wtheta(k,bi,bj) = ratio/wtheta(k,bi,bj)/wtheta(k,bi,bj)/2.
1129 else
1130 wtheta(k,bi,bj) = 0.0 _d 0
1131 endif
1132 if (wsalt(k,bi,bj) .ne. 0.) then
1133 wsalt(k,bi,bj) = ratio/wsalt(k,bi,bj)/wsalt(k,bi,bj)/2.
1134 else
1135 wsalt(k,bi,bj) = 0.0 _d 0
1136 endif
1137 enddo
1138 cph)
1139
1140 #ifdef ALLOW_OBCS_COST_CONTRIBUTION
1141 do iobcs = 1,nobcs
1142 do k = 1,nr
1143 #ifdef ALLOW_OBCSN_COST_CONTRIBUTION
1144 if (wobcsn(k,iobcs) .ne. 0.) then
1145 wobcsn(k,iobcs) =
1146 & ratio/wobcsn(k,iobcs)/wobcsn(k,iobcs)
1147 else
1148 wobcsn(k,iobcs) = 0.0 _d 0
1149 endif
1150 #endif
1151 #ifdef ALLOW_OBCSS_COST_CONTRIBUTION
1152 if (wobcss(k,iobcs) .ne. 0.) then
1153 wobcss(k,iobcs) =
1154 & ratio/wobcss(k,iobcs)/wobcss(k,iobcs)
1155 else
1156 wobcss(k,iobcs) = 0.0 _d 0
1157 endif
1158 #endif
1159 #ifdef ALLOW_OBCSW_COST_CONTRIBUTION
1160 if (wobcsw(k,iobcs) .ne. 0.) then
1161 wobcsw(k,iobcs) =
1162 & ratio/wobcsw(k,iobcs)/wobcsw(k,iobcs)
1163 else
1164 wobcsw(k,iobcs) = 0.0 _d 0
1165 endif
1166 #endif
1167 #ifdef ALLOW_OBCSE_COST_CONTRIBUTION
1168 if (wobcse(k,iobcs) .ne. 0.) then
1169 wobcse(k,iobcs) =
1170 & ratio/wobcse(k,iobcs)/wobcse(k,iobcs)
1171 else
1172 wobcse(k,iobcs) = 0.0 _d 0
1173 endif
1174 #endif
1175 enddo
1176 enddo
1177 #endif
1178
1179 enddo
1180 enddo
1181
1182 #if (defined (ALLOW_HFLUX_COST_CONTRIBUTION))
1183 call active_write_xy_loc( 'whflux', whflux, 1, 0, mythid, dummy)
1184 call active_write_xy_loc( 'whflux2', whflux2, 1, 0, mythid, dummy)
1185 #elif (defined (ALLOW_ATEMP_COST_CONTRIBUTION))
1186 call active_write_xy_loc( 'watemp', watemp, 1, 0, mythid, dummy)
1187 #endif
1188
1189 #if (defined (ALLOW_SFLUX_COST_CONTRIBUTION))
1190 call active_write_xy_loc( 'wsflux', wsflux, 1, 0, mythid, dummy)
1191 call active_write_xy_loc( 'wsflux2', wsflux2, 1, 0, mythid, dummy)
1192 #elif (defined (ALLOW_AQH_COST_CONTRIBUTION))
1193 call active_write_xy_loc( 'waqh', waqh, 1, 0, mythid, dummy)
1194 #endif
1195
1196 #if (defined (ALLOW_USTRESS_COST_CONTRIBUTION))
1197 call active_write_xy_loc( 'wtauu', wtauu, 1, 0, mythid, dummy)
1198 call active_write_xy_loc( 'wtauu2', wtauu2, 1, 0, mythid, dummy)
1199 #elif (defined (ALLOW_UWIND_COST_CONTRIBUTION))
1200 call active_write_xy_loc( 'wuwind', wuwind, 1, 0, mythid, dummy)
1201 #endif
1202
1203 #if (defined (ALLOW_VSTRESS_COST_CONTRIBUTION))
1204 call active_write_xy_loc( 'wtauv', wtauv, 1, 0, mythid, dummy)
1205 call active_write_xy_loc( 'wtauv2', wtauv2, 1, 0, mythid, dummy)
1206 #elif (defined (ALLOW_VWIND_COST_CONTRIBUTION))
1207 call active_write_xy_loc( 'wvwind', wvwind, 1, 0, mythid, dummy)
1208 #endif
1209
1210 #ifdef ALLOW_OBCSN_COST_CONTRIBUTION
1211 #endif
1212
1213 end
  ViewVC Help
Powered by ViewVC 1.1.22