/[MITgcm]/MITgcm_contrib/verification_other/global1x1_tot/code/ecco_cost_weights.F
ViewVC logotype

Annotation of /MITgcm_contrib/verification_other/global1x1_tot/code/ecco_cost_weights.F

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


Revision 1.1 - (hide annotations) (download)
Sat Feb 4 02:45:01 2012 UTC (13 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64k, checkpoint64x, checkpoint64z, checkpoint64o, checkpoint64p, checkpoint64r, checkpoint64w, checkpoint64v, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint64m, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65u, checkpoint65j, checkpoint67a, checkpoint67b, checkpoint65i, checkpoint67d, checkpoint65m, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65e, checkpoint64i, checkpoint64h, checkpoint65h, checkpoint65, checkpoint64j, checkpoint65n, HEAD
move un-tested set-up from MITgcm/verification/global1x1_tot to Contrib/verification_other

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

  ViewVC Help
Powered by ViewVC 1.1.22