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

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

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

Revision 1.1 - (show annotations) (download)
Fri Apr 23 19:55:12 2010 UTC (15 years, 3 months ago) by mmazloff
Branch: MAIN
CVS Tags: HEAD 
original files
1 C $Header: /u/gcmpack/MITgcm/pkg/ecco/ecco_cost_weights.F,v 1.36 2009/11/03 03:32:25 gforget 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 integer num_var
60
61 _RL factor
62 _RL wti(nr)
63 _RL wsi(nr)
64 _RL wui(nr)
65 _RL wvi(nr)
66 _RL whflux0m
67 _RL wsflux0m
68 _RL wtau0m
69 _RL ratio
70 _RL dummy
71 _RL wsshv4tmp ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
72
73 logical lwtheta2InUse
74 logical lwsalt2InUse
75 logical lwthetaLevInUse
76 logical lwsaltLevInUse
77
78 c == external ==
79
80 integer ifnblnk
81 external ifnblnk
82 integer ilnblnk
83 external ilnblnk
84
85 c == end of interface ==
86
87 lwtheta2InUse = .false.
88 lwsalt2InUse = .false.
89 lwthetaLevInUse = .false.
90 lwsaltLevInUse = .false.
91
92 jtlo = mybylo(mythid)
93 jthi = mybyhi(mythid)
94 itlo = mybxlo(mythid)
95 ithi = mybxhi(mythid)
96 jmin = 1-oly
97 jmax = sny+oly
98 imin = 1-olx
99 imax = snx+olx
100
101 c-- Initialize background weights
102 whflux0m = whflux0
103 wsflux0m = wsflux0
104 wtau0m = wtau0
105
106 c-- Initialize variance (weight) fields.
107 do k = 1,nr
108 wti(k) = 0. _d 0
109 wsi(k) = 0. _d 0
110 wui(k) = 0. _d 0
111 wvi(k) = 0. _d 0
112 enddo
113 do bj = jtlo,jthi
114 do bi = itlo,ithi
115 do j = jmin,jmax
116 do i = imin,imax
117 whflux (i,j,bi,bj) = 0. _d 0
118 whfluxm (i,j,bi,bj) = 0. _d 0
119 wsflux (i,j,bi,bj) = 0. _d 0
120 wsfluxm (i,j,bi,bj) = 0. _d 0
121 wtauu (i,j,bi,bj) = 0. _d 0
122 wtauum (i,j,bi,bj) = 0. _d 0
123 wtauv (i,j,bi,bj) = 0. _d 0
124 wtauvm (i,j,bi,bj) = 0. _d 0
125 watemp (i,j,bi,bj) = 0. _d 0
126 waqh (i,j,bi,bj) = 0. _d 0
127 wprecip (i,j,bi,bj) = 0. _d 0
128 wswflux (i,j,bi,bj) = 0. _d 0
129 wswdown (i,j,bi,bj) = 0. _d 0
130 wsnowprecip (i,j,bi,bj) = 0. _d 0
131 wlwflux (i,j,bi,bj) = 0. _d 0
132 wlwdown (i,j,bi,bj) = 0. _d 0
133 wevap (i,j,bi,bj) = 0. _d 0
134 wapressure(i,j,bi,bj) = 0. _d 0
135 wrunoff (i,j,bi,bj) = 0. _d 0
136 wuwind (i,j,bi,bj) = 0. _d 0
137 wvwind (i,j,bi,bj) = 0. _d 0
138 wsst (i,j,bi,bj) = 0. _d 0
139 wsss (i,j,bi,bj) = 0. _d 0
140 wtp (i,j,bi,bj) = 0. _d 0
141 wers (i,j,bi,bj) = 0. _d 0
142 wgfo (i,j,bi,bj) = 0. _d 0
143 do num_var=1,NSSHV4COST
144 wsshv4 (i,j,num_var,bi,bj) = 0. _d 0
145 enddo
146 wp (i,j,bi,bj) = 0. _d 0
147 wudrift (i,j,bi,bj) = 0. _d 0
148 wvdrift (i,j,bi,bj) = 0. _d 0
149 cph(
150 whflux2 (i,j,bi,bj) = 0. _d 0
151 wsflux2 (i,j,bi,bj) = 0. _d 0
152 wtauu2 (i,j,bi,bj) = 0. _d 0
153 wtauv2 (i,j,bi,bj) = 0. _d 0
154 cph)
155 wbottomdrag (i,j,bi,bj) = wbottomdrag0
156 enddo
157 enddo
158 enddo
159 enddo
160 do bj = jtlo,jthi
161 do bi = itlo,ithi
162 do k = 1,Nr
163 wtheta (k,bi,bj) = 0. _d 0
164 wsalt (k,bi,bj) = 0. _d 0
165 wctdt (k,bi,bj) = 0. _d 0
166 wctds (k,bi,bj) = 0. _d 0
167 wdiffkr(k,bi,bj) = wdiffkr0
168 wkapgm (k,bi,bj) = wkapgm0
169 wkapredi (k,bi,bj) = wkapredi0
170 wedtaux(k,bi,bj) = wedtau0
171 wedtauy(k,bi,bj) = wedtau0
172 do j = jmin,jmax
173 do i = imin,imax
174 wtheta2 (i,j,k,bi,bj) = 0. _d 0
175 wsalt2 (i,j,k,bi,bj) = 0. _d 0
176 wdiffkr2(i,j,k,bi,bj) = wdiffkr0
177 wkapgm2 (i,j,k,bi,bj) = wkapgm0
178 wkapredi2 (i,j,k,bi,bj) = wkapredi0
179 wedtaux2(i,j,k,bi,bj) = wedtau0
180 wedtauy2(i,j,k,bi,bj) = wedtau0
181 wthetaLev (i,j,k,bi,bj) = 0. _d 0
182 wsaltLev (i,j,k,bi,bj) = 0. _d 0
183 wdiffkrFld(i,j,k,bi,bj) = wdiffkr0
184 wkapgmFld (i,j,k,bi,bj) = wkapgm0
185 wkaprediFld (i,j,k,bi,bj) = wkapredi0
186 wedtauxFld(i,j,k,bi,bj) = wedtau0
187 wedtauyFld(i,j,k,bi,bj) = wedtau0
188 enddo
189 enddo
190 enddo
191 enddo
192 enddo
193
194 #if (defined (ALLOW_OBCS_COST_CONTRIBUTION) || \
195 defined (ALLOW_OBCS_CONTROL))
196 do iobcs = 1,nobcs
197 do k = 1,Nr
198 #if (defined (ALLOW_OBCSN_CONTROL) || \
199 defined (ALLOW_OBCSN_COST_CONTRIBUTION))
200 wobcsn(k,iobcs) = 0. _d 0
201 #endif
202 #if (defined (ALLOW_OBCSS_CONTROL) || \
203 defined (ALLOW_OBCSS_COST_CONTRIBUTION))
204 wobcss(k,iobcs) = 0. _d 0
205 #endif
206 #if (defined (ALLOW_OBCSW_CONTROL) || \
207 defined (ALLOW_OBCSW_COST_CONTRIBUTION))
208 wobcsw(k,iobcs) = 0. _d 0
209 #endif
210 #if (defined (ALLOW_OBCSE_CONTROL) || \
211 defined (ALLOW_OBCSE_COST_CONTRIBUTION))
212 wobcse(k,iobcs) = 0. _d 0
213 #endif
214 enddo
215 enddo
216 #endif
217
218 c-- Build area weighting matrix used in the cost function
219 c-- contributions.
220
221 c-- Define frame.
222 do j = jmin,jmax
223 do i = imin,imax
224 c-- North/South and West/East edges set to zero.
225 cph if ( (j .lt. 1) .or. (j .gt. sny) .or.
226 cph & (i .lt. 1) .or. (i .gt. snx) ) then
227 cph frame(i,j) = 0. _d 0
228 cph else
229 frame(i,j) = 1. _d 0
230 cph endif
231 enddo
232 enddo
233
234 c-- First account for the grid used.
235 if (usingCartesianGrid) then
236 factor = 0. _d 0
237 else if (usingSphericalPolarGrid) then
238 factor = 1. _d 0
239 endif
240
241 do bj = jtlo,jthi
242 do bi = itlo,ithi
243 do j = jmin,jmax
244 do i = imin,imax
245 cds cosphi(i,j,bi,bj) = cos(yc(i,j,bi,bj)*deg2rad*factor)*
246 cds & frame(i,j)
247 cosphi(i,j,bi,bj) = frame(i,j)
248 enddo
249 enddo
250 enddo
251 enddo
252
253 c-- Read error information and set up weight matrices.
254 _BEGIN_MASTER(myThid)
255 ilo = ifnblnk(data_errfile)
256 ihi = ilnblnk(data_errfile)
257 CALL OPEN_COPY_DATA_FILE(
258 I data_errfile(ilo:ihi),
259 I 'ECCO_COST_WEIGHTS',
260 O gwunit,
261 I myThid )
262
263 read(gwunit,*) ratio
264 #if (defined (ALLOW_OBCS_COST_CONTRIBUTION) || defined (ALLOW_OBCS_CONTROL))
265 & , wbaro
266 #endif
267 do k = 1,nr
268 read(gwunit,*) wti(k), wsi(k)
269 #if (defined (ALLOW_OBCS_COST_CONTRIBUTION) || defined (ALLOW_OBCS_CONTROL))
270 & , wvi(k)
271 #endif
272 end do
273 close(gwunit)
274
275 _END_MASTER(myThid)
276
277 _BARRIER
278
279 jmin = 1
280 jmax = sny
281 imin = 1
282 imax = snx
283
284 do bj = jtlo,jthi
285 do bi = itlo,ithi
286
287 wsfluxmm(bi,bj) = 1.
288 whfluxmm(bi,bj) = 1.
289
290 c-- The "classic" state estimation tool wastes memory here;
291 c-- as long as there is not more information available there
292 c-- is no need to add the zonal and meridional directions.
293 do k = 1,nr
294 wtheta(k,bi,bj) = wti(k)
295 wsalt (k,bi,bj) = wsi(k)
296 wcurrent(k,bi,bj) = wvi(k)
297 c--
298 if (wtheta(k,bi,bj) .ne. 0.) then
299 wtheta(k,bi,bj) = ratio/wtheta(k,bi,bj)/wtheta(k,bi,bj)
300 else
301 wtheta(k,bi,bj) = 0.0 _d 0
302 endif
303 if (wsalt(k,bi,bj) .ne. 0.) then
304 wsalt(k,bi,bj) = ratio/wsalt(k,bi,bj)/wsalt(k,bi,bj)
305 else
306 wsalt(k,bi,bj) = 0.0 _d 0
307 endif
308 enddo
309
310 do k = 1,nr
311 #ifdef ALLOW_OBCSN_COST_CONTRIBUTION
312 wobcsn(k,1) = wti(k)
313 wobcsn(k,2) = wsi(k)
314 wobcsn(k,3) = wvi(k)
315 wobcsn(k,4) = wvi(k)
316 #endif
317 #ifdef ALLOW_OBCSS_COST_CONTRIBUTION
318 wobcss(k,1) = wti(k)
319 wobcss(k,2) = wsi(k)
320 wobcss(k,3) = wvi(k)
321 wobcss(k,4) = wvi(k)
322 #endif
323 #ifdef ALLOW_OBCSW_COST_CONTRIBUTION
324 wobcsw(k,1) = wti(k)
325 wobcsw(k,2) = wsi(k)
326 wobcsw(k,3) = wvi(k)
327 wobcsw(k,4) = wvi(k)
328 #endif
329 #ifdef ALLOW_OBCSE_COST_CONTRIBUTION
330 wobcse(k,1) = wti(k)
331 wobcse(k,2) = wsi(k)
332 wobcse(k,3) = wvi(k)
333 wobcse(k,4) = wvi(k)
334 #endif
335 enddo
336 enddo
337 enddo
338
339 #if (defined (ALLOW_SALT0_COST_CONTRIBUTION) || defined (ALLOW_SALT0_CONTROL))
340 lwsaltLevInUse = .true.
341 if ( salt0errfile .NE. ' ' ) then
342 call mdsreadfield( salt0errfile, 32, 'RL', Nr,
343 & wsaltLev, 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 ( wsaltLev(i,j,k,bi,bj).eq.0 ) then
351 wsaltLev(i,j,k,bi,bj) = 0. _d 0
352 else
353 wsaltLev(i,j,k,bi,bj)=frame(i,j)*maskC(i,j,k,bi,bj)/
354 $ ( wsaltLev(i,j,k,bi,bj)*wsaltLev(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 wsaltLev(i,j,k,bi,bj)=
368 $ wsalt(k,bi,bj)*frame(i,j)*maskC(i,j,k,bi,bj)
369 enddo
370 enddo
371 enddo
372 enddo
373 enddo
374 endif
375 call active_write_xyz( 'wsaltLev', wsaltLev,
376 & 1, 0, mythid, dummy)
377 _EXCH_XYZ_RL( wsaltLev, myThid )
378 #endif
379
380 #if (defined (ALLOW_THETA0_COST_CONTRIBUTION) || defined (ALLOW_THETA0_CONTROL))
381 lwthetaLevInUse = .true.
382 if ( temp0errfile .NE. ' ' ) then
383 call mdsreadfield( temp0errfile, 32, 'RL', Nr,
384 & wthetaLev, 1, mythid)
385 do bj = jtlo,jthi
386 do bi = itlo,ithi
387 do k = 1,nr
388 do j = jmin,jmax
389 do i = imin,imax
390 c-- Test for missing values.
391 if ( wthetaLev(i,j,k,bi,bj).eq.0 ) then
392 wthetaLev(i,j,k,bi,bj) = 0. _d 0
393 else
394 wthetaLev(i,j,k,bi,bj)=frame(i,j)*maskC(i,j,k,bi,bj)/
395 $ ( wthetaLev(i,j,k,bi,bj)*wthetaLev(i,j,k,bi,bj) )
396 endif
397 enddo
398 enddo
399 enddo
400 enddo
401 enddo
402 else
403 do bj = jtlo,jthi
404 do bi = itlo,ithi
405 do k = 1,nr
406 do j = jmin,jmax
407 do i = imin,imax
408 wthetaLev(i,j,k,bi,bj)=
409 $ wtheta(k,bi,bj)*frame(i,j)*maskC(i,j,k,bi,bj)
410 enddo
411 enddo
412 enddo
413 enddo
414 enddo
415 endif
416 call active_write_xyz( 'wthetaLev', wthetaLev,
417 & 1, 0, mythid, dummy)
418 _EXCH_XYZ_RL( wthetaLev, myThid )
419 #endif
420
421 #if (defined (ALLOW_ARGO_SALT_COST_CONTRIBUTION) || \
422 defined (ALLOW_SSS_COST_CONTRIBUTION) || \
423 defined (ALLOW_CTDS_COST_CONTRIBUTION)|| \
424 defined (ALLOW_CTDSCLIM_COST_CONTRIBUTION))
425
426 lwsalt2InUse = .true.
427 if ( salterrfile .NE. ' ' ) then
428 call mdsreadfield( salterrfile, 32, 'RL', Nr, wsalt2, 1,
429 & mythid)
430
431 do bj = jtlo,jthi
432 do bi = itlo,ithi
433 do k = 1,nr
434 do j = jmin,jmax
435 do i = imin,imax
436 c-- Test for missing values.
437 if (wsalt(k,bi,bj).eq.0. .or.
438 $ wsalt2(i,j,k,bi,bj).eq.0.) then
439 wsalt2(i,j,k,bi,bj) = 0. _d 0
440 else
441 wsalt2(i,j,k,bi,bj)=frame(i,j)*maskC(i,j,k,bi,bj)/
442 $ ( wsalt2(i,j,k,bi,bj)*wsalt2(i,j,k,bi,bj) )
443 endif
444 enddo
445 enddo
446 enddo
447 enddo
448 enddo
449 else
450 do bj = jtlo,jthi
451 do bi = itlo,ithi
452 do k = 1,nr
453 do j = jmin,jmax
454 do i = imin,imax
455 wsalt2(i,j,k,bi,bj)=
456 $ wsalt(k,bi,bj)*frame(i,j)*maskC(i,j,k,bi,bj)
457 enddo
458 enddo
459 enddo
460 enddo
461 enddo
462 endif
463 _EXCH_XYZ_RL( wsalt2, myThid )
464 #endif
465
466 #if (defined (ALLOW_ARGO_THETA_COST_CONTRIBUTION) || \
467 defined (ALLOW_SST_COST_CONTRIBUTION) || \
468 defined (ALLOW_TMI_SST_COST_CONTRIBUTION) || \
469 defined (ALLOW_DAILYSST_COST_CONTRIBUTION) || \
470 defined (ALLOW_CTDT_COST_CONTRIBUTION) || \
471 defined (ALLOW_CTDTCLIM_COST_CONTRIBUTION) || \
472 defined (ALLOW_XBT_COST_CONTRIBUTION))
473
474 lwtheta2InUse = .true.
475 if ( temperrfile .NE. ' ' ) then
476 call mdsreadfield( temperrfile, 32, 'RL', Nr, wtheta2, 1,
477 & mythid)
478 do bj = jtlo,jthi
479 do bi = itlo,ithi
480 do k = 1,nr
481 do j = jmin,jmax
482 do i = imin,imax
483 c-- Test for missing values.
484 if (wtheta(k,bi,bj).eq.0. .or.
485 $ wtheta2(i,j,k,bi,bj).eq.0.) then
486 wtheta2(i,j,k,bi,bj) = 0. _d 0
487 else
488 wtheta2(i,j,k,bi,bj)= frame(i,j)*maskC(i,j,k,bi,bj)/
489 $ ( wtheta2(i,j,k,bi,bj)*wtheta2(i,j,k,bi,bj) )
490 endif
491 enddo
492 enddo
493 enddo
494 enddo
495 enddo
496 else
497 do bj = jtlo,jthi
498 do bi = itlo,ithi
499 do k = 1,nr
500 do j = jmin,jmax
501 do i = imin,imax
502 if (wtheta(k,bi,bj).eq.0 ) then
503 wtheta2(i,j,k,bi,bj) = 0. _d 0
504 else
505 wtheta2(i,j,k,bi,bj) =
506 $ wtheta(k,bi,bj)*frame(i,j)*maskC(i,j,k,bi,bj)
507 endif
508 enddo
509 enddo
510 enddo
511 enddo
512 enddo
513 endif
514 _EXCH_XYZ_RL( wtheta2, myThid )
515 #endif
516
517 #if (defined (ALLOW_SST_COST_CONTRIBUTION) || defined (ALLOW_SST_CONTROL))
518 if ( ssterrfile .NE. ' ' )
519 & call mdsreadfield( ssterrfile, 32, 'RL', 1, wsst, 1,
520 & mythid)
521 #endif
522 #if (defined (ALLOW_SSS_COST_CONTRIBUTION) || defined (ALLOW_SSS_CONTROL))
523 if ( ssserrfile .NE. ' ' )
524 & call mdsreadfield( ssserrfile, 32, 'RL', 1, wsss, 1,
525 & mythid)
526 #endif
527
528 k = 1
529 do bj = jtlo,jthi
530 do bi = itlo,ithi
531 do j = jmin,jmax
532 do i = imin,imax
533 #if (defined (ALLOW_SST_COST_CONTRIBUTION) || \
534 defined (ALLOW_DAILYSST_COST_CONTRIBUTION) || \
535 defined (ALLOW_SST_CONTROL))
536
537 if ( ssterrfile .NE. ' ' ) then
538 cgf use specific weights for sst
539 if (wsst(i,j,bi,bj).ne.0)
540 & wsst(i,j,bi,bj)= frame(i,j)*maskC(i,j,k,bi,bj)/
541 & ( wsst(i,j,bi,bj)*wsst(i,j,bi,bj) )
542 else
543 cgf use general hydrography weights
544 if ( lwtheta2InUse ) then
545 wsst(i,j,bi,bj) = wtheta2(i,j,k,bi,bj)
546 elseif ( lwthetaLevInUse ) then
547 wsst(i,j,bi,bj) = wthetaLev(i,j,k,bi,bj)
548 else
549 wsst(i,j,bi,bj) =
550 & wtheta(k,bi,bj)*frame(i,j)*maskC(i,j,k,bi,bj)
551 endif
552 endif
553 #endif
554
555 #if (defined (ALLOW_SSS_COST_CONTRIBUTION) || defined (ALLOW_SSS_CONTROL))
556 if ( ssserrfile .NE. ' ' ) then
557 cgf use specific weights for sss
558 if (wsss(i,j,bi,bj).ne.0)
559 & wsss(i,j,bi,bj)= frame(i,j)*maskC(i,j,k,bi,bj)/
560 & ( wsss(i,j,bi,bj)*wsss(i,j,bi,bj) )
561
562 else
563 cgf use general hydrography weights
564 if ( lwsalt2InUse ) then
565 wsss(i,j,bi,bj) = wsalt2(i,j,k,bi,bj)
566 elseif ( lwsaltLevInUse ) then
567 wsss(i,j,bi,bj) = wsaltLev(i,j,k,bi,bj)
568 else
569 wsss(i,j,bi,bj) =
570 & wsalt(k,bi,bj)*frame(i,j)*maskC(i,j,k,bi,bj)
571 endif
572 endif
573 #endif
574 enddo
575 enddo
576 enddo
577 enddo
578 CMM(came in here and changed all write_loc to write)
579 #if (defined (ALLOW_SST_COST_CONTRIBUTION) || \
580 defined (ALLOW_DAILYSST_COST_CONTRIBUTION) || \
581 defined (ALLOW_SST_CONTROL))
582 call active_write_xy( 'wsst', wsst, 1, 0, mythid, dummy)
583 #endif
584 #if (defined (ALLOW_SSS_COST_CONTRIBUTION) || defined (ALLOW_SSS_CONTROL))
585 call active_write_xy( 'wsss', wsss, 1, 0, mythid, dummy)
586 #endif
587
588 #ifdef ALLOW_EGM96_ERROR_DIAG
589 c-- Read egm-96 geoid covariance. Data in units of meters.
590 nnz = 1
591 irec = 1
592 call mdsreadfield( geoid_errfile, cost_iprec, cost_yftype,
593 & nnz, wp, irec, mythid )
594 c-- Set all tile edges to zero.
595 do bj = jtlo,jthi
596 do bi = itlo,ithi
597 do j = jmin,jmax
598 do i = imin,imax
599 wp(i,j,bi,bj) = wp(i,j,bi,bj)*frame(i,j)
600 cph-indonesian(
601 if ( xC(i,j,bi,bj) .GT. 120. .AND.
602 & xC(i,j,bi,bj) .LT. 130. .AND.
603 & yC(i,j,bi,bj) .GT. -10. .AND.
604 & yC(i,j,bi,bj) .LT. 10. ) then
605 wp(i,j,bi,bj) = wp(i,j,bi,bj)*100.
606 endif
607 cph-indonesian)
608 enddo
609 enddo
610 enddo
611 enddo
612 #else
613 do bj = jtlo,jthi
614 do bi = itlo,ithi
615 do j = jmin,jmax
616 do i = imin,imax
617 wp(i,j,bi,bj) = frame(i,j)
618 enddo
619 enddo
620 enddo
621 enddo
622 #endif
623
624 #ifdef ALLOW_SSH_COST_CONTRIBUTION
625 c-- Read SSH anomaly rms field. Data in units of centimeters.
626 nnz = 1
627 irec = 1
628 if ( ssh_errfile .NE. ' ' ) then
629 call mdsreadfield( ssh_errfile, cost_iprec, cost_yftype,
630 & nnz, wtp, irec, mythid )
631
632 do bj = jtlo,jthi
633 do bi = itlo,ithi
634 k = 1
635 do j = jmin,jmax
636 do i = imin,imax
637 c-- Unit conversion to meters. ERS error is set to
638 c-- T/P error + 5cm
639 if (maskC(i,j,k,bi,bj) .eq. 0.) then
640 wtp (i,j,bi,bj) = 0. _d 0
641 wers(i,j,bi,bj) = 0. _d 0
642 wgfo(i,j,bi,bj) = 0. _d 0
643 else
644 wtp (i,j,bi,bj) = ( wtp(i,j,bi,bj) * 0.01 * 0.5 )
645 & *frame(i,j)
646 wers(i,j,bi,bj) = ( wtp(i,j,bi,bj) + 0.05 )
647 & *frame(i,j)
648 wgfo(i,j,bi,bj) = wers(i,j,bi,bj)
649 endif
650 enddo
651 enddo
652 enddo
653 enddo
654 endif
655
656 c-- overwrite T/P error field, if available:
657 if ( tp_errfile .NE. ' ' )
658 & call mdsreadfield( tp_errfile, cost_iprec, cost_yftype, nnz,
659 & wtp, irec, mythid )
660
661 c-- overwrite ERS error field, if available:
662 if ( ers_errfile .NE. ' ' )
663 & call mdsreadfield( ers_errfile, cost_iprec, cost_yftype, nnz,
664 & wers, irec, mythid )
665
666 c-- overwrite GFO error field, if available:
667 if ( gfo_errfile .NE. ' ' )
668 & call mdsreadfield( gfo_errfile, cost_iprec, cost_yftype, nnz,
669 & wgfo, irec, mythid )
670
671 do bj = jtlo,jthi
672 do bi = itlo,ithi
673 k = 1
674 do j = jmin,jmax
675 do i = imin,imax
676 if (maskC(i,j,k,bi,bj) .eq. 0.) then
677 if ( tp_errfile .NE. ' ' )
678 & wtp (i,j,bi,bj) = 0. _d 0
679 if ( ers_errfile .NE. ' ' )
680 & wers(i,j,bi,bj) = 0. _d 0
681 if ( gfo_errfile .NE. ' ' )
682 & wgfo(i,j,bi,bj) = 0. _d 0
683 else
684 c-- convert from cm to m and set to 0.1m for missing values.
685 if ( tp_errfile .NE. ' ' ) then
686 wtp (i,j,bi,bj) = wtp (i,j,bi,bj) * 0.01 * frame(i,j)
687 cph should not be necessary for T/P and Jason
688 cph if ( wtp (i,j,bi,bj) .EQ. 0. )
689 cph & wtp (i,j,bi,bj) = 0.1 * frame(i,j)
690 endif
691 if ( ers_errfile .NE. ' ' ) then
692 wers(i,j,bi,bj) = wers(i,j,bi,bj) * 0.01 * frame(i,j)
693 if ( wers(i,j,bi,bj) .EQ. 0. )
694 & wers(i,j,bi,bj) = 0.1 * frame(i,j)
695 endif
696 if ( gfo_errfile .NE. ' ' ) then
697 wgfo(i,j,bi,bj) = wgfo(i,j,bi,bj) * 0.01 * frame(i,j)
698 if ( wgfo(i,j,bi,bj) .EQ. 0. )
699 & wgfo(i,j,bi,bj) = 0.1 * frame(i,j)
700 endif
701 endif
702 enddo
703 enddo
704 enddo
705 enddo
706
707 cph-indonesian(
708 do bj = jtlo,jthi
709 do bi = itlo,ithi
710 k = 1
711 do j = jmin,jmax
712 do i = imin,imax
713 if ( xC(i,j,bi,bj) .GT. 120. .AND.
714 & xC(i,j,bi,bj) .LT. 130. .AND.
715 & yC(i,j,bi,bj) .GT. -10. .AND.
716 & yC(i,j,bi,bj) .LT. 10. ) then
717 wtp(i,j,bi,bj) = wtp(i,j,bi,bj)*100.
718 wers(i,j,bi,bj) = wers(i,j,bi,bj)*100.
719 wgfo(i,j,bi,bj) = wgfo(i,j,bi,bj)*100.
720 endif
721 enddo
722 enddo
723 enddo
724 enddo
725 cph-indonesian)
726 CMM -- polar oceans -- SLA is questionable
727 CMM TP turning angle is -66.2 but GRACE ends at 65.8 so use
728 CMM is index 1 to 73 -- doesnt get all points -- need to address
729 do bj = jtlo,jthi
730 do bi = itlo,ithi
731 k = 1
732 do j = jmin,jmax
733 do i = imin,imax
734 if ( yC(i,j,bi,bj) .LT. -65.8 ) then
735 wp(i,j,bi,bj) = 0.
736 wtp(i,j,bi,bj) = 0.
737 wers(i,j,bi,bj) = 0.
738 wgfo(i,j,bi,bj) = 0.
739 endif
740 enddo
741 enddo
742 enddo
743 enddo
744 CMM)
745
746
747 #endif /* ALLOW_SSH_COST_CONTRIBUTION */
748
749 #ifdef ALLOW_SSHV4_COST
750
751 do num_var=1,NSSHV4COST
752
753 if ( sshv4cost_errfile(num_var) .NE. ' ' ) then
754 c-- read error standard deviation
755 call mdsreadfield( sshv4cost_errfile(num_var),
756 & 32, 'RL', 1, wsshv4tmp, 1, mythid)
757 c-- convert to units of meters
758 do bj = jtlo,jthi
759 do bi = itlo,ithi
760 do j = jmin,jmax
761 do i = imin,imax
762 wsshv4tmp(i,j,bi,bj)=wsshv4tmp(i,j,bi,bj)
763 & *sshv4cost_errfactor(num_var)
764 enddo
765 enddo
766 enddo
767 enddo
768 else
769 do bj = jtlo,jthi
770 do bi = itlo,ithi
771 do j = jmin,jmax
772 do i = imin,imax
773 wsshv4tmp(i,j,bi,bj)=0. _d 0
774 enddo
775 enddo
776 enddo
777 enddo
778 endif
779
780 do bj = jtlo,jthi
781 do bi = itlo,ithi
782 do j = jmin,jmax
783 do i = imin,imax
784 if (wsshv4tmp(i,j,bi,bj).ne.0) then
785 wsshv4tmp(i,j,bi,bj)= frame(i,j)*maskC(i,j,k,bi,bj)/
786 & ( wsshv4tmp(i,j,bi,bj)* wsshv4tmp(i,j,bi,bj) )
787 wsshv4(i,j,num_var,bi,bj)=wsshv4tmp(i,j,bi,bj)
788 endif
789 enddo
790 enddo
791 enddo
792 enddo
793
794 call active_write_xy( 'wsshv4', wsshv4tmp,
795 & num_var, 0, mythid, dummy)
796
797 enddo
798
799 #endif
800
801 #ifdef ALLOW_BP_COST_CONTRIBUTION
802 if ( bperrfile .NE. ' ' )
803 & call mdsreadfield( bperrfile, 32, 'RL', 1, wbp, 1,
804 & mythid)
805
806 do bj = jtlo,jthi
807 do bi = itlo,ithi
808 do j = jmin,jmax
809 do i = imin,imax
810 if (wbp(i,j,bi,bj).ne.0)
811 & wbp(i,j,bi,bj)= frame(i,j)*maskC(i,j,k,bi,bj)/
812 & ( wbp(i,j,bi,bj)* wbp(i,j,bi,bj) )
813 enddo
814 enddo
815 enddo
816 enddo
817
818 call active_write_xy( 'wbp', wbp, 1, 0, mythid, dummy)
819 #endif
820
821 c-- Read zonal wind stress variance.
822 #if (defined (ALLOW_SCAT_COST_CONTRIBUTION))
823
824 nnz = 1
825 irec = 1
826 call mdsreadfield( scatx_errfile, cost_iprec, cost_yftype, nnz,
827 & wscatx, irec, mythid )
828 call mdsreadfield( scaty_errfile, cost_iprec, cost_yftype, nnz,
829 & wscaty, irec, mythid )
830
831 do bj = jtlo,jthi
832 do bi = itlo,ithi
833 k = 1
834 do j = jmin,jmax
835 do i = imin,imax
836 c-- Test for missing values.
837 if (wscatx(i,j,bi,bj) .lt. -9900.) then
838 wscatx(i,j,bi,bj) = 0. _d 0
839 endif
840 wscatx(i,j,bi,bj) = wscatx(i,j,bi,bj)
841 wscatx(i,j,bi,bj) = max(wscatx(i,j,bi,bj),wtau0)
842 wscatx(i,j,bi,bj) = wscatx(i,j,bi,bj)*maskw(i,j,k,bi,bj)
843 & *frame(i,j)
844 if (wscaty(i,j,bi,bj) .lt. -9900.) then
845 wscaty(i,j,bi,bj) = 0. _d 0
846 endif
847 wscaty(i,j,bi,bj) = wscaty(i,j,bi,bj)
848 wscaty(i,j,bi,bj) = max(wscaty(i,j,bi,bj),wtau0)
849 wscaty(i,j,bi,bj) = wscaty(i,j,bi,bj)*masks(i,j,k,bi,bj)
850 & *frame(i,j)
851 enddo
852 enddo
853 enddo
854 enddo
855
856 #endif
857
858 c-- Read zonal wind stress variance.
859 #if (defined (ALLOW_STRESS_MEAN_COST_CONTRIBUTION))
860 nnz = 1
861 irec = 1
862 cph call mdsreadfield( tauum_errfile, cost_iprec, cost_yftype,
863 cph & nnz, wtauum, irec, mythid )
864 cph call mdsreadfield( tauvm_errfile, cost_iprec, cost_yftype,
865 cph & nnz, wtauvm, irec, mythid )
866
867 do bj = jtlo,jthi
868 do bi = itlo,ithi
869 k = 1
870 do j = jmin,jmax
871 do i = imin,imax
872 c-- Test for missing values.
873 if (wtauum(i,j,bi,bj) .lt. -9900.) then
874 wtauum(i,j,bi,bj) = 0. _d 0
875 endif
876 wtauum(i,j,bi,bj) = max(wtauum(i,j,bi,bj),wtau0m)
877 & *frame(i,j)
878 c-- Test for missing values.
879 if (wtauvm(i,j,bi,bj) .lt. -9900.) then
880 wtauvm(i,j,bi,bj) = 0. _d 0
881 endif
882 wtauvm(i,j,bi,bj) = max(wtauvm(i,j,bi,bj),wtau0m)
883 & *frame(i,j)
884 enddo
885 enddo
886 enddo
887 enddo
888 #endif
889
890 #if (defined (ALLOW_USTRESS_COST_CONTRIBUTION) || defined (ALLOW_USTRESS_CONTROL))
891 nnz = 1
892 ce irec = 2
893 ce( due to Patrick's processing:
894 irec = 1
895 ce)
896 if ( tauu_errfile .NE. ' ' ) then
897 call mdsreadfield( tauu_errfile, cost_iprec, cost_yftype,
898 & nnz, wtauu, irec, mythid )
899 endif
900
901 do bj = jtlo,jthi
902 do bi = itlo,ithi
903 k = 1
904 do j = jmin,jmax
905 do i = imin,imax
906 c-- Test for missing values.
907 if (wtauu(i,j,bi,bj) .lt. -9900.) then
908 wtauu(i,j,bi,bj) = 0. _d 0
909 endif
910 wtauu(i,j,bi,bj) = wtauu(i,j,bi,bj)*0.1
911 wtauu(i,j,bi,bj) = max(wtauu(i,j,bi,bj),wtau0)
912 wtauu(i,j,bi,bj) = wtauu(i,j,bi,bj)*maskw(i,j,k,bi,bj)
913 & *frame(i,j)
914 cph(
915 wtauu2(i,j,bi,bj) = wtau0*maskW(i,j,k,bi,bj)*frame(i,j)
916 cph)
917 enddo
918 enddo
919 enddo
920 enddo
921
922 #elif (defined (ALLOW_UWIND_COST_CONTRIBUTION) || defined (ALLOW_UWIND_CONTROL))
923
924 nnz = 1
925 ce irec = 2
926 ce( due to Patrick's processing:
927 irec = 1
928 ce)
929 if ( uwind_errfile .NE. ' ' ) then
930 call mdsreadfield( uwind_errfile, cost_iprec, cost_yftype,
931 & nnz, wuwind, irec, mythid )
932 endif
933
934 do bj = jtlo,jthi
935 do bi = itlo,ithi
936 k = 1
937 do j = jmin,jmax
938 do i = imin,imax
939 c-- Test for missing values.
940 if (wuwind(i,j,bi,bj) .lt. -9900.) then
941 wuwind(i,j,bi,bj) = 0. _d 0
942 endif
943 wuwind(i,j,bi,bj) = wuwind(i,j,bi,bj)
944 wuwind(i,j,bi,bj) = max(wuwind(i,j,bi,bj),wwind0)
945 wuwind(i,j,bi,bj) = wuwind(i,j,bi,bj)*maskc(i,j,k,bi,bj)
946 & *frame(i,j)
947 enddo
948 enddo
949 enddo
950 enddo
951 #endif
952
953 c-- Read meridional wind stress variance.
954 #if (defined (ALLOW_VSTRESS_COST_CONTRIBUTION) || defined (ALLOW_VSTRESS_CONTROL))
955 nnz = 1
956 ce irec = 2
957 ce( due to Patrick's processing:
958 irec = 1
959 ce)
960
961 if ( tauv_errfile .NE. ' ' ) then
962 call mdsreadfield( tauv_errfile, cost_iprec, cost_yftype, nnz,
963 & wtauv, irec, mythid )
964 endif
965
966 do bj = jtlo,jthi
967 do bi = itlo,ithi
968 do j = jmin,jmax
969 do i = imin,imax
970 c-- Test for missing values.
971 if (wtauv(i,j,bi,bj) .lt. -9900.) then
972 wtauv(i,j,bi,bj) = 0. _d 0
973 endif
974 wtauv(i,j,bi,bj) = wtauv(i,j,bi,bj)*0.1
975 wtauv(i,j,bi,bj) = max(wtauv(i,j,bi,bj),wtau0)
976 wtauv(i,j,bi,bj) = wtauv(i,j,bi,bj)*masks(i,j,k,bi,bj)
977 & *frame(i,j)
978 cph(
979 wtauv2(i,j,bi,bj) = wtau0*maskS(i,j,k,bi,bj)*frame(i,j)
980 cph)
981 enddo
982 enddo
983 enddo
984 enddo
985
986 #elif (defined (ALLOW_VWIND_COST_CONTRIBUTION) || defined (ALLOW_VWIND_CONTROL))
987
988 nnz = 1
989 ce irec = 2
990 ce( due to Patrick's processing:
991 irec = 1
992 ce)
993
994 if ( vwind_errfile .NE. ' ' ) then
995 call mdsreadfield( vwind_errfile, cost_iprec, cost_yftype,
996 & nnz, wvwind, irec, mythid )
997 endif
998
999 do bj = jtlo,jthi
1000 do bi = itlo,ithi
1001 do j = jmin,jmax
1002 do i = imin,imax
1003 c-- Test for missing values.
1004 if (wvwind(i,j,bi,bj) .lt. -9900.) then
1005 wvwind(i,j,bi,bj) = 0. _d 0
1006 endif
1007 wvwind(i,j,bi,bj) = wvwind(i,j,bi,bj)
1008 wvwind(i,j,bi,bj) = max(wvwind(i,j,bi,bj),wwind0)
1009 wvwind(i,j,bi,bj) = wvwind(i,j,bi,bj)*maskc(i,j,k,bi,bj)
1010 & *frame(i,j)
1011 enddo
1012 enddo
1013 enddo
1014 enddo
1015 #endif
1016
1017 #if (defined (ALLOW_HFLUX_COST_CONTRIBUTION) || defined (ALLOW_HFLUX_CONTROL))
1018 c-- Read heat flux flux variance.
1019 nnz = 1
1020 c-- First record in data file: mean field.
1021 c-- Second record in data file: rms field.
1022 ce irec = 2
1023 ce( due to Patrick's processing:
1024 irec = 1
1025 ce)
1026 if ( hflux_errfile .NE. ' ' ) then
1027 call mdsreadfield( hflux_errfile, cost_iprec, cost_yftype,
1028 & nnz, whflux, irec, mythid )
1029 endif
1030
1031 do bj = jtlo,jthi
1032 do bi = itlo,ithi
1033 do j = jmin,jmax
1034 do i = imin,imax
1035 c-- Test for missing values.
1036 if (whflux(i,j,bi,bj) .lt. -9900.) then
1037 whflux(i,j,bi,bj) = 0. _d 0
1038 endif
1039 c-- Data are in units of W/m**2.
1040 whflux(i,j,bi,bj) = whflux(i,j,bi,bj)/3.
1041 whflux(i,j,bi,bj) = max(whflux(i,j,bi,bj),whflux0)
1042 & *frame(i,j)
1043 whfluxm(i,j,bi,bj) = max(whfluxm(i,j,bi,bj),whflux0m)
1044 & *frame(i,j)
1045 cph(
1046 whflux2(i,j,bi,bj) = whflux0*frame(i,j)
1047 cph)
1048 enddo
1049 enddo
1050 enddo
1051 enddo
1052 #elif (defined (ALLOW_ATEMP_COST_CONTRIBUTION) || defined (ALLOW_ATEMP_CONTROL))
1053 c-- Read atmos. temp. variance.
1054 nnz = 1
1055 c-- First record in data file: mean field.
1056 c-- Second record in data file: rms field.
1057 ce irec = 2
1058 ce( due to Patrick's processing:
1059 irec = 1
1060 ce)
1061 if ( atemp_errfile .NE. ' ' ) then
1062 call mdsreadfield( atemp_errfile, cost_iprec, cost_yftype,
1063 & nnz, watemp, irec, mythid )
1064 endif
1065
1066 do bj = jtlo,jthi
1067 do bi = itlo,ithi
1068 do j = jmin,jmax
1069 do i = imin,imax
1070 c-- Test for missing values.
1071 if (watemp(i,j,bi,bj) .lt. -9900.) then
1072 watemp(i,j,bi,bj) = 0. _d 0
1073 endif
1074 c-- Data are in units of W/m**2.
1075 watemp(i,j,bi,bj) = watemp(i,j,bi,bj)
1076 watemp(i,j,bi,bj) = max(watemp(i,j,bi,bj),watemp0)
1077 & *frame(i,j)
1078 enddo
1079 enddo
1080 enddo
1081 enddo
1082
1083
1084 #endif
1085
1086 #if (defined (ALLOW_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SFLUX_CONTROL))
1087 c-- Read salt flux variance. Second read: data in units of m/s.
1088 nnz = 1
1089 c-- First record in data file: mean field.
1090 c-- Second record in data file: rms field.
1091 ce irec = 2
1092 ce( due to Patrick's processing:
1093 irec = 1
1094 ce)
1095 if ( sflux_errfile .NE. ' ' ) then
1096 call mdsreadfield( sflux_errfile, cost_iprec, cost_yftype,
1097 & nnz, wsflux, irec, mythid )
1098 endif
1099
1100 do bj = jtlo,jthi
1101 do bi = itlo,ithi
1102 do j = jmin,jmax
1103 do i = imin,imax
1104 c-- Test for missing values.
1105 if (wsflux(i,j,bi,bj) .lt. -9900.) then
1106 wsflux(i,j,bi,bj) = 0. _d 0
1107 endif
1108 c-- Data are in units of m/s.
1109 wsflux(i,j,bi,bj) = wsflux(i,j,bi,bj) / 3.
1110 wsflux(i,j,bi,bj) = max(wsflux(i,j,bi,bj),wsflux0)
1111 & *frame(i,j)
1112 wsfluxm(i,j,bi,bj) = max(wsfluxm(i,j,bi,bj),wsflux0m)
1113 & *frame(i,j)
1114 cph(
1115 wsflux2(i,j,bi,bj) = wsflux0*frame(i,j)
1116 cph)
1117 enddo
1118 enddo
1119 enddo
1120 enddo
1121 #elif (defined (ALLOW_AQH_COST_CONTRIBUTION) || defined (ALLOW_AQH_CONTROL))
1122 c-- Secific humid. variance. Second read: data in units of m/s.
1123 nnz = 1
1124 c-- First record in data file: mean field.
1125 c-- Second record in data file: rms field.
1126 ce irec = 2
1127 ce( due to Patrick's processing:
1128 irec = 1
1129 ce)
1130 if ( aqh_errfile .NE. ' ' ) then
1131 call mdsreadfield( aqh_errfile, cost_iprec, cost_yftype, nnz,
1132 & waqh, irec, mythid )
1133 endif
1134
1135 do bj = jtlo,jthi
1136 do bi = itlo,ithi
1137 do j = jmin,jmax
1138 do i = imin,imax
1139 c-- Test for missing values.
1140 if (waqh(i,j,bi,bj) .lt. -9900.) then
1141 waqh(i,j,bi,bj) = 0. _d 0
1142 endif
1143 c-- Data are in units of m/s.
1144 waqh(i,j,bi,bj) = waqh(i,j,bi,bj)
1145 waqh(i,j,bi,bj) = max(waqh(i,j,bi,bj),waqh0)
1146 & *frame(i,j)
1147 enddo
1148 enddo
1149 enddo
1150 enddo
1151 #endif
1152
1153 #if (defined (ALLOW_PRECIP_COST_CONTRIBUTION) || defined (ALLOW_PRECIP_CONTROL))
1154 c-- Atmos. precipitation variance. Second read: data in units of m/s.
1155 nnz = 1
1156 c-- First record in data file: mean field.
1157 c-- Second record in data file: rms field.
1158 ce irec = 2
1159 ce( due to Patrick's processing:
1160 irec = 1
1161 ce)
1162 if ( precip_errfile .NE. ' ' ) then
1163 call mdsreadfield( precip_errfile, cost_iprec, cost_yftype,
1164 & nnz, wprecip, irec, mythid )
1165 endif
1166
1167 do bj = jtlo,jthi
1168 do bi = itlo,ithi
1169 do j = jmin,jmax
1170 do i = imin,imax
1171 c-- Test for missing values.
1172 if (wprecip(i,j,bi,bj) .lt. -9900.) then
1173 wprecip(i,j,bi,bj) = 0. _d 0
1174 endif
1175 c-- Data are in units of m/s.
1176 wprecip(i,j,bi,bj) = wprecip(i,j,bi,bj)
1177 wprecip(i,j,bi,bj) = max(wprecip(i,j,bi,bj),wprecip0)
1178 & *frame(i,j)
1179 enddo
1180 enddo
1181 enddo
1182 enddo
1183 #endif
1184
1185 #if (defined (ALLOW_SWFLUX_COST_CONTRIBUTION) || defined (ALLOW_SWFLUX_CONTROL))
1186 c-- Atmos. swfluxitation variance. Second read: data in units of m/s.
1187 nnz = 1
1188 c-- First record in data file: mean field.
1189 c-- Second record in data file: rms field.
1190 ce irec = 2
1191 ce( due to Patrick's processing:
1192 irec = 1
1193 ce)
1194 if ( swflux_errfile .NE. ' ' ) then
1195 call mdsreadfield( swflux_errfile, cost_iprec, cost_yftype,
1196 & nnz, wswflux, irec, mythid )
1197 endif
1198
1199 do bj = jtlo,jthi
1200 do bi = itlo,ithi
1201 do j = jmin,jmax
1202 do i = imin,imax
1203 c-- Test for missing values.
1204 if (wswflux(i,j,bi,bj) .lt. -9900.) then
1205 wswflux(i,j,bi,bj) = 0. _d 0
1206 endif
1207 c-- Data are in units of m/s.
1208 wswflux(i,j,bi,bj) = wswflux(i,j,bi,bj)
1209 wswflux(i,j,bi,bj) = max(wswflux(i,j,bi,bj),wswflux0)
1210 & *frame(i,j)
1211 enddo
1212 enddo
1213 enddo
1214 enddo
1215 #endif
1216
1217 #if (defined (ALLOW_SWDOWN_COST_CONTRIBUTION) || defined (ALLOW_SWDOWN_CONTROL))
1218 c-- Atmos. swdownitation variance. Second read: data in units of m/s.
1219 nnz = 1
1220 c-- First record in data file: mean field.
1221 c-- Second record in data file: rms field.
1222 ce irec = 2
1223 ce( due to Patrick's processing:
1224 irec = 1
1225 ce)
1226 if ( swdown_errfile .NE. ' ' ) then
1227 call mdsreadfield( swdown_errfile, cost_iprec, cost_yftype,
1228 & nnz, wswdown, irec, mythid )
1229 endif
1230
1231 do bj = jtlo,jthi
1232 do bi = itlo,ithi
1233 do j = jmin,jmax
1234 do i = imin,imax
1235 c-- Test for missing values.
1236 if (wswdown(i,j,bi,bj) .lt. -9900.) then
1237 wswdown(i,j,bi,bj) = 0. _d 0
1238 endif
1239 c-- Data are in units of m/s.
1240 wswdown(i,j,bi,bj) = wswdown(i,j,bi,bj)
1241 wswdown(i,j,bi,bj) = max(wswdown(i,j,bi,bj),wswdown0)
1242 & *frame(i,j)
1243 enddo
1244 enddo
1245 enddo
1246 enddo
1247 #endif
1248
1249 #if (defined (ALLOW_LWFLUX_COST_CONTRIBUTION) || defined (ALLOW_LWFLUX_CONTROL))
1250 c-- Atmos. lwfluxitation variance. Second read: data in units of m/s.
1251 nnz = 1
1252 c-- First record in data file: mean field.
1253 c-- Second record in data file: rms field.
1254 ce irec = 2
1255 ce( due to Patrick's processing:
1256 irec = 1
1257 ce)
1258 if ( lwflux_errfile .NE. ' ' ) then
1259 call mdsreadfield( lwflux_errfile, cost_iprec, cost_yftype,
1260 & nnz, wlwflux, irec, mythid )
1261 endif
1262
1263 do bj = jtlo,jthi
1264 do bi = itlo,ithi
1265 do j = jmin,jmax
1266 do i = imin,imax
1267 c-- Test for missing values.
1268 if (wlwflux(i,j,bi,bj) .lt. -9900.) then
1269 wlwflux(i,j,bi,bj) = 0. _d 0
1270 endif
1271 c-- Data are in units of m/s.
1272 wlwflux(i,j,bi,bj) = wlwflux(i,j,bi,bj)
1273 wlwflux(i,j,bi,bj) = max(wlwflux(i,j,bi,bj),wlwflux0)
1274 & *frame(i,j)
1275 enddo
1276 enddo
1277 enddo
1278 enddo
1279 #endif
1280
1281 #if (defined (ALLOW_LWDOWN_COST_CONTRIBUTION) || defined (ALLOW_LWDOWN_CONTROL))
1282 c-- Atmos. lwdownitation variance. Second read: data in units of m/s.
1283 nnz = 1
1284 c-- First record in data file: mean field.
1285 c-- Second record in data file: rms field.
1286 ce irec = 2
1287 ce( due to Patrick's processing:
1288 irec = 1
1289 ce)
1290 if ( lwdown_errfile .NE. ' ' ) then
1291 call mdsreadfield( lwdown_errfile, cost_iprec, cost_yftype,
1292 & nnz, wlwdown, irec, mythid )
1293 endif
1294
1295 do bj = jtlo,jthi
1296 do bi = itlo,ithi
1297 do j = jmin,jmax
1298 do i = imin,imax
1299 c-- Test for missing values.
1300 if (wlwdown(i,j,bi,bj) .lt. -9900.) then
1301 wlwdown(i,j,bi,bj) = 0. _d 0
1302 endif
1303 c-- Data are in units of m/s.
1304 wlwdown(i,j,bi,bj) = wlwdown(i,j,bi,bj)
1305 wlwdown(i,j,bi,bj) = max(wlwdown(i,j,bi,bj),wlwdown0)
1306 & *frame(i,j)
1307 enddo
1308 enddo
1309 enddo
1310 enddo
1311 #endif
1312
1313 #if (defined (ALLOW_SNOWPRECIP_COST_CONTRIBUTION) || defined (ALLOW_SNOWPRECIP_CONTROL))
1314 c-- Atmos. snowprecipitation variance. Second read: data in units of m/s.
1315 nnz = 1
1316 c-- First record in data file: mean field.
1317 c-- Second record in data file: rms field.
1318 ce irec = 2
1319 ce( due to Patrick's processing:
1320 irec = 1
1321 ce)
1322 if ( snowprecip_errfile .NE. ' ' ) then
1323 call mdsreadfield( snowprecip_errfile, cost_iprec, cost_yftype,
1324 & nnz, wsnowprecip, irec, mythid )
1325 endif
1326
1327 do bj = jtlo,jthi
1328 do bi = itlo,ithi
1329 do j = jmin,jmax
1330 do i = imin,imax
1331 c-- Test for missing values.
1332 if (wsnowprecip(i,j,bi,bj) .lt. -9900.) then
1333 wsnowprecip(i,j,bi,bj) = 0. _d 0
1334 endif
1335 c-- Data are in units of m/s.
1336 wsnowprecip(i,j,bi,bj) = wsnowprecip(i,j,bi,bj)
1337 wsnowprecip(i,j,bi,bj) =
1338 & max(wsnowprecip(i,j,bi,bj),wsnowprecip0)
1339 & *frame(i,j)
1340 enddo
1341 enddo
1342 enddo
1343 enddo
1344 #endif
1345
1346 #if (defined (ALLOW_EVAP_COST_CONTRIBUTION) || defined (ALLOW_EVAP_CONTROL))
1347 c-- Atmos. evapitation variance. Second read: data in units of m/s.
1348 nnz = 1
1349 c-- First record in data file: mean field.
1350 c-- Second record in data file: rms field.
1351 ce irec = 2
1352 ce( due to Patrick's processing:
1353 irec = 1
1354 ce)
1355 if ( evap_errfile .NE. ' ' ) then
1356 call mdsreadfield( evap_errfile, cost_iprec, cost_yftype,
1357 & nnz, wevap, irec, mythid )
1358 endif
1359
1360 do bj = jtlo,jthi
1361 do bi = itlo,ithi
1362 do j = jmin,jmax
1363 do i = imin,imax
1364 c-- Test for missing values.
1365 if (wevap(i,j,bi,bj) .lt. -9900.) then
1366 wevap(i,j,bi,bj) = 0. _d 0
1367 endif
1368 c-- Data are in units of m/s.
1369 wevap(i,j,bi,bj) = wevap(i,j,bi,bj)
1370 wevap(i,j,bi,bj) = max(wevap(i,j,bi,bj),wevap0)
1371 & *frame(i,j)
1372 enddo
1373 enddo
1374 enddo
1375 enddo
1376 #endif
1377
1378 #if (defined (ALLOW_APRESSURE_COST_CONTRIBUTION) || defined (ALLOW_APRESSURE_CONTROL))
1379 c-- Atmos. apressureitation variance. Second read: data in units of m/s.
1380 nnz = 1
1381 c-- First record in data file: mean field.
1382 c-- Second record in data file: rms field.
1383 ce irec = 2
1384 ce( due to Patrick's processing:
1385 irec = 1
1386 ce)
1387 if ( apressure_errfile .NE. ' ' ) then
1388 call mdsreadfield( apressure_errfile, cost_iprec, cost_yftype,
1389 & nnz, wapressure, irec, mythid )
1390 endif
1391
1392 do bj = jtlo,jthi
1393 do bi = itlo,ithi
1394 do j = jmin,jmax
1395 do i = imin,imax
1396 c-- Test for missing values.
1397 if (wapressure(i,j,bi,bj) .lt. -9900.) then
1398 wapressure(i,j,bi,bj) = 0. _d 0
1399 endif
1400 c-- Data are in units of m/s.
1401 wapressure(i,j,bi,bj) = wapressure(i,j,bi,bj)
1402 wapressure(i,j,bi,bj) =
1403 & max(wapressure(i,j,bi,bj),wapressure0)
1404 & *frame(i,j)
1405 enddo
1406 enddo
1407 enddo
1408 enddo
1409 #endif
1410
1411 #if (defined (ALLOW_RUNOFF_COST_CONTRIBUTION) || defined (ALLOW_RUNOFF_CONTROL))
1412 c-- Atmos. runoffitation variance. Second read: data in units of m/s.
1413 nnz = 1
1414 c-- First record in data file: mean field.
1415 c-- Second record in data file: rms field.
1416 ce irec = 2
1417 ce( due to Patrick's processing:
1418 irec = 1
1419 ce)
1420 if ( runoff_errfile .NE. ' ' ) then
1421 call mdsreadfield( runoff_errfile, cost_iprec, cost_yftype,
1422 & nnz, wrunoff, irec, mythid )
1423 endif
1424
1425 do bj = jtlo,jthi
1426 do bi = itlo,ithi
1427 do j = jmin,jmax
1428 do i = imin,imax
1429 c-- Test for missing values.
1430 if (wrunoff(i,j,bi,bj) .lt. -9900.) then
1431 wrunoff(i,j,bi,bj) = 0. _d 0
1432 endif
1433 c-- Data are in units of m/s.
1434 wrunoff(i,j,bi,bj) = wrunoff(i,j,bi,bj)
1435 wrunoff(i,j,bi,bj) = max(wrunoff(i,j,bi,bj),wrunoff0)
1436 & *frame(i,j)
1437 enddo
1438 enddo
1439 enddo
1440 enddo
1441 #endif
1442
1443 #if (defined (ALLOW_BOTTOMDRAG_COST_CONTRIBUTION) || defined (ALLOW_BOTTOMDRAG_CONTROL))
1444 if ( bottomdrag_errfile .NE. ' ' ) then
1445
1446 call mdsreadfield( bottomdrag_errfile, cost_iprec, cost_yftype,
1447 & nnz, wbottomdrag, irec, mythid )
1448
1449 do bj = jtlo,jthi
1450 do bi = itlo,ithi
1451 do j = jmin,jmax
1452 do i = imin,imax
1453 c-- Test for missing values.
1454 if (wbottomdrag(i,j,bi,bj) .lt. -9900.) then
1455 wbottomdrag(i,j,bi,bj) = 0. _d 0
1456 endif
1457 enddo
1458 enddo
1459 enddo
1460 enddo
1461
1462 endif
1463 #endif
1464
1465 #if (defined (ALLOW_DIFFKR_COST_CONTRIBUTION) || defined (ALLOW_DIFFKR_CONTROL))
1466 if ( diffkr_errfile .NE. ' ' ) then
1467
1468 call mdsreadfield( diffkr_errfile, cost_iprec, cost_yftype,
1469 & Nr, wdiffkr2, 1, mythid )
1470
1471 do bj = jtlo,jthi
1472 do bi = itlo,ithi
1473 do k = 1,nr
1474 do j = jmin,jmax
1475 do i = imin,imax
1476 c-- Test for missing values.
1477 if (wdiffkr2(i,j,k,bi,bj) .lt. -9900.) then
1478 wdiffkr2(i,j,k,bi,bj) = 0. _d 0
1479 endif
1480 enddo
1481 enddo
1482 enddo
1483 enddo
1484 enddo
1485
1486 endif
1487 #endif
1488
1489 #if (defined (ALLOW_KAPGM_COST_CONTRIBUTION) || defined (ALLOW_KAPGM_CONTROL))
1490 if ( kapgm_errfile .NE. ' ' ) then
1491
1492 call mdsreadfield( kapgm_errfile, cost_iprec, cost_yftype,
1493 & Nr, wkapgm2, 1, mythid )
1494
1495 do bj = jtlo,jthi
1496 do bi = itlo,ithi
1497 do k = 1,nr
1498 do j = jmin,jmax
1499 do i = imin,imax
1500 c-- Test for missing values.
1501 if (wkapgm2(i,j,k,bi,bj) .lt. -9900.) then
1502 wkapgm2(i,j,k,bi,bj) = 0. _d 0
1503 endif
1504 enddo
1505 enddo
1506 enddo
1507 enddo
1508 enddo
1509
1510 endif
1511 #endif
1512
1513 #if (defined (ALLOW_KAPREDI_COST_CONTRIBUTION) || defined (ALLOW_KAPREDI_CONTROL))
1514 if ( kapredi_errfile .NE. ' ' ) then
1515
1516 call mdsreadfield( kapredi_errfile, cost_iprec, cost_yftype,
1517 & Nr, wkapredi2, 1, mythid )
1518
1519 do bj = jtlo,jthi
1520 do bi = itlo,ithi
1521 do k = 1,nr
1522 do j = jmin,jmax
1523 do i = imin,imax
1524 c-- Test for missing values.
1525 if (wkapredi2(i,j,k,bi,bj) .lt. -9900.) then
1526 wkapredi2(i,j,k,bi,bj) = 0. _d 0
1527 endif
1528 enddo
1529 enddo
1530 enddo
1531 enddo
1532 enddo
1533
1534 endif
1535 #endif
1536
1537 #if ( defined (ALLOW_EDDYPSI_COST_CONTRIBUTION) || defined (ALLOW_EDDYPSI_CONTROL) )
1538 if ( edtau_errfile .NE. ' ' ) then
1539
1540 call mdsreadfield( edtau_errfile, cost_iprec, cost_yftype,
1541 & Nr, wedtaux2, 1, mythid )
1542
1543 do bj = jtlo,jthi
1544 do bi = itlo,ithi
1545 do k = 1,nr
1546 do j = jmin,jmax
1547 do i = imin,imax
1548 c-- Test for missing values.
1549 if (wedtaux2(i,j,k,bi,bj) .lt. -9900.) then
1550 wedtaux2(i,j,k,bi,bj) = 0. _d 0
1551 endif
1552 wedtauy2(i,j,k,bi,bj)=wedtaux2(i,j,k,bi,bj)
1553 enddo
1554 enddo
1555 enddo
1556 enddo
1557 enddo
1558
1559 endif
1560 #endif
1561
1562 c-- Units have to be checked!
1563 do bj = jtlo,jthi
1564 do bi = itlo,ithi
1565 do j = jmin,jmax
1566 do i = imin,imax
1567 if (wtp(i,j,bi,bj) .ne. 0.) then
1568 wtp (i,j,bi,bj) = 1./wtp(i,j,bi,bj)/wtp(i,j,bi,bj)
1569 endif
1570 if (wers(i,j,bi,bj) .ne. 0.) then
1571 wers(i,j,bi,bj) = 1./wers(i,j,bi,bj)/wers(i,j,bi,bj)
1572 endif
1573 if (wgfo(i,j,bi,bj) .ne. 0.) then
1574 wgfo(i,j,bi,bj) = 1./wgfo(i,j,bi,bj)/wgfo(i,j,bi,bj)
1575 endif
1576 if (wp(i,j,bi,bj) .ne. 0.) then
1577 wp(i,j,bi,bj) = 1./wp(i,j,bi,bj)/wp(i,j,bi,bj)
1578 endif
1579 if (wtauu(i,j,bi,bj) .ne. 0.) then
1580 wtauu(i,j,bi,bj) = 1./wtauu(i,j,bi,bj)/wtauu(i,j,bi,bj)
1581 else
1582 wtauu(i,j,bi,bj) = 0.0 _d 0
1583 endif
1584 if (wtauum(i,j,bi,bj) .ne. 0.) then
1585 wtauum(i,j,bi,bj) =
1586 & 1./wtauum(i,j,bi,bj)/wtauum(i,j,bi,bj)
1587 else
1588 wtauum(i,j,bi,bj) = 0.0 _d 0
1589 endif
1590 if (wscatx(i,j,bi,bj) .ne. 0.) then
1591 wscatx(i,j,bi,bj) =
1592 & 1./wscatx(i,j,bi,bj)/wscatx(i,j,bi,bj)
1593 else
1594 wscatx(i,j,bi,bj) = 0.0 _d 0
1595 endif
1596 if (wtauv(i,j,bi,bj) .ne. 0.) then
1597 wtauv(i,j,bi,bj) = 1./wtauv(i,j,bi,bj)/wtauv(i,j,bi,bj)
1598 else
1599 wtauv(i,j,bi,bj) = 0.0 _d 0
1600 endif
1601 if (wtauvm(i,j,bi,bj) .ne. 0.) then
1602 wtauvm(i,j,bi,bj) =
1603 & 1./wtauvm(i,j,bi,bj)/wtauvm(i,j,bi,bj)
1604 else
1605 wtauvm(i,j,bi,bj) = 0.0 _d 0
1606 endif
1607 if (wscaty(i,j,bi,bj) .ne. 0.) then
1608 wscaty(i,j,bi,bj) =
1609 & 1./wscaty(i,j,bi,bj)/wscaty(i,j,bi,bj)
1610 else
1611 wscaty(i,j,bi,bj) = 0.0 _d 0
1612 endif
1613 if (whflux(i,j,bi,bj) .ne. 0.) then
1614 whflux(i,j,bi,bj) =
1615 & 1./whflux(i,j,bi,bj)/whflux(i,j,bi,bj)
1616 else
1617 whflux(i,j,bi,bj) = 0.0 _d 0
1618 endif
1619 if (whfluxm(i,j,bi,bj) .ne. 0.) then
1620 whfluxm(i,j,bi,bj) =
1621 & 1./whfluxm(i,j,bi,bj)/whfluxm(i,j,bi,bj)
1622 else
1623 whfluxm(i,j,bi,bj) = 0.0 _d 0
1624 endif
1625 if (wsflux(i,j,bi,bj) .ne. 0.) then
1626 wsflux(i,j,bi,bj) =
1627 & 1./wsflux(i,j,bi,bj)/wsflux(i,j,bi,bj)
1628 else
1629 wsflux(i,j,bi,bj) = 0.0 _d 0
1630 endif
1631 if (wsfluxm(i,j,bi,bj) .ne. 0.) then
1632 wsfluxm(i,j,bi,bj) =
1633 & 1./wsfluxm(i,j,bi,bj)/wsfluxm(i,j,bi,bj)
1634 else
1635 wsfluxm(i,j,bi,bj) = 0.0 _d 0
1636 endif
1637 if (wuwind(i,j,bi,bj) .ne. 0.) then
1638 wuwind(i,j,bi,bj) =
1639 & 1./wuwind(i,j,bi,bj)/wuwind(i,j,bi,bj)
1640 else
1641 wuwind(i,j,bi,bj) = 0.0 _d 0
1642 endif
1643 if (wvwind(i,j,bi,bj) .ne. 0.) then
1644 wvwind(i,j,bi,bj) =
1645 & 1./wvwind(i,j,bi,bj)/wvwind(i,j,bi,bj)
1646 else
1647 wvwind(i,j,bi,bj) = 0.0 _d 0
1648 endif
1649 if (watemp(i,j,bi,bj) .ne. 0.) then
1650 watemp(i,j,bi,bj) =
1651 & 1./watemp(i,j,bi,bj)/watemp(i,j,bi,bj)
1652 else
1653 watemp(i,j,bi,bj) = 0.0 _d 0
1654 endif
1655 if (waqh(i,j,bi,bj) .ne. 0.) then
1656 waqh(i,j,bi,bj) =
1657 & 1./waqh(i,j,bi,bj)/waqh(i,j,bi,bj)
1658 else
1659 waqh(i,j,bi,bj) = 0.0 _d 0
1660 endif
1661 if (wprecip(i,j,bi,bj) .ne. 0.) then
1662 wprecip(i,j,bi,bj) =
1663 & 1./wprecip(i,j,bi,bj)/wprecip(i,j,bi,bj)
1664 else
1665 wprecip(i,j,bi,bj) = 0.0 _d 0
1666 endif
1667 if (wswflux(i,j,bi,bj) .ne. 0.) then
1668 wswflux(i,j,bi,bj) =
1669 & 1./wswflux(i,j,bi,bj)/wswflux(i,j,bi,bj)
1670 else
1671 wswflux(i,j,bi,bj) = 0.0 _d 0
1672 endif
1673 if (wswdown(i,j,bi,bj) .ne. 0.) then
1674 wswdown(i,j,bi,bj) =
1675 & 1./wswdown(i,j,bi,bj)/wswdown(i,j,bi,bj)
1676 else
1677 wswdown(i,j,bi,bj) = 0.0 _d 0
1678 endif
1679 if (wlwflux(i,j,bi,bj) .ne. 0.) then
1680 wlwflux(i,j,bi,bj) =
1681 & 1./wlwflux(i,j,bi,bj)/wlwflux(i,j,bi,bj)
1682 else
1683 wlwflux(i,j,bi,bj) = 0.0 _d 0
1684 endif
1685 if (wlwdown(i,j,bi,bj) .ne. 0.) then
1686 wlwdown(i,j,bi,bj) =
1687 & 1./wlwdown(i,j,bi,bj)/wlwdown(i,j,bi,bj)
1688 else
1689 wlwdown(i,j,bi,bj) = 0.0 _d 0
1690 endif
1691 if (wevap(i,j,bi,bj) .ne. 0.) then
1692 wevap(i,j,bi,bj) =
1693 & 1./wevap(i,j,bi,bj)/wevap(i,j,bi,bj)
1694 else
1695 wevap(i,j,bi,bj) = 0.0 _d 0
1696 endif
1697 if (wsnowprecip(i,j,bi,bj) .ne. 0.) then
1698 wsnowprecip(i,j,bi,bj) =
1699 & 1./wsnowprecip(i,j,bi,bj)/wsnowprecip(i,j,bi,bj)
1700 else
1701 wsnowprecip(i,j,bi,bj) = 0.0 _d 0
1702 endif
1703 if (wapressure(i,j,bi,bj) .ne. 0.) then
1704 wapressure(i,j,bi,bj) =
1705 & 1./wapressure(i,j,bi,bj)/wapressure(i,j,bi,bj)
1706 else
1707 wapressure(i,j,bi,bj) = 0.0 _d 0
1708 endif
1709 if (wrunoff(i,j,bi,bj) .ne. 0.) then
1710 wrunoff(i,j,bi,bj) =
1711 & 1./wrunoff(i,j,bi,bj)/wrunoff(i,j,bi,bj)
1712 else
1713 wrunoff(i,j,bi,bj) = 0.0 _d 0
1714 endif
1715 if (wbottomdrag(i,j,bi,bj) .ne. 0.) then
1716 wbottomdrag(i,j,bi,bj) =
1717 & 1./wbottomdrag(i,j,bi,bj)/wbottomdrag(i,j,bi,bj)
1718 else
1719 wbottomdrag(i,j,bi,bj) = 0.0 _d 0
1720 endif
1721 if (wsfluxmm(bi,bj).ne.0.)
1722 & wsfluxmm(bi,bj) = 1./wsfluxmm(bi,bj)*wsfluxmm(bi,bj)
1723 if (whfluxmm(bi,bj).ne.0.)
1724 & whfluxmm(bi,bj) = 1./whfluxmm(bi,bj)*whfluxmm(bi,bj)
1725
1726 cph(
1727 if (whflux2(i,j,bi,bj) .ne. 0.) then
1728 whflux2(i,j,bi,bj) =
1729 & 1./whflux2(i,j,bi,bj)/whflux2(i,j,bi,bj)
1730 else
1731 whflux2(i,j,bi,bj) = 0.0 _d 0
1732 endif
1733 if (wsflux2(i,j,bi,bj) .ne. 0.) then
1734 wsflux2(i,j,bi,bj) =
1735 & 1./wsflux2(i,j,bi,bj)/wsflux2(i,j,bi,bj)
1736 else
1737 wsflux2(i,j,bi,bj) = 0.0 _d 0
1738 endif
1739 if (wtauu2(i,j,bi,bj) .ne. 0.) then
1740 wtauu2(i,j,bi,bj) =
1741 & 1./wtauu2(i,j,bi,bj)/wtauu2(i,j,bi,bj)
1742 else
1743 wtauu2(i,j,bi,bj) = 0.0 _d 0
1744 endif
1745 if (wtauv2(i,j,bi,bj) .ne. 0.) then
1746 wtauv2(i,j,bi,bj) =
1747 & 1./wtauv2(i,j,bi,bj)/wtauv2(i,j,bi,bj)
1748 else
1749 wtauv2(i,j,bi,bj) = 0.0 _d 0
1750 endif
1751 cph)
1752 enddo
1753 enddo
1754
1755 #ifdef ALLOW_OBCS_COST_CONTRIBUTION
1756 do iobcs = 1,nobcs
1757 do k = 1,nr
1758 #ifdef ALLOW_OBCSN_COST_CONTRIBUTION
1759 if (wobcsn(k,iobcs) .ne. 0.) then
1760 wobcsn(k,iobcs) =
1761 & ratio/wobcsn(k,iobcs)/wobcsn(k,iobcs)
1762 else
1763 wobcsn(k,iobcs) = 0.0 _d 0
1764 endif
1765 #endif
1766 #ifdef ALLOW_OBCSS_COST_CONTRIBUTION
1767 if (wobcss(k,iobcs) .ne. 0.) then
1768 wobcss(k,iobcs) =
1769 & ratio/wobcss(k,iobcs)/wobcss(k,iobcs)
1770 else
1771 wobcss(k,iobcs) = 0.0 _d 0
1772 endif
1773 #endif
1774 #ifdef ALLOW_OBCSW_COST_CONTRIBUTION
1775 if (wobcsw(k,iobcs) .ne. 0.) then
1776 wobcsw(k,iobcs) =
1777 & ratio/wobcsw(k,iobcs)/wobcsw(k,iobcs)
1778 else
1779 wobcsw(k,iobcs) = 0.0 _d 0
1780 endif
1781 #endif
1782 #ifdef ALLOW_OBCSE_COST_CONTRIBUTION
1783 if (wobcse(k,iobcs) .ne. 0.) then
1784 wobcse(k,iobcs) =
1785 & ratio/wobcse(k,iobcs)/wobcse(k,iobcs)
1786 else
1787 wobcse(k,iobcs) = 0.0 _d 0
1788 endif
1789 #endif
1790 enddo
1791 enddo
1792 #endif
1793
1794 enddo
1795 enddo
1796
1797 do bj = jtlo,jthi
1798 do bi = itlo,ithi
1799 do k = 1,nr
1800 if (wdiffkr(k,bi,bj) .ne. 0.) then
1801 wdiffkr(k,bi,bj) = 1./wdiffkr(k,bi,bj)/wdiffkr(k,bi,bj)
1802 else
1803 wdiffkr(k,bi,bj) = 0.0 _d 0
1804 endif
1805 if (wkapgm(k,bi,bj) .ne. 0.) then
1806 wkapgm(k,bi,bj) = 1./wkapgm(k,bi,bj)/wkapgm(k,bi,bj)
1807 else
1808 wkapgm(k,bi,bj) = 0.0 _d 0
1809 endif
1810 if (wkapredi(k,bi,bj) .ne. 0.) then
1811 wkapredi(k,bi,bj) = 1./wkapredi(k,bi,bj)/wkapredi(k,bi,bj)
1812 else
1813 wkapredi(k,bi,bj) = 0.0 _d 0
1814 endif
1815 if (wedtaux(k,bi,bj) .ne. 0.) then
1816 wedtaux(k,bi,bj) = 1./wedtaux(k,bi,bj)/wedtaux(k,bi,bj)
1817 else
1818 wedtaux(k,bi,bj) = 0.0 _d 0
1819 endif
1820 if (wedtauy(k,bi,bj) .ne. 0.) then
1821 wedtauy(k,bi,bj) = 1./wedtauy(k,bi,bj)/wedtauy(k,bi,bj)
1822 else
1823 wedtauy(k,bi,bj) = 0.0 _d 0
1824 endif
1825 do j = jmin,jmax
1826 do i = imin,imax
1827 if (wdiffkr2(i,j,k,bi,bj) .ne. 0.) then
1828 wdiffkr2(i,j,k,bi,bj) = frame(i,j)/
1829 & wdiffkr2(i,j,k,bi,bj)/wdiffkr2(i,j,k,bi,bj)
1830 else
1831 wdiffkr2(i,j,k,bi,bj) = 0.0 _d 0
1832 endif
1833 wdiffkrFld(i,j,k,bi,bj) = wdiffkr2(i,j,k,bi,bj)
1834 c
1835 if (wkapgm2(i,j,k,bi,bj) .ne. 0.) then
1836 wkapgm2(i,j,k,bi,bj) = frame(i,j)/
1837 & wkapgm2(i,j,k,bi,bj)/wkapgm2(i,j,k,bi,bj)
1838 else
1839 wkapgm2(i,j,k,bi,bj) = 0.0 _d 0
1840 endif
1841 wkapgmFld(i,j,k,bi,bj) = wkapgm2(i,j,k,bi,bj)
1842 c
1843 if (wkapredi2(i,j,k,bi,bj) .ne. 0.) then
1844 wkapredi2(i,j,k,bi,bj) = frame(i,j)/
1845 & wkapredi2(i,j,k,bi,bj)/wkapredi2(i,j,k,bi,bj)
1846 else
1847 wkapredi2(i,j,k,bi,bj) = 0.0 _d 0
1848 endif
1849 wkaprediFld(i,j,k,bi,bj) = wkapredi2(i,j,k,bi,bj)
1850 c
1851 if (wedtaux2(i,j,k,bi,bj) .ne. 0.) then
1852 wedtaux2(i,j,k,bi,bj) = frame(i,j)/
1853 & wedtaux2(i,j,k,bi,bj)/wedtaux2(i,j,k,bi,bj)
1854 else
1855 wedtaux2(i,j,k,bi,bj) = 0.0 _d 0
1856 endif
1857 wedtauxFld(i,j,k,bi,bj) = wedtaux2(i,j,k,bi,bj)
1858 c
1859 if (wedtauy2(i,j,k,bi,bj) .ne. 0.) then
1860 wedtauy2(i,j,k,bi,bj) = frame(i,j)/
1861 & wedtauy2(i,j,k,bi,bj)/wedtauy2(i,j,k,bi,bj)
1862 else
1863 wedtauy2(i,j,k,bi,bj) = 0.0 _d 0
1864 endif
1865 wedtauyFld(i,j,k,bi,bj) = wedtauy2(i,j,k,bi,bj)
1866 enddo
1867 enddo
1868 enddo
1869 enddo
1870 enddo
1871
1872 c
1873 c write masks and weights to files to be read by a master process
1874 c
1875 call active_write_xyz( 'maskCtrlC', maskC,
1876 & 1, 0, mythid, dummy)
1877 call active_write_xyz( 'maskCtrlW', maskW,
1878 & 1, 0, mythid, dummy)
1879 call active_write_xyz( 'maskCtrlS', maskS,
1880 & 1, 0, mythid, dummy)
1881
1882 #if (defined (ALLOW_HFLUX_COST_CONTRIBUTION) || defined (ALLOW_HFLUX_CONTROL))
1883 call active_write_xy( 'whflux', whflux, 1, 0, mythid, dummy)
1884 call active_write_xy( 'whflux2', whflux2, 1, 0, mythid, dummy)
1885 #elif (defined (ALLOW_ATEMP_COST_CONTRIBUTION) || defined (ALLOW_ATEMP_CONTROL))
1886 call active_write_xy( 'watemp', watemp, 1, 0, mythid, dummy)
1887 #endif
1888
1889 #if (defined (ALLOW_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SFLUX_CONTROL))
1890 call active_write_xy( 'wsflux', wsflux, 1, 0, mythid, dummy)
1891 call active_write_xy( 'wsflux2', wsflux2, 1, 0, mythid, dummy)
1892 #elif (defined (ALLOW_AQH_COST_CONTRIBUTION) || defined (ALLOW_AQH_CONTROL))
1893 call active_write_xy( 'waqh', waqh, 1, 0, mythid, dummy)
1894 #endif
1895
1896 #if (defined (ALLOW_PRECIP_COST_CONTRIBUTION) || defined (ALLOW_PRECIP_CONTROL))
1897 call active_write_xy( 'wprecip', wprecip, 1, 0, mythid, dummy)
1898 #endif
1899
1900 #if (defined (ALLOW_SWFLUX_COST_CONTRIBUTION) || defined (ALLOW_SWFLUX_CONTROL))
1901 call active_write_xy( 'wswflux', wswflux, 1, 0, mythid, dummy)
1902 #endif
1903
1904 #if (defined (ALLOW_SWDOWN_COST_CONTRIBUTION) || defined (ALLOW_SWDOWN_CONTROL))
1905 call active_write_xy( 'wswdown', wswdown, 1, 0, mythid, dummy)
1906 #endif
1907
1908 #if (defined (ALLOW_LWFLUX_COST_CONTRIBUTION) || defined (ALLOW_LWFLUX_CONTROL))
1909 call active_write_xy( 'wlwflux', wlwflux, 1, 0, mythid, dummy)
1910 #endif
1911
1912 #if (defined (ALLOW_LWDOWN_COST_CONTRIBUTION) || defined (ALLOW_LWDOWN_CONTROL))
1913 call active_write_xy( 'wlwdown', wlwdown, 1, 0, mythid, dummy)
1914 #endif
1915
1916 #if (defined (ALLOW_EVAP_COST_CONTRIBUTION) || defined (ALLOW_EVAP_CONTROL))
1917 call active_write_xy( 'wevap', wevap, 1, 0, mythid, dummy)
1918 #endif
1919
1920 #if (defined (ALLOW_SNOWPRECIP_COST_CONTRIBUTION) || defined (ALLOW_SNOWPRECIP_CONTROL))
1921 call active_write_xy( 'wsnowprecip', wsnowprecip,
1922 & 1, 0, mythid, dummy)
1923 #endif
1924
1925 #if (defined (ALLOW_APRESSURE_COST_CONTRIBUTION) || defined (ALLOW_APRESSURE_CONTROL))
1926 call active_write_xy( 'wapressure', wapressure,
1927 & 1, 0, mythid, dummy)
1928 #endif
1929
1930 #if (defined (ALLOW_RUNOFF_COST_CONTRIBUTION) || defined (ALLOW_RUNOFF_CONTROL))
1931 call active_write_xy( 'wrunoff', wrunoff, 1, 0, mythid, dummy)
1932 #endif
1933
1934 #if (defined (ALLOW_USTRESS_COST_CONTRIBUTION) || defined (ALLOW_USTRESS_CONTROL))
1935 call active_write_xy( 'wtauu', wtauu, 1, 0, mythid, dummy)
1936 call active_write_xy( 'wtauu2', wtauu2, 1, 0, mythid, dummy)
1937 #elif (defined (ALLOW_UWIND_COST_CONTRIBUTION) || defined (ALLOW_UWIND_CONTROL))
1938 call active_write_xy( 'wuwind', wuwind, 1, 0, mythid, dummy)
1939 #endif
1940
1941 #if (defined (ALLOW_VSTRESS_COST_CONTRIBUTION) || defined (ALLOW_VSTRESS_CONTROL))
1942 call active_write_xy( 'wtauv', wtauv, 1, 0, mythid, dummy)
1943 call active_write_xy( 'wtauv2', wtauv2, 1, 0, mythid, dummy)
1944 #elif (defined (ALLOW_VWIND_COST_CONTRIBUTION) || defined (ALLOW_VWIND_CONTROL))
1945 call active_write_xy( 'wvwind', wvwind, 1, 0, mythid, dummy)
1946 #endif
1947
1948 #ifdef ALLOW_OBCSN_COST_CONTRIBUTION
1949 #endif
1950
1951 #if (defined (ALLOW_KAPGM_COST_CONTRIBUTION) || defined (ALLOW_KAPGM_CONTROL))
1952 call active_write_xyz( 'wkapgmFld',wkapgmFld,
1953 & 1, 0, mythid, dummy)
1954 #endif
1955
1956 #if (defined (ALLOW_KAPREDI_COST_CONTRIBUTION) || defined (ALLOW_KAPREDI_CONTROL))
1957 call active_write_xyz( 'wkaprediFld',wkaprediFld,
1958 & 1, 0, mythid, dummy)
1959 #endif
1960
1961 #if (defined (ALLOW_DIFFKR_COST_CONTRIBUTION) || defined (ALLOW_DIFFKR_CONTROL))
1962 call active_write_xyz( 'wdiffkrFld',wdiffkrFld,
1963 & 1, 0, mythid, dummy)
1964 #endif
1965
1966 #if ( defined (ALLOW_EDDYPSI_COST_CONTRIBUTION) || defined (ALLOW_EDDYPSI_CONTROL) )
1967 call active_write_xyz( 'wedtauxFld',wedtauxFld,
1968 & 1, 0, mythid, dummy)
1969 call active_write_xyz( 'wedtauyFld',wedtauyFld,
1970 & 1, 0, mythid, dummy)
1971 #endif
1972
1973 #if (defined (ALLOW_BOTTOMDRAG_COST_CONTRIBUTION) || defined (ALLOW_BOTTOMDRAG_CONTROL))
1974 call active_write_xy( 'wbottomdrag', wbottomdrag
1975 & , 1, 0, mythid, dummy)
1976 #endif
1977 CMM(
1978 #ifdef ALLOW_EGM96_ERROR_DIAG
1979 call active_write_xy( 'wp_DOT', wp
1980 & , 1, 0, mythid, dummy)
1981 #endif
1982 CMM)
1983 end
  ViewVC Help
Powered by ViewVC 1.1.22