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

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

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

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