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

Annotation of /MITgcm_contrib/arctic40km/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)
Wed May 3 23:59:39 2006 UTC (19 years, 2 months ago) by heimbach
Branch: MAIN
Adding state estimation setup for arctic50km/

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

  ViewVC Help
Powered by ViewVC 1.1.22