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