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

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

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


Revision 1.1 - (hide annotations) (download)
Fri Apr 23 19:55:12 2010 UTC (15 years, 3 months ago) by mmazloff
Branch: MAIN
CVS Tags: HEAD
original files

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

  ViewVC Help
Powered by ViewVC 1.1.22