/[MITgcm]/MITgcm_contrib/arctic40km/code_ad/ecco_cost_weights.F
ViewVC logotype

Contents of /MITgcm_contrib/arctic40km/code_ad/ecco_cost_weights.F

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

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